Skip to content

Commit

Permalink
checking CI for modifications to prior only
Browse files Browse the repository at this point in the history
  • Loading branch information
V-Rang authored and V-Rang committed Mar 26, 2024
1 parent 684e678 commit 164f928
Show file tree
Hide file tree
Showing 4 changed files with 50 additions and 146 deletions.
9 changes: 3 additions & 6 deletions hippylibX/algorithms/NewtonCG.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,9 +212,8 @@ def _solve_ls(self,x):
tolcg = min(cg_coarse_tolerance, math.sqrt(gradnorm/gradnorm_ini))

HessApply = ReducedHessian(self.model)
HessApply_wrap = HessApply.as_petsc_wrapper()
solver = CGSolverSteihaug(comm = self.model.prior.Vh.mesh.comm)
solver.set_operator(HessApply_wrap)
solver.set_operator(HessApply)
solver.set_preconditioner(self.model.Rsolver())
solver.parameters["rel_tolerance"] = tolcg
solver.parameters["max_iter"] = cg_max_iter
Expand Down Expand Up @@ -339,9 +338,8 @@ def _solve_tr(self,x):
tolcg = min(cg_coarse_tolerance, math.sqrt(gradnorm/gradnorm_ini))

HessApply = ReducedHessian(self.model)
HessApply_wrap = HessApply.as_petsc_wrapper()
solver = CGSolverSteihaug(comm = self.model.prior.Vh.mesh.comm)
solver.set_operator(HessApply_wrap)
solver.set_operator(HessApply)
solver.set_preconditioner(self.model.Rsolver())
if self.it > 1:
solver.set_TR(delta_TR, self.model.prior.R)
Expand Down Expand Up @@ -369,8 +367,7 @@ def _solve_tr(self,x):
#Calculate Predicted Reduction
H_mhat = self.model.generate_vector(PARAMETER)
H_mhat.array[:] = 0.
# HessApply.mult(mhat,H_mhat)
HessApply_wrap.mult(mhat,H_mhat)
HessApply.mult(mhat,H_mhat)
mg_mhat = mg.inner(mhat)
PRED_RED = -0.5*mhat.inner(H_mhat) - mg_mhat
# print( "PREDICTED REDUCTION", PRED_RED, "ACTUAL REDUCTION", ACTUAL_RED)
Expand Down
112 changes: 36 additions & 76 deletions hippylibX/algorithms/cgsolverSteihaug.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,14 +97,10 @@ def set_operator(self, A):
"""

self.A = A
# self.r = self.A.init_vector(0)
# self.z = self.A.init_vector(0)
# self.d = self.A.init_vector(0)
# self.Ad = self.A.init_vector(0)
self.r = self.A.createVecLeft()
self.z = self.A.createVecLeft()
self.d = self.A.createVecLeft()
self.Ad = self.A.createVecLeft()
self.r = self.A.init_vector(0)
self.z = self.A.init_vector(0)
self.d = self.A.init_vector(0)
self.Ad = self.A.init_vector(0)


def set_preconditioner(self, B_solver):
Expand Down Expand Up @@ -169,45 +165,27 @@ def solve(self,x,b):


if self.parameters["zero_initial_guess"]:
# self.r.array[:] = b.array
temp_petsc_vec_b = dlx.la.create_petsc_vector_wrap(b)
self.r.scale(0.)
self.r.axpy(1.,temp_petsc_vec_b)
temp_petsc_vec_b.destroy()
self.r.array[:] = b.array

x.array[:] = 0.

else:
assert self.TR_radius_2==None
temp_petsc_vec_x = dlx.la.create_petsc_vector_wrap(x)
# self.A.mult(x,self.r)
self.A.mult(temp_petsc_vec_x,self.r)
temp_petsc_vec_x.destroy()

# self.r.array[:] *= -1.
# self.r.array[:] += b.array
self.r.scale(-1.)
temp_petsc_vec_b = dlx.la.create_petsc_vector_wrap(b)
self.r.axpy(1.,temp_petsc_vec_b)
temp_petsc_vec_b.destroy()
self.A.mult(x,self.r)
self.r.array[:] *= -1.
self.r.array[:] += b.array

# self.z.array[:] = 0.
self.z.scale(0.)

# temp_self_z = dlx.la.create_petsc_vector_wrap(self.z)
# temp_self_r = dlx.la.create_petsc_vector_wrap(self.r)
# self.B_solver.solve(temp_self_r, temp_self_z) #expects (b,x)
# temp_self_z.destroy()
# temp_self_r.destroy()

self.B_solver.solve(self.r, self.z) #expects (b,x)
self.z.array[:] = 0.

temp_self_z = dlx.la.create_petsc_vector_wrap(self.z)
temp_self_r = dlx.la.create_petsc_vector_wrap(self.r)
self.B_solver.solve(temp_self_r, temp_self_z) #expects (b,x)
temp_self_z.destroy()
temp_self_r.destroy()

# self.d.array[:] = self.z.array
self.d.scale(0.)
self.d.axpy(1.,self.z)
self.d.array[:] = self.z.array

# nom0 = inner(self.d,self.r)
nom0 = self.d.dot(self.r)
nom0 = inner(self.d,self.r)

nom = nom0

Expand All @@ -228,30 +206,21 @@ def solve(self,x,b):

return


self.A.mult(self.d, self.Ad)
# den = inner(self.Ad, self.d)
den = self.Ad.dot(self.d)
den = inner(self.Ad, self.d)
if den <= 0.0:
self.converged = True
self.reasonid = 2
# x.array[:] += self.d.array
# self.r.array[:] -= self.Ad.array

temp_petsc_vec_x = dlx.la.create_petsc_vector_wrap(x)
temp_petsc_vec_x.axpy(1., self.d)
self.r.axpy(-1., self.Ad)
x.array[:] += self.d.array
self.r.array[:] -= self.Ad.array

# temp_self_z = dlx.la.create_petsc_vector_wrap(self.z)
# temp_self_r = dlx.la.create_petsc_vector_wrap(self.r)
# self.B_solver.solve(temp_self_r, temp_self_z) #expects (b,x)
# temp_self_z.destroy()
# temp_self_r.destroy()
temp_self_z = dlx.la.create_petsc_vector_wrap(self.z)
temp_self_r = dlx.la.create_petsc_vector_wrap(self.r)
self.B_solver.solve(temp_self_r, temp_self_z) #expects (b,x)
temp_self_z.destroy()
temp_self_r.destroy()

self.B_solver.solve(self.r, self.z) #expects (b,x)

# nom = inner(self.r, self.z)
nom = self.r.inner(self.z)
nom = inner(self.r, self.z)
self.final_norm = math.sqrt(nom)
if(self.parameters["print_level"] >= 0):
print( self.reason[self.reasonid])
Expand All @@ -273,19 +242,15 @@ def solve(self,x,b):
print( "Converged in ", self.iter, " iterations with final norm ", self.final_norm)
break

# self.r.array[:] -= alpha*self.Ad.array
self.r.axpy(-alpha,self.Ad)

# temp_self_z = dlx.la.create_petsc_vector_wrap(self.z)
# temp_self_r = dlx.la.create_petsc_vector_wrap(self.r)
# self.B_solver.solve(temp_self_r, temp_self_z) #expects (b,x)
# temp_self_z.destroy()
# temp_self_r.destroy()
self.r.array[:] -= alpha*self.Ad.array

self.B_solver.solve(self.r, self.z) #expects (b,x)
temp_self_z = dlx.la.create_petsc_vector_wrap(self.z)
temp_self_r = dlx.la.create_petsc_vector_wrap(self.r)
self.B_solver.solve(temp_self_r, temp_self_z) #expects (b,x)
temp_self_z.destroy()
temp_self_r.destroy()

# betanom = inner(self.r, self.z)
betanom = self.r.dot(self.z)
betanom = inner(self.r, self.z)
if self.parameters["print_level"] == 1:
print( " Iteration : ", self.iter, " (B r, r) = ", betanom)

Expand All @@ -311,16 +276,11 @@ def solve(self,x,b):
break

beta = betanom/nom
# self.d.array[:] = beta*self.d.array
# self.d.array[:] += self.z.array

self.d.scale(beta)
self.d.axpy(1., self.z)

self.d.array[:] = beta*self.d.array
self.d.array[:] += self.z.array
self.A.mult(self.d,self.Ad)

# den = inner(self.d, self.Ad)
den = self.d.dot(self.Ad)
den = inner(self.d, self.Ad)
if den <= 0.0:
self.converged = True
self.reasonid = 2
Expand Down
45 changes: 6 additions & 39 deletions hippylibX/modeling/modelVerify.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,23 +46,9 @@ def modelVerify(model, m0 : dlx.la.Vector, is_quadratic = False, misfit_only=Fal

model.setPointForHessianEvaluations(x)

H_object = ReducedHessian(model, misfit_only=misfit_only)
H = H_object.as_petsc_wrapper()


# H_obj = ReducedHessian(model, misfit_only=misfit_only)
# H = petsc4py.PETSc.Mat().createPython(model.prior.M.getSizes(),comm=model.prior.Vh.mesh.comm)
# H.setPythonContext(H_obj)
# H.setUp()

H = ReducedHessian(model, misfit_only=misfit_only)
Hh = model.generate_vector(PARAMETER)
temp_petsc_vec_h = dlx.la.create_petsc_vector_wrap(h)
temp_petsc_vec_Hh = dlx.la.create_petsc_vector_wrap(Hh)

H.mult(temp_petsc_vec_h, temp_petsc_vec_Hh)
# H.mult(h,Hh)
temp_petsc_vec_h.destroy()
temp_petsc_vec_Hh.destroy()
H.mult(h, Hh)

if eps is None:
n_eps = 32
Expand Down Expand Up @@ -114,32 +100,13 @@ def modelVerify(model, m0 : dlx.la.Vector, is_quadratic = False, misfit_only=Fal
xx = model.generate_vector(PARAMETER)
parRandom.normal(1., xx)


yy = model.generate_vector(PARAMETER)
parRandom.normal(1., yy)

# ytHx = H.inner(yy,xx)
# xtHy = H.inner(xx,yy)
ytHx = H.inner(yy,xx)

#######################################
temp_petsc_vec_xx = dlx.la.create_petsc_vector_wrap(xx)
temp_petsc_vec_yy = dlx.la.create_petsc_vector_wrap(yy)

# #in two separate operations:
Ay = model.generate_vector(PARAMETER)
temp_petsc_vec_Ay = dlx.la.create_petsc_vector_wrap(Ay)
temp_petsc_vec_Ay.scale(0.)
H.mult(temp_petsc_vec_xx,temp_petsc_vec_Ay)
ytHx = temp_petsc_vec_yy.dot(temp_petsc_vec_Ay)

temp_petsc_vec_Ay.scale(0.)
H.mult(temp_petsc_vec_yy, temp_petsc_vec_Ay)
xtHy = temp_petsc_vec_xx.dot(temp_petsc_vec_Ay)

temp_petsc_vec_Ay.destroy()
temp_petsc_vec_xx.destroy()
temp_petsc_vec_yy.destroy()
# #######################################

xtHy = H.inner(xx,yy)

if np.abs(ytHx + xtHy) > 0.:
rel_symm_error = 2*abs(ytHx - xtHy)/(ytHx + xtHy)
Expand Down Expand Up @@ -177,4 +144,4 @@ def modelVerifyPlotErrors(is_quadratic : bool, eps : np.ndarray, err_grad : np.n
plt.title("FD Gradient Check")
plt.subplot(122)
plt.loglog(eps, err_H, "-ob", eps, eps*(err_H[0]/eps[0]), "-.k")
plt.title("FD Hessian Check")
plt.title("FD Hessian Check")
30 changes: 5 additions & 25 deletions hippylibX/modeling/reducedHessian.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,6 @@
from ..algorithms import linalg
import numpy as np
from ..utils import vector2Function
import petsc4py

#decorator for functions in classes that are not used -> may not be needed in the final
#version of X
Expand Down Expand Up @@ -69,32 +68,18 @@ def init_vector(self, dim: int) -> dlx.la.Vector:
return self.model.init_parameter(dim)


def mult(self,mat, x : petsc4py.PETSc.Vec ,y : petsc4py.PETSc.Vec ) -> None:
def mult(self,x : dlx.la.Vector ,y : dlx.la.Vector ) -> None:

"""
Apply the reduced Hessian (or the Gauss-Newton approximation) to the vector :code:`x`. Return the result in :code:`y`.
"""

x_dlx = self.model.generate_vector(PARAMETER)
y_dlx = self.model.generate_vector(PARAMETER)

x_tmp_dlx = dlx.la.create_petsc_vector_wrap(x_dlx)
x_tmp_dlx.axpy(1.,x)

if self.gauss_newton_approx:
self.GNHessian(x_dlx,y_dlx)
self.GNHessian(x,y)
else:
self.TrueHessian(x_dlx,y_dlx)

tmp = dlx.la.create_petsc_vector_wrap(y_dlx)

y.axpby(1.,0.,tmp) #y = 1. tmp + 0.*y
tmp.destroy()
x_tmp_dlx.destroy()
self.TrueHessian(x,y)

self.ncalls += 1



def inner(self, x : dlx.la.Vector, y : dlx.la.Vector):
"""
Perform the inner product between :code:`x` and :code:`y` in the norm induced by the reduced
Expand Down Expand Up @@ -150,9 +135,4 @@ def TrueHessian(self, x : dlx.la.Vector , y : dlx.la.Vector) -> None:
self.model.applyR(x,self.yhelp)
y.array[:] += self.yhelp.array


def as_petsc_wrapper(self) -> petsc4py.PETSc.Mat:
H = petsc4py.PETSc.Mat().createPython(self.model.prior.M.getSizes(),comm=self.model.prior.Vh.mesh.comm)
H.setPythonContext(self)
H.setUp()
return H

0 comments on commit 164f928

Please sign in to comment.