Reaction and ODE updates

Author: cwlanyon

advectionGP/models.py

@@ -323,10 +323,8 @@ def computeConcentration(self,source,enforce_nonnegative=False):
 
         self.conc = concentration
         return c      
-   
-
 
-class AdjointAdvectionDiffusionReactionModel(AdjointAdvectionDiffusionModel,AdvectionDiffusionReactionModel):
+class AdjointAdvectionDiffusionReactionModel(AdvectionDiffusionReactionModel,AdjointAdvectionDiffusionModel):
     def computeAdjoint(self,H):
         """
         Runs the backward PDE (adjoint problem)
@@ -359,6 +357,69 @@ def computeAdjoint(self,H):
             v[-i-1,1:Nx-1,1:Ny-1]=v[-i,1:Nx-1,1:Ny-1] +dt*( H[-i,1:Nx-1,1:Ny-1]+u[0][-i,1:Nx-1,1:Ny-1]*(v[-i,2:Nx,1:Ny-1]-v[-i,0:Nx-2,1:Ny-1])/(2*dx) +u[1][-i,1:Nx-1,1:Ny-1]*(v[-i,1:Nx-1,2:Ny]-v[-i,1:Nx-1,0:Ny-2] )/(2*dy)+k_0*(v[-i,2:Nx,1:Ny-1]-2*v[-i,1:Nx-1,1:Ny-1]  +v[-i,0:Nx-2,1:Ny-1])/dx2+k_0*(v[-i,1:Nx-1,2:Ny]-2*v[-i,1:Nx-1,1:Ny-1]  +v[-i,1:Nx-1,0:Ny-2])/dy2 -R*v[-i,1:Nx-1,1:Ny-1])
         return v
 
+class SimpleODEModel(AdvectionDiffusionModel):
+    def __init__(self,boundary,resolution,kernel,noiseSD,sensormodel,N_feat=25,spatial_averaging=1.0):
+        """
+        The Advection Diffusion Reaction Model.
+        
+        At the moment we assume a 3d grid [time, x, y].
+        
+        Parameters:
+            boundary = a two element tuple of the corners of the grid. e.g. ([0,0,0],[10,10,10])        
+            resolution = a list of the grid size in each dimension. e.g. [10,20,20]
+            kernel = the kernel to use
+            noiseSD = the noise standard deviation
+            sensormodel = an instatiation of a SensorModel class that implements the getHs method.
+            N_feat = number of fourier features
+            spatial_averaging = how big the volume the sensor measures (default 0.001).
+            u = wind speed
+            k_0 = diffusion constant
+            R = reaction constant
+        """
+        super().__init__(boundary,resolution,kernel,noiseSD,sensormodel,N_feat,spatial_averaging)        
+        
+                
+    def computeConcentration(self,source,enforce_nonnegative=False):
+        """
+        Computes concentrations.
+        Arguments:
+         source == forcing function (shape: Nt x Nx x Ny). Can either be generated by ... or determine manually.
+         enforce_nonnegative = default False,. Setting to true will force concentration to be non-negative each iteration.
+        returns array of concentrations (shape: Nt x Nx x Ny), given source. (also saved it in self.concentration)
+        """
+        #source = self.source
+        
+        #get the grid step sizes, their squares and the size of the grid
+        dt,dx,dy,dx2,dy2,Nt,Nx,Ny = self.getGridStepSize()
+        
+        c=np.zeros(((self.resolution)))
+        
+        c[0,:,:]=0
+
+        for i in range(0,Nt-1):
+            # Internal Calc
+            c[i+1,:,:]=c[i,:,:] +dt*(source[i,:,:])
+            if enforce_nonnegative: c[c<0]=0
+        concentration = c 
+        
+        self.conc = concentration
+        return c      
+
+class AdjointSimpleODEModel(SimpleODEModel,AdjointAdvectionDiffusionModel):
+    def computeAdjoint(self,H):
+        """
+        Runs the backward PDE (adjoint problem)
+        Gets called for an observation instance (H).
+        (v is the result of the adjoint operation)
+        """
+        dt,dx,dy,dx2,dy2,Nt,Nx,Ny = self.getGridStepSize()
+
+        v=np.zeros(((Nt,Nx,Ny)))
+        v[-1,:,:]=0.0
+        for i in range(1,Nt): #TODO might be better to rewrite as range(Nt-1,1,-1)...
+    #Internal calculation (not on the boundary)
+            v[-i-1,:,:]=v[-i,:,:] +dt*( H[-i,:,:])
+        return v
 
 class MCMCAdvectionDiffusionModel(AdvectionDiffusionModel):
     def computeLikelihood(self,pred_y,act_y):