Modifications using numpy operations, CI

.github/workflows/CI_testing.yml

@@ -0,0 +1,36 @@
+name: CI
+on:
+  push:
+    branches:
+      - numpy_debug_v2
+
+defaults:
+  run:
+    shell: bash -l {0}
+jobs:
+  style:
+    runs-on: ubuntu-22.04
+    container:
+      image: dolfinx/dolfinx:stable
+      options: --user 1001 --privileged 
+    name: CI test
+    steps:
+      - name: run code
+        uses: actions/checkout@v2
+      - name: run serial check
+        run: |
+          cd ./hippylibX/test && 
+          output=$(mpirun -n 1 python3 testing_suite_file.py) &&
+          if [[ "$output" == "1"  ]]; then
+            exit 1
+          fi
+
+      - name: run parallel check
+        run: |
+          cd ./hippylibX/test && 
+          output=$(mpirun -n 2 python3 testing_suite_file.py) &&
+          if [[ "$output" == "1"  ]]; then
+            exit 1
+          fi
+
+

Dockerfile

@@ -1,3 +1,4 @@
 FROM  dolfinx/dolfinx:stable
 MAINTAINER Venugopal Ranganathan <[email protected]>
+RUN pip install memory_profiler
 RUN pip install psutil

example/poisson_example.py

@@ -8,11 +8,9 @@
 import os
 import dolfinx.fem.petsc
 
-
 from matplotlib import pyplot as plt
 
 sys.path.append( os.environ.get('HIPPYLIBX_BASE_DIR', "../") )
-
 import hippylibX as hpx
 
 def master_print(comm, *args, **kwargs):
@@ -94,7 +92,7 @@ def run_inversion(nx : int, ny : int, noise_variance : float, prior_param : dict
     d.x.scatter_forward()
 
     misfit_form = PoissonMisfitForm(d,noise_variance)
-    misfit = hpx.NonGaussianContinuousMisfit(msh, Vh, misfit_form)
+    misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form)
 
     prior_mean = dlx.fem.Function(Vh_m)
     prior_mean.x.array[:] = 0.01
@@ -108,17 +106,16 @@ def run_inversion(nx : int, ny : int, noise_variance : float, prior_param : dict
     hpx.parRandom.normal(1.,noise)
     prior.sample(noise,m0)
 
-    hpx.modelVerify(comm,model,m0,is_quadratic=False,misfit_only=True,verbose=(rank == 0))
-    hpx.modelVerify(comm,model,m0,is_quadratic=False,misfit_only=False,verbose=(rank == 0))
 
+    eps, err_grad, err_H,rel_symm_error = hpx.modelVerify(model,m0,is_quadratic=False,misfit_only=True,verbose=(rank == 0))
 
-    # eps, err_grad, err_H = hpx.modelVerify(comm,model,m0,is_quadratic=False,misfit_only=True,verbose=(rank == 0))
+    data_misfit_True = {"eps":eps,"err_grad":err_grad, "err_H": err_H, "sym_Hessian_value":rel_symm_error}
 
-    # if(rank == 0):
-    #     print(err_grad,'\n')    
-    #     print(err_H)
-    #     plt.show()  
 
+    eps, err_grad, err_H,rel_symm_error = hpx.modelVerify(model,m0,is_quadratic=False,misfit_only=False,verbose=(rank == 0))
+
+    data_misfit_False = {"eps":eps,"err_grad":err_grad, "err_H": err_H, "sym_Hessian_value":rel_symm_error}
+
     # # #######################################
 
     prior_mean_copy = prior.generate_parameter(0)
@@ -151,18 +148,30 @@ def run_inversion(nx : int, ny : int, noise_variance : float, prior_param : dict
     master_print (comm, "Termination reason: ", solver.termination_reasons[solver.reason])
     master_print (comm, "Final gradient norm: ", solver.final_grad_norm)
     master_print (comm, "Final cost: ", solver.final_cost)
+
+    optimizer_results = {}
+    if(solver.termination_reasons[solver.reason] == 'Norm of the gradient less than tolerance'):
+        optimizer_results['optimizer']  = True
+    else:
+        optimizer_results['optimizer'] = False
 
-    #######################################
 
+    final_results = {"data_misfit_True":data_misfit_True,
+                     "data_misfit_False":data_misfit_False,
+                     "optimizer_results":optimizer_results}
 
 
+    return final_results
+
+
+    #######################################
 
 if __name__ == "__main__":    
     nx = 64
     ny = 64
     noise_variance = 1e-4
     prior_param = {"gamma": 0.1, "delta": 1.}
     run_inversion(nx, ny, noise_variance, prior_param)
-
-
+    plt.savefig("poisson_result_FD_Gradient_Hessian_Check")
+    plt.show()
 

example/sfsi_toy_gaussian.py

@@ -8,7 +8,6 @@
 import os
 import dolfinx.fem.petsc
 
-
 from matplotlib import pyplot as plt
 
 sys.path.append( os.environ.get('HIPPYLIBX_BASE_DIR', "../") )
@@ -61,7 +60,6 @@ def run_inversion(nx : int, ny : int, noise_variance : float, prior_param : dict
     comm = MPI.COMM_WORLD
     rank  = comm.rank
     nproc = comm.size
-
     fname = 'meshes/circle.xdmf'
     fid = dlx.io.XDMFFile(comm,fname,"r")
     msh = fid.read_mesh(name='mesh')
@@ -110,7 +108,7 @@ def run_inversion(nx : int, ny : int, noise_variance : float, prior_param : dict
     d.x.scatter_forward()
 
     misfit_form = PACTMisfitForm(d, noise_variance)
-    misfit = hpx.NonGaussianContinuousMisfit(msh, Vh, misfit_form)
+    misfit = hpx.NonGaussianContinuousMisfit(Vh, misfit_form)
 
     prior_mean = dlx.fem.Function(Vh_m)
     prior_mean.x.array[:] = 0.01
@@ -124,14 +122,15 @@ def run_inversion(nx : int, ny : int, noise_variance : float, prior_param : dict
     hpx.parRandom.normal(1.,noise)
     prior.sample(noise,m0)
 
-    hpx.modelVerify(comm,model,m0,is_quadratic=False,misfit_only=True,verbose=(rank == 0))
-    hpx.modelVerify(comm,model,m0,is_quadratic=False,misfit_only=False,verbose=(rank == 0))
 
-    # if(rank == 0):
-    #     print(err_grad,'\n')
-    #     print(err_H)
-    #     plt.show()  
+    eps, err_grad, err_H,rel_symm_error = hpx.modelVerify(model,m0,is_quadratic=False,misfit_only=True,verbose=(rank == 0))
 
+    data_misfit_True = {"eps":eps,"err_grad":err_grad, "err_H": err_H, "sym_Hessian_value":rel_symm_error}
+
+    eps, err_grad, err_H,rel_symm_error = hpx.modelVerify(model,m0,is_quadratic=False,misfit_only=False,verbose=(rank == 0))
+
+    data_misfit_False = {"eps":eps,"err_grad":err_grad, "err_H": err_H, "sym_Hessian_value":rel_symm_error}
+
     # #######################################
 
     prior_mean_copy = prior.generate_parameter(0)
@@ -166,6 +165,21 @@ def run_inversion(nx : int, ny : int, noise_variance : float, prior_param : dict
     master_print (comm, "Final gradient norm: ", solver.final_grad_norm)
     master_print (comm, "Final cost: ", solver.final_cost)
 
+
+    optimizer_results = {}
+    if(solver.termination_reasons[solver.reason] == 'Norm of the gradient less than tolerance'):
+        optimizer_results['optimizer']  = True
+    else:
+        optimizer_results['optimizer'] = False
+
+
+    final_results = {"data_misfit_True":data_misfit_True,
+                     "data_misfit_False":data_misfit_False,
+                     "optimizer_results":optimizer_results}
+
+    return final_results
+
+
     #######################################
 
 if __name__ == "__main__":    
@@ -174,5 +188,7 @@ def run_inversion(nx : int, ny : int, noise_variance : float, prior_param : dict
     noise_variance = 1e-6
     prior_param = {"gamma": 0.05, "delta": 1.}
     run_inversion(nx, ny, noise_variance, prior_param)
+    plt.savefig("qpact_result_FD_Gradient_Hessian_Check")
+    plt.show()
 
 

example/testing_folder/test_script.py

@@ -0,0 +1,36 @@
+import numpy as np
+import pickle
+
+with open('outputs.pickle', 'rb') as f:
+    data = pickle.load(f)
+
+xvals = data["xvals"]
+arr_1 = data['arr_1']
+arr_2 = data['arr_2']
+value = data['sym_Hessian_value']
+
+if(value > 1e-10):
+    print("1")
+    exit
+
+#values from arr_1
+relevant_values_1, relevant_values_2 = arr_1[10:], arr_2[10:]
+
+x =  xvals[10:]
+
+data = relevant_values_1
+y = data
+coefficients = np.polyfit(np.log(x), np.log(y), 1)
+slope_1 = coefficients[0]
+
+data = relevant_values_2
+y = data
+coefficients = np.polyfit(np.log(x), np.log(y), 1)
+slope_2 = coefficients[0]
+
+
+if(np.abs(slope_1 - 1) > 1e-1 or np.abs(slope_2 - 1) > 1e-1):
+    print("1")
+    exit
+
+print("0")

hippylibX/algorithms/NewtonCG.py

@@ -200,7 +200,6 @@ def _solve_ls(self,x):
                 gradnorm_ini = gradnorm
                 tol = max(abs_tol, gradnorm_ini*rel_tol)
 
-            # print(tol)
 
             # check if solution is reached
             if (gradnorm < tol) and (self.it > 0):

hippylibX/algorithms/linSolvers.py

@@ -2,7 +2,7 @@
 
 
 #I don't use either of the following 2 functions for now,
-# currently just using this in the _createLUSolver method in
+# currently just using the following: in the _createLUSolver method in
 # the PDEVariationalProblem class in modeling/PDEProblem.py file
 
 # def _createLUSolver(self):

hippylibX/modeling/PDEProblem.py

@@ -238,9 +238,19 @@ def setLinearizationPoint(self,x : list, gauss_newton_approx) -> None:
         self.solver_adj_inc.setOperators(self.At)
 
         if gauss_newton_approx:
+            if(self.Wuu is not None):
+                self.Wuu.destroy()
+
             self.Wuu = None
+
+            if(self.Wmu is not None):
+                self.Wmu.destroy()
             self.Wmu = None
+
+            if(self.Wmm is not None):
+                self.Wmm.destroy()
             self.Wmm = None
+
         else:
 
             if self.Wuu is None:
@@ -265,8 +275,6 @@ def setLinearizationPoint(self,x : list, gauss_newton_approx) -> None:
             self.Wmu = self.Wum.copy()
             self.Wmu.transpose()
 
-            # print(self.Wum.isTranspose(self.Wmu)) #True
-
             if self.Wmm is None:
                 self.Wmm = dlx.fem.petsc.create_matrix(dlx.fem.form(ufl.derivative(g_form[PARAMETER],x_fun[PARAMETER],x_fun_trial[PARAMETER])))
 

hippylibX/modeling/Regularization.py

@@ -8,69 +8,38 @@
 import numpy as np
 
 class VariationalRegularization:
-    def __init__(self, mesh, Vh, functional_handler, isQuadratic=False):
+    def __init__(self, Vh, functional_handler, isQuadratic=False):
         self.Vh = Vh #Function space of the parameter.
         self.functional_handler = functional_handler #a function or a functor that takes as input m (as dl.Function) and evaluate the regularization functional
         self.isQuadratic = isQuadratic #Whether the functional is a quadratic form (i.e. the Hessian is constant) or not (the Hessian dependes on m
-        self.mesh = mesh #to allreduce over the entire mesh in the cost function
 
         self.xfun = dlx.fem.Function(self.Vh)
-
-    # def cost(self, m):
-    # 1. Cast the petsc4py vector m to a dlx.Function mfun
-    # 2. Compute the cost by calling assemble on self.functional_handler(mfun)
-    # 3. Return the value of the cost (Make sure to call a AllReduce for parallel computations
-
-    def cost(self,m):
-        # mfun = vector2Function(m,self.Vh)
 
+    def cost(self,m):    
         updateFromVector(self.xfun, m)
         mfun = self.xfun
-
         loc_cost = self.functional_handler(mfun)
         glb_cost_proc = dlx.fem.assemble_scalar(dlx.fem.form(loc_cost))
-        return self.mesh.comm.allreduce(glb_cost_proc, op=MPI.SUM )
+        return self.Vh[STATE].mesh.comm.allreduce(glb_cost_proc, op=MPI.SUM )
 
-
-    # def grad(self, m, out):
-    #  1. Cast the petsc4py vector m to a dlx.Function mfun
-    #  2. call symbolic differentation of self.functional_handler(mfun) wrt mfun
-    #  3. call assemble, update ghosts, and store the result in out
 
     def grad(self, x : dlx.la.Vector, out: dlx.la.Vector) -> None:
-        # mfun = vector2Function(x,self.Vh)
 
         updateFromVector(self.xfun, x)
         mfun = self.xfun
-
-
         L = dlx.fem.form(ufl.derivative(self.functional_handler(mfun),mfun,ufl.TestFunction(self.Vh)))
-
         out.array[:] = 0.
         tmp_out = dlx.la.create_petsc_vector_wrap(out)
         tmp_out.ghostUpdate(addv=petsc4py.PETSc.InsertMode.ADD_VALUES, mode=petsc4py.PETSc.ScatterMode.REVERSE)
         tmp_out.destroy()
 
-    # def setLinearizationPoint(self, m):
-    #   1. Cast the petsc4py vector m to a dlx.Function mfun
-    #   2. call symbolic differentiation (twice) to get the second variation of self.functional_handler(mfun) wrt mfun
-    #   3. assemble the Hessian operator (it's a sparse matrix!) in the attribute self.R
-    #   4. set up a linearsolver self.Rsolver that uses CG as Krylov method, gamg as preconditioner and self.R as operator
-
     def setLinearizationPoint(self, m, rel_tol=1e-12, max_iter=1000):
-        # mfun = vector2Function(m,self.Vh)
 
         updateFromVector(self.xfun, m)
         mfun = self.xfun
-
         L = ufl.derivative(ufl.derivative(self.functional_handler(mfun),mfun), mfun)
         self.R = dlx.fem.petsc.assemble_matrix(dlx.fem.form(L))
         self.R.assemble()
-
-        #for Rsolver:
-        #have to call _BilaplacianRsolver(self.Asolver,self.M)
-        #so have to construct self.Asolver, self.M first
-
         self.Rsolver = petsc4py.PETSc.KSP().create()
         self.Rsolver.getPC().setType(petsc4py.PETSc.PC.Type.GAMG)
         self.Rsolver.setType(petsc4py.PETSc.KSP.Type.CG)

hippylibX/modeling/misfit.py

@@ -6,9 +6,7 @@
 import numpy as np
 
 class NonGaussianContinuousMisfit(object):
-    def __init__(self, mesh : dlx.mesh.Mesh, Vh : list, form):
-        #mesh needed for comm.allreduce in cost function
-        self.mesh = mesh
+    def __init__(self,Vh : list, form):
         self.Vh = Vh
         self.form = form
 
@@ -28,7 +26,7 @@ def cost(self,x : list) -> float:
 
         loc_cost = self.form(u_fun,m_fun)
         glb_cost_proc = dlx.fem.assemble_scalar(dlx.fem.form(loc_cost))
-        return self.mesh.comm.allreduce(glb_cost_proc, op=MPI.SUM )
+        return self.Vh[hpx.STATE].mesh.comm.allreduce(glb_cost_proc, op=MPI.SUM )
 
     def grad(self, i : int, x : list, out: dlx.la.Vector) -> dlx.la.Vector: