From f4c63d47a3ae58aae1a5a3520b16e86e8cc71de6 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 24 Oct 2013 14:12:37 -0700 Subject: [PATCH] took for loop out of inversion framework. This has to be dealt with in the specific code, and is more flexible. --- SimPEG/forward/DCProblem.py | 91 +++++++++++++++++++++++++------------ SimPEG/inverse/Inversion.py | 11 ++--- SimPEG/inverse/Optimize.py | 10 +++- 3 files changed, 75 insertions(+), 37 deletions(-) diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index c2494daf..48f25f12 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -1,7 +1,8 @@ from SimPEG.mesh import TensorMesh from SimPEG.forward import Problem, SyntheticProblem from SimPEG.tests import checkDerivative -from SimPEG.utils import ModelBuilder, sdiag +from SimPEG.utils import ModelBuilder, sdiag, mkvc +from SimPEG import Solver import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as linalg @@ -48,7 +49,7 @@ class DCProblem(Problem): return phi - def J(self, m, v, u=None, solve=None): + def J(self, m, v, u=None): """ :param numpy.array m: model :param numpy.array v: vector to multiply @@ -70,6 +71,9 @@ class DCProblem(Problem): J(v) = - P ( A(m)^{-1} ( G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) v ) ) """ + if u is None: + u = self.field(m) + P = self.P D = self.mesh.faceDiv G = self.mesh.cellGrad @@ -78,15 +82,19 @@ class DCProblem(Problem): mT_dm = self.modelTransformDeriv(m) dCdu = A - dCdm = D * ( sdiag( G * u ) * ( Av_dm * ( mT_dm * v ) ) ) - if solve is None: - solve = linalg.factorized(dCdu) + dCdm = np.empty_like(u) + for i, ui in enumerate(u.T): # loop over each column + dCdm[:, i] = D * ( sdiag( G * ui ) * ( Av_dm * ( mT_dm * v ) ) ) - Jv = - P * solve(dCdm) + solve = Solver(dCdu) + # solve = linalg.factorized(dCdu) + + Jv = - P * solve.solve(dCdm) return Jv - def Jt(self, m, v, u=None, solve=None): + def Jt(self, m, v, u=None): + """Takes data, turns it into a model..ish""" P = self.P D = self.mesh.faceDiv G = self.mesh.cellGrad @@ -95,12 +103,15 @@ class DCProblem(Problem): mT_dm = self.modelTransformDeriv(m) dCdu = A.T + solve = Solver(dCdu) - if solve is None: - solve = linalg.factorized(dCdu.tocsc()) - w = solve(P.T*v) + w = solve.solve(P.T*v) - Jtv = - mT_dm.T * ( Av_dm.T * ( sdiag( G * u ) * ( D.T * w ) ) ) + Jtv = 0 + for i, ui in enumerate(u.T): # loop over each column + Jtv += sdiag( G * ui ) * ( D.T * w[:,i] ) + + Jtv = - mT_dm.T * ( Av_dm.T * Jtv ) return Jtv @@ -138,6 +149,7 @@ if __name__ == '__main__': from SimPEG.regularization import Regularization from SimPEG import inverse + import matplotlib.pyplot as plt # Create the mesh h1 = np.ones(100) @@ -145,16 +157,16 @@ if __name__ == '__main__': mesh = TensorMesh([h1,h2]) # Create some parameters for the model - sig1 = 1 - sig2 = 0.01 + sig1 = np.log(1) + sig2 = np.log(0.01) # Create a synthetic model from a block in a half-space p0 = [20, 20] p1 = [50, 50] condVals = [sig1, sig2] mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals) - mesh.plotImage(mSynth, showIt=False) - + plt.colorbar(mesh.plotImage(mSynth)) + # plt.show() # Set up the projection nelec = 50 @@ -185,30 +197,51 @@ if __name__ == '__main__': problem.RHS = q problem.dobs = dobs problem.std = dobs*0 + 0.05 - m0 = mesh.gridCC[:,0]*0+sig1 + m0 = mesh.gridCC[:,0]*0+sig2 - # print problem.misfit(m0) - # print problem.misfit(mSynth) - opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=1) + + # Adjoint Test + u = np.random.rand(mesh.nC, problem.RHS.shape[1]) + v = np.random.rand(mesh.nC) + w = np.random.rand(*dobs.shape) + Jv = mkvc(problem.J(mSynth, v, u=u)) + print mkvc(w).dot(Jv) + print v.dot(problem.Jt(mSynth, w, u=u)) + + # Check Derivative + dm = np.random.randn(*m0.shape) + for alp in np.logspace(-2,-6, 5): + a = problem.dpred(m0) + b = problem.dpred(m0 + alp*dm) + c = problem.J(m0, alp*dm) + print np.linalg.norm(a-b), np.linalg.norm(a-b+c) + + + # derChk = lambda m: [problem.dpred(m), problem.J(mSynth,m)] + # checkDerivative(derChk, mSynth) + + + opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=3) reg = Regularization(mesh) inv = inverse.Inversion(problem, reg, opt) - m = inv.run(m0) - - mesh.plotImage(m,showIt=True) - # Check Derivative derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)] checkDerivative(derChk, mSynth) - # Adjoint Test - # u = np.random.rand(mesh.nC) - # v = np.random.rand(mesh.nC) - # w = np.random.rand(dobs.shape[0]) - # print w.dot(problem.J(mSynth, v, u=u)) - # print v.dot(problem.Jt(mSynth, w, u=u)) + + print inv.dataObj(m0) + print inv.dataObj(mSynth) + + m = inv.run(m0) + + plt.colorbar(mesh.plotImage(m)) + print m + plt.show() + + diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index f217be83..fd166318 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -11,6 +11,7 @@ class Inversion(object): self.prob = prob self.reg = reg self.opt = opt + self.opt.parent = self @property def Wd(self): @@ -42,7 +43,7 @@ class Inversion(object): return m def getBeta(self): - return 1e3 + return 1e2 def stoppingCriteria(self): self._STOP = np.zeros(2,dtype=bool) @@ -138,9 +139,7 @@ class Inversion(object): R = self.Wd*self.prob.misfit(m, u=u) - dmisfit = 0 - for i in range(self.prob.RHS.shape[1]): # Loop over each right hand side - dmisfit += self.prob.Jt(m, self.Wd[:,i]*R[:,i], u=u[:,i]) + dmisfit = self.prob.Jt(m, self.Wd * R, u=u) return dmisfit @@ -182,8 +181,6 @@ class Inversion(object): R = self.Wd*self.prob.misfit(m, u=u) - dmisfit = 0 - for i in range(self.prob.RHS.shape[1]): # Loop over each right hand side - dmisfit += self.prob.Jt(m, self.Wd[:,i] * self.Wd[:,i] * self.prob.J(m, v, u=u[:,i]), u=u[:,i]) + dmisfit = self.prob.Jt(m, self.Wd * self.Wd * self.prob.J(m, v, u=u), u=u) return dmisfit diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 1bfda72c..2545a381 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -63,6 +63,14 @@ class Minimize(object): return self.xc + @property + def parent(self): + """This is the parent of the optimization routine.""" + return getattr(self, '_parent', None) + @parent.setter + def parent(self, value): + self._parent = value + def startup(self, x0): self._iter = 0 self._iterLS = 0 @@ -150,7 +158,7 @@ class GaussNewton(Minimize): class InexactGaussNewton(Minimize): name = 'InexactGaussNewton' def findSearchDirection(self): - p, info = sp.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=10) + p, info = sp.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=5) return p