Merge pull request #13 from OpenDA-Association/12-integrating-neural-…

…networks

12 adding neural networks

openda/algorithms/PINN.py

@@ -0,0 +1,196 @@
+import numpy as np
+import torch
+from torch import nn
+
+
+class NN(nn.Module):
+    """
+    Interface of the (ordinary) Neural Network.
+
+    Input (x_train) is a tensor containing states [state, previous_state].
+    Ouptut (y_train) is a tensor containing estimated parameter(s) [f].
+
+    @author Aron Schouten
+
+    """
+    def __init__(self, device, layers, enkf, data):
+        """
+        Constructor.
+        """
+        super().__init__()
+        self.D = enkf.model_factory.model_attributes[0]['D']
+        self.g = enkf.model_factory.model_attributes[0]['g']
+        self.dx = enkf.model_factory.model_attributes[0]['L']/(enkf.model_factory.model_attributes[0]['n']+0.5)
+        self.dt = enkf.model_factory.model_attributes[4][1]/np.timedelta64(1,'s')
+
+        self.depth = len(layers)-2
+
+        self.y_max = max( [data[1].max(), data[3].max()] )
+        self.y_min = min( [data[1].min(), data[3].min()] )
+
+        self.x_train = data[0].to(device)
+        self.y_train = (data[1].to(device) - self.y_min) / (self.y_max - self.y_min)
+
+        self.x_test = data[2].to(device)
+        self.y_test = (data[3].to(device) - self.y_min) / (self.y_max - self.y_min)
+
+        # self.u_b = max([self.x_train.max(), self.x_test.max()])
+        # self.l_b = min([self.x_train.min(), self.x_test.min()])
+
+        self.activation = nn.ReLU()
+        self.linears = nn.ModuleList([nn.Linear(layers[i], layers[i+1]) for i in range(self.depth+1)])
+
+        self.device = device
+
+        for i in range(self.depth+1):
+            nn.init.kaiming_normal_(self.linears[i].weight, nonlinearity='relu')
+            nn.init.constant_(self.linears[i].bias, 0.01)
+
+    def forward(self, x):
+        """
+        Forward x through the neural network.
+
+        :param x: States.
+        :return: Forwarded states.
+        """
+        # x = (x - self.l_b)/(self.u_b - self.l_b)
+
+        for i in range(self.depth):
+            x = self.linears[i](x)
+            x = self.activation(x)
+        x = self.linears[-1](x)
+
+        return x
+
+    def loss(self, x, y):
+        """
+        Calculates loss of estimated parameter(s).
+
+        :param x: States.
+        :param y: True value of parameter(s).
+        :return: MSE of estimated parameter(s), 1/N*∑|f_est - f_real|^2.
+        """
+        loss = ( (self.forward(x) - y)**2 ).mean()
+        return loss
+
+    def train_model(self, optimizer, n_epochs, batch_size):
+        """
+        Train the model.
+
+        :param optimizer: The (torch.)optimizer that is used for training.
+        :param n_epochs: Number of epochs.
+        :param batch_size: The batch size to train with.
+        :return: Lists with epoch numbers, the in-sample MSE's, out-sample MSE's, and (out-sample) biases.
+        """
+        output = [[], [], []]
+        for epoch in range(n_epochs):
+            for i in range(0, len(self.x_train), batch_size):
+                batch_X = self.x_train[i:i+batch_size]
+                batch_y = self.y_train[i:i+batch_size]
+
+                optimizer.zero_grad()
+                loss = self.loss(batch_X, batch_y)
+                loss.backward()
+                # print(loss)
+                optimizer.step()
+
+            val_loss = self.test_model()
+            output[0].append(epoch)
+            output[1].append(loss.item())
+            output[2].append(val_loss.item())
+            print(f'Epoch {epoch}: Loss = {loss}, Validaton loss = {val_loss}')
+
+        return output
+
+    def test_model(self):
+        """
+        Temporary puts the model in evaluation mode and calculates validation loss.
+
+        :return: Validation loss (loss on test data).
+        """
+        self.eval() # Necessary when using Dropout or BatchNorm layers
+        with torch.no_grad(): # To make sure we don't train when testing
+            val_loss = self.loss(self.x_test, self.y_test)
+        self.train()
+
+        return val_loss
+
+    def predict(self, x):
+        """
+        Temporary puts the model in evaluation mode and calculates estimates parameters.
+
+        :return: Estimated parameters [f]
+        """
+        self.eval() # Necessary when using Dropout or BatchNorm layers
+        with torch.no_grad(): # To make sure we don't train when testing
+            y_pred = self(x)
+        self.train()
+
+        return y_pred*self.y_max + (1-y_pred)*self.y_min
+
+
+class PINN(NN):
+    """
+    Interface of the Physics-Informed Neural Network.
+
+    Input (x_train) is a tensor containing states [state, previous_state].
+    Ouptut (y_train) is a tensor containing estimated parameter(s) [f].
+
+    @author Aron Schouten
+
+    """
+
+    def loss_param(self, x, y):
+        """
+        Calculates loss of estimated parameter(s).
+
+        :param x: States.
+        :param y: True value of parameter(s).
+        :return: MSE of estimated parameter(s), 1/N*∑|f_est - f_real|^2.
+        """
+        loss = ( (self.forward(x)- y)**2 ).mean()
+        return loss
+
+    def loss_PDE(self, x):
+        """
+        Calculates loss of PDE's.
+
+        :param x: States.
+        :return: MSE of PDE's, 1/N*∑|h_t + g * u_x|^2 + 1/N*∑|u_t + D * h_x + f_est * u|^2.
+        """
+        f = self.y_min + self.forward(x)*(self.y_max - self.y_min)
+
+        n = x.shape[1]//2
+
+        h = x[:,:n][:,0:-1:2]
+        u = x[:,:n][:,1::2]
+        prev_h = x[:,n:][:,0:-1:2]
+        prev_u = x[:,n:][:,1::2]
+
+        h_t = (h - prev_h)/self.dt
+        u_t = (u - prev_u)/self.dt
+
+        # Extending u and h to keep the same dimensions after taking derivative
+        h_ext = torch.concat((torch.zeros(h.shape[0], device=self.device).view(-1,1), h),-1)
+        u_ext = torch.concat((torch.zeros(u.shape[0], device=self.device).view(-1,1), u),-1)
+
+        h_x = (h_ext[:,1:] - h_ext[:,:-1])/self.dx
+        u_x = (u_ext[:,1:] - u_ext[:,:-1])/self.dx
+
+        loss = ( (h_t + self.g*u_x)**2 ).mean()
+        loss += ( (u_t + self.D*h_x + f*u)**2 ).mean()
+
+        return loss
+
+    def loss(self, x, y):
+        """
+        Calculates total loss.
+
+        :param x: States.
+        :param y: True value of parameter(s).
+        :return: Sum of all losses.
+        """
+        loss_param = self.loss_param(x, y)
+        loss_PDE = self.loss_PDE(x)
+        # print(f'Loss for |f-f*|^2 is {loss_param}. Loss for |PDE|^2 is {loss_PDE}. Ratio = {loss_param/loss_PDE}')
+        return loss_param + loss_PDE

openda/models/SaintVenantStochModelFactory.py

@@ -9,12 +9,15 @@ class SaintVenantModelFactory:
     @author Aron Schouten
     """
 
-    def __init__(self):
+    def __init__(self, f=None):
         """
         Constructor
         """
+        if not f: f = np.random.uniform(1e-5, 1e-3)
+        print(f'Initial f = {f}')
+
         names = ["D", "f", "g", "L", "n"]
-        param_values = [15, 2.3e-4, 9.81, 55.e3, 11]
+        param_values = [15, f, 9.81, 55.e3, 11]
         param_uncertainty = [0, 0, 0, 0, 0]
 
         param = dict(zip(names, param_values))

openda/models/SaintVenantStochModelInstance.py

@@ -42,7 +42,9 @@ def __init__(self, model_attributes, noise_config, main_or_ens=None):
 
         self.current_time = PyTime(self.span[0])
         self.state = np.array(self.state)
+        self.prev_state = np.zeros_like(self.state)
         self.t = 0
+        self.f = self.param['f']
 
     def get_time_horizon(self):
         """
@@ -80,8 +82,9 @@ def compute(self, time):
         :param time: Time to compute to.
         :return:
         """
+        self.prev_state = self.state.copy()
         end_time = time.get_start()
-        A, B, phi = _get_model(self.param, self.span)
+        A, B, phi = self._get_model()
         std = np.sqrt(1-phi**2) * 0.2 # Std of model noise chosen according to desired std of AR(1)
         newx = self.state
         t_now = self.current_time.get_start()
@@ -99,6 +102,9 @@ def compute(self, time):
         self.current_time = PyTime(end_time)
         self.state = newx
 
+    def set_f(self, f):
+        self.f = f
+
     def get_observations(self, description):
         """
         Get model values corresponding to the descriptions.
@@ -140,50 +146,57 @@ def get_state(self):
         :return: State vector.
         """
         return self.state
-
-
-def _get_model(param, span):
-    """
-    Returns model matrices A and B such that A*x_new=B*x
-    A and B are tri-diagonal sparce matrices, and have the order h[0], u[0], ..., h[n], u[n]  
-    """
-    n=param['n']
-    Adata=np.zeros((3,2*n+1))
-    Bdata=np.zeros((3,2*n+1))
-    Adata[1,0]=1. # Left boundary
-    Adata[1,2*n-1]=1. # Right boundary
-    Adata[1,2*n] = 1
-    phi = np.exp( -(span[1]/np.timedelta64(1,'s'))/(6.*60.*60.) )
-    Bdata[1,2*n] = phi
-    # i=1,3,5,... du/dt  + g dh/sx + f u = 0
-    #  u[n+1,m] + 0.5 g dt/dx ( h[n+1,m+1/2] - h[n+1,m-1/2]) + 0.5 dt f u[n+1,m] 
-    #= u[n  ,m] - 0.5 g dt/dx ( h[n  ,m+1/2] - h[n  ,m-1/2]) - 0.5 dt f u[n  ,m]
-    dt = span[1]/np.timedelta64(1,'s')
-    dx = param['L']/(n+0.5)
-    temp1=0.5*param['g']*dt/dx
-    temp2=0.5*param['f']*dt
-    for i in np.arange(1,2*n-1,2):
-        Adata[0,i-1]= -temp1
-        Adata[1,i  ]= 1.0 + temp2
-        Adata[2,i+1]= +temp1
-        Bdata[0,i-1]= +temp1
-        Bdata[1,i  ]= 1.0 - temp2
-        Bdata[2,i+1]= -temp1
-    # i=2,4,6,... dh/dt + D du/dx = 0
-    #  h[n+1,m] + 0.5 D dt/dx ( u[n+1,m+1/2] - u[n+1,m-1/2])  
-    #= h[n  ,m] - 0.5 D dt/dx ( u[n  ,m+1/2] - u[n  ,m-1/2])
-    temp1=0.5*param['D']*dt/dx
-    for i in np.arange(2,2*n,2):
-        Adata[0,i-1]= -temp1
-        Adata[1,i  ]= 1.0
-        Adata[2,i+1]= +temp1
-        Bdata[0,i-1]= +temp1
-        Bdata[1,i  ]= 1.0
-        Bdata[2,i+1]= -temp1    
-    # Build sparse matrix
-    A=spdiags(Adata,np.array([-1,0,1]),2*n+1,2*n+1).tolil()
-    A[0, -1]=-1
-
-    B=spdiags(Bdata,np.array([-1,0,1]),2*n+1,2*n+1)
 
-    return A.tocsr(), B.tocsr(), phi
+    def get_prev_state(self):
+        """
+        Returns the previous state of the model.
+
+        :return: Previous state vector.
+        """
+        return self.prev_state
+
+    def _get_model(self):
+        """
+        Returns model matrices A and B such that A*x_new=B*x
+        A and B are tri-diagonal sparce matrices, and have the order h[0], u[0], ..., h[n], u[n]  
+        """
+        n=self.param['n']
+        Adata=np.zeros((3,2*n+1))
+        Bdata=np.zeros((3,2*n+1))
+        Adata[1,0]=1. # Left boundary
+        Adata[1,2*n-1]=1. # Right boundary
+        Adata[1,2*n] = 1
+        phi = np.exp( -(self.span[1]/np.timedelta64(1,'s'))/(6.*60.*60.) )
+        Bdata[1,2*n] = phi
+        # i=1,3,5,... du/dt  + g dh/sx + f u = 0
+        #  u[n+1,m] + 0.5 g dt/dx ( h[n+1,m+1/2] - h[n+1,m-1/2]) + 0.5 dt f u[n+1,m] 
+        #= u[n  ,m] - 0.5 g dt/dx ( h[n  ,m+1/2] - h[n  ,m-1/2]) - 0.5 dt f u[n  ,m]
+        dt = self.span[1]/np.timedelta64(1,'s')
+        dx = self.param['L']/(n+0.5)
+        temp1=0.5*self.param['g']*dt/dx
+        temp2=0.5*self.f*dt
+        for i in np.arange(1,2*n-1,2):
+            Adata[0,i-1]= -temp1
+            Adata[1,i  ]= 1.0 + temp2
+            Adata[2,i+1]= +temp1
+            Bdata[0,i-1]= +temp1
+            Bdata[1,i  ]= 1.0 - temp2
+            Bdata[2,i+1]= -temp1
+        # i=2,4,6,... dh/dt + D du/dx = 0
+        #  h[n+1,m] + 0.5 D dt/dx ( u[n+1,m+1/2] - u[n+1,m-1/2])  
+        #= h[n  ,m] - 0.5 D dt/dx ( u[n  ,m+1/2] - u[n  ,m-1/2])
+        temp1=0.5*self.param['D']*dt/dx
+        for i in np.arange(2,2*n,2):
+            Adata[0,i-1]= -temp1
+            Adata[1,i  ]= 1.0
+            Adata[2,i+1]= +temp1
+            Bdata[0,i-1]= +temp1
+            Bdata[1,i  ]= 1.0
+            Bdata[2,i+1]= -temp1    
+        # Build sparse matrix
+        A=spdiags(Adata,np.array([-1,0,1]),2*n+1,2*n+1).tolil()
+        A[0, -1]=-1
+
+        B=spdiags(Bdata,np.array([-1,0,1]),2*n+1,2*n+1)
+
+        return A.tocsr(), B.tocsr(), phi

setup.cfg

@@ -39,6 +39,7 @@ install_requires =
     scipy
     pandas
     matplotlib
+    torch
     xmlschema
 
 [options.data_files]

observations/generate_observations.py → tests/observations/generate_observations.py

@@ -1,9 +1,10 @@
-import matplotlib.pyplot as plt
 import matplotlib.dates as mdates
-import pandas as pd
+import matplotlib.pyplot as plt
 import numpy as np
-from openda.models.SaintVenantStochModelFactory import SaintVenantModelFactory
+import pandas as pd
 from openda.costFunctions.JObjects import PyTime
+from openda.models.SaintVenantStochModelFactory import SaintVenantModelFactory
+
 
 def plot_series(t, series_data, xlocs_waterlevel, xlocs_velocity):
     titles=[]
@@ -63,7 +64,7 @@ def test():
     for i in description:
         header.append(i)
     df.columns = header
-    df.to_csv(r".\observations\obs (simulated).csv", sep=';', index=False)
+    df.to_csv(r"./tests/observations/obs_simulated_5min.csv", sep=';', index=False)
 
 
 if __name__ == '__main__':