From 2a8f43aa1015cb27e242684ac51fd3efb3def9c3 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 22 Oct 2013 19:42:38 -0700 Subject: [PATCH] Inversion Framework - a start.. --- SimPEG/forward/Problem.py | 139 +++--------------- SimPEG/inverse/Inversion.py | 179 ++++++++++++++++++++++++ SimPEG/inverse/Optimize.py | 54 ++++--- SimPEG/regularization/Regularization.py | 113 +++++++++++++++ SimPEG/tests/TestUtils.py | 10 +- 5 files changed, 341 insertions(+), 154 deletions(-) create mode 100644 SimPEG/inverse/Inversion.py create mode 100644 SimPEG/regularization/Regularization.py diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 5b716f1f..04f9771a 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -49,16 +49,6 @@ class Problem(object): def RHS(self, value): self._RHS = value - @property - def W(self): - """ - Standard deviation weighting matrix. - """ - return self._W - @W.setter - def W(self, value): - self._W = value - @property def P(self): """ @@ -83,16 +73,24 @@ class Problem(object): def dobs(self, value): self._dobs = value - def evalFunction(self, m, doDerivative=True): + def misfit(self, m, u=None): """ - :param numpy.array m: model - :param bool doDerivative: do you want to compute the derivative? - :rtype: numpy.array - :return: Jv - """ - f = self.misfit(m) + :param numpy.array m: geophysical model + :param numpy.array u: fields + :rtype: float + :return: data misfit - return f, g, H + The data misfit: + + .. math:: + + \mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs} + + Where P is a projection matrix that brings the field on the full domain to the data measurement locations; + u is the field of interest; d_obs is the observed data. + """ + + return self.dpred(m, u=u) - self.dobs def J(self, m, v, u=None): """ @@ -201,112 +199,7 @@ class Problem(object): """ return sdiag(np.exp(mkvc(m))) - def misfit(self, m, u=None): - """ - :param numpy.array m: geophysical model - :param numpy.array u: fields - :rtype: float - :return: data misfit - The data misfit using an l_2 norm is: - - .. math:: - - \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 - - Where P is a projection matrix that brings the field on the full domain to the data measurement locations; - u is the field of interest; d_obs is the observed data; and W is the weighting matrix. - """ - - R = self.W*(self.dpred(m, u=u) - self.dobs) - R = mkvc(R) - return 0.5*R.dot(R) - - def misfitDeriv(self, m, u=None): - """ - :param numpy.array m: geophysical model - :param numpy.array u: fields - :rtype: numpy.array - :return: data misfit derivative - - The data misfit using an l_2 norm is: - - .. math:: - - \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 - - If the field, u, is provided, the calculation of the data is fast: - - .. math:: - - \mathbf{d}_\\text{pred} = \mathbf{Pu(m)} - - \mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) - - Where P is a projection matrix that brings the field on the full domain to the data measurement locations; - u is the field of interest; d_obs is the observed data; and W is the weighting matrix. - - The derivative of this, with respect to the model, is: - - .. math:: - - \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R} - - """ - if u is None: - u = self.field(m) - - R = self.W*(self.dpred(m, u=u) - self.dobs) - - dmisfit = 0 - for i in range(self.RHS.shape[1]): # Loop over each right hand side - dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i]) - - return dmisfit - - def misfitDerivDeriv(self, m, u=None): - """ - :param numpy.array m: geophysical model - :param numpy.array u: fields - :rtype: numpy.array - :return: data misfit derivative - - The data misfit using an l_2 norm is: - - .. math:: - - \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 - - If the field, u, is provided, the calculation of the data is fast: - - .. math:: - - \mathbf{d}_\\text{pred} = \mathbf{Pu(m)} - - \mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) - - Where P is a projection matrix that brings the field on the full domain to the data measurement locations; - u is the field of interest; d_obs is the observed data; and W is the weighting matrix. - - The derivative of this, with respect to the model, is: - - .. math:: - - \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R} - - \\frac{\partial^2 \mu_\\text{data}}{\partial^2 \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ W J} - - """ - if u is None: - u = self.field(m) - - R = self.W*(self.dpred(m, u=u) - self.dobs) - - dmisfit = 0 - for i in range(self.RHS.shape[1]): # Loop over each right hand side - dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i]) - - return dmisfit class SyntheticProblem(object): diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py new file mode 100644 index 00000000..99173af0 --- /dev/null +++ b/SimPEG/inverse/Inversion.py @@ -0,0 +1,179 @@ +import numpy as np + +class Inversion(object): + """docstring for Inversion""" + + maxIter = 10 + + def __init__(self, prob, reg, opt): + self.prob = prob + self.reg = reg + self.opt = opt + + @property + def W(self): + """ + Standard deviation weighting matrix. + """ + return self._W + @W.setter + def W(self, value): + self._W = value + + def run(self, m0): + self._iter = 0 + while True: + self._beta = self.getBeta() + self.opt.minimize(self.evalFunction,m) + if self.stoppingCriteria(): break + self._iter += 1 + + def getBeta(self): + return 1 + + def stoppingCriteria(self): + self._STOP = np.zeros(2,dtype=bool) + self._STOP[0] = self._iter >= maxIter + self._STOP[1] = self._phi_d_last <= self.phi_d_target + return np.any(self._STOP) + + + def evalFunction(self, m, return_g=True, return_H=True): + + u = self.prob.field(m) + phi_d = self.dataObj(m, u) + phi_m = self.modelObj(m) + + self._phi_d_last = phi_d + self._phi_m_last = phi_m + + f = phi_d + self._beta * phi_m + + out = (f,) + if return_g: + phi_dDeriv = self.dataObjDeriv(m, u) + phi_mDeriv = self.modelObjDeriv(m) + + g = phi_dDeriv + self._beta * phi_mDeriv + out += (g,) + + if return_H: + def H_fun(v): + phi_d2Deriv = self.dataObj2Deriv(m, u, v) + phi_m2Deriv = self.modelObj2Deriv(m)*v + + return phi_d2Deriv + self._beta * phi_m2Deriv + + out += (H_fun,) + return out + + + def modelObj(self, m, u=None): + self.reg.misfit(m) + + + def dataObj(self, m, u=None): + """ + :param numpy.array m: geophysical model + :param numpy.array u: fields + :rtype: float + :return: data misfit + + The data misfit using an l_2 norm is: + + .. math:: + + \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 + + Where P is a projection matrix that brings the field on the full domain to the data measurement locations; + u is the field of interest; d_obs is the observed data; and W is the weighting matrix. + """ + R = self.Wd*self.prob.misfit(u=u) + R = mkvc(R) + return 0.5*R.dot(R) + + def dataObjDeriv(self, m, u=None): + """ + :param numpy.array m: geophysical model + :param numpy.array u: fields + :rtype: numpy.array + :return: data misfit derivative + + The data misfit using an l_2 norm is: + + .. math:: + + \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 + + If the field, u, is provided, the calculation of the data is fast: + + .. math:: + + \mathbf{d}_\\text{pred} = \mathbf{Pu(m)} + + \mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) + + Where P is a projection matrix that brings the field on the full domain to the data measurement locations; + u is the field of interest; d_obs is the observed data; and W is the weighting matrix. + + The derivative of this, with respect to the model, is: + + .. math:: + + \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R} + + """ + if u is None: + u = self.field(m) + + R = self.W*(self.dpred(m, u=u) - self.dobs) + + dmisfit = 0 + for i in range(self.RHS.shape[1]): # Loop over each right hand side + dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i]) + + return dmisfit + + def dataObj2Deriv(self, m, u=None): + """ + :param numpy.array m: geophysical model + :param numpy.array u: fields + :rtype: numpy.array + :return: data misfit derivative + + The data misfit using an l_2 norm is: + + .. math:: + + \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 + + If the field, u, is provided, the calculation of the data is fast: + + .. math:: + + \mathbf{d}_\\text{pred} = \mathbf{Pu(m)} + + \mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) + + Where P is a projection matrix that brings the field on the full domain to the data measurement locations; + u is the field of interest; d_obs is the observed data; and W is the weighting matrix. + + The derivative of this, with respect to the model, is: + + .. math:: + + \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R} + + \\frac{\partial^2 \mu_\\text{data}}{\partial^2 \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ W J} + + """ + if u is None: + u = self.field(m) + + R = self.W*(self.dpred(m, u=u) - self.dobs) + + dmisfit = 0 + for i in range(self.RHS.shape[1]): # Loop over each right hand side + dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i]) + + return dmisfit diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 965fb308..c3ad8dd2 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -23,8 +23,9 @@ class Minimize(object): tolG = 1e-4 eps = 1e-16 - def __init__(self, problem, **kwargs): - self.problem = problem + printIter = [] # push to here if you want to print these on iter + + def __init__(self, **kwargs): self.setKwargs(**kwargs) def setKwargs(self, **kwargs): @@ -35,13 +36,20 @@ class Minimize(object): else: raise Exception('%s attr is not recognized' % attr) - def minimize(self, x0): + def minimize(self, evalFunction, x0): + """ + evalFunction is a function handle:: + + evalFunction(x, return_g=True, return_H=True ) + + """ + self.evalFunction = evalFunction self.startup(x0) self.printInit() while True: - self.f, self.g, self.H = self.evalFunction(self.xc) + self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) self.printIter() if self.stoppingCriteria(): break p = self.findSearchDirection() @@ -67,31 +75,17 @@ class Minimize(object): """ printIter is called at the beginning of the optimization routine. - If the problem object has a printInit function it will be called here:: - - self.problem.printInit(self) - """ - if hasattr(self.problem, 'printInit'): - self.problem.printInit(self) - else: - print "%s %s %s" % ('='*22, self.name, '='*22) - print "iter\tJc\t\tnorm(dJ)\tLS" - print "%s" % '-'*57 + print "%s %s %s" % ('='*22, self.name, '='*22) + print "iter\tJc\t\tnorm(dJ)\tLS" + print "%s" % '-'*57 def printIter(self): """ printIter is called directly after function evaluations. - If the problem object has a printIter function it will be called here:: - - self.problem.printIter(self) - """ - if hasattr(self.problem, 'printIter'): - self.problem.printIter(self) - else: - print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS) + print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS) def printDone(self): print "%s STOP! %s" % ('-'*25,'-'*25) @@ -102,10 +96,6 @@ class Minimize(object): print "%d : iter = %3d\t <= maxIter\t = %3d" % (self._STOP[4], self._iter, self.maxIter) print "%s DONE! %s\n" % ('='*25,'='*25) - def evalFunction(self, x, doDerivative=True): - f, g, H = self.problem(x) - return f, g, H - def findSearchDirection(self): return -self.g @@ -128,7 +118,7 @@ class Minimize(object): iterLS = 0 while iterLS < self.maxIterLS: xt = self.xc + t*p - ft, temp, temp = self.evalFunction(xt, doDerivative=False) + ft = self.evalFunction(xt, return_g=False, return_H=False) if ft < self.f + t*self.LSreduction*descent: break iterLS += 1 @@ -153,6 +143,12 @@ class GaussNewton(Minimize): return np.linalg.solve(self.H,-self.g) +class InexactGaussNewton(Minimize): + name = 'InexactGaussNewton' + def findSearchDirection(self): + return sparse.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=10) + + class SteepestDescent(Minimize): name = 'SteepestDescent' def findSearchDirection(self): @@ -162,9 +158,9 @@ if __name__ == '__main__': from SimPEG.tests import Rosenbrock, checkDerivative x0 = np.array([2.6, 3.7]) checkDerivative(Rosenbrock, x0, plotIt=False) - xOpt = GaussNewton(Rosenbrock, maxIter=20).minimize(x0) + xOpt = GaussNewton(maxIter=20).minimize(Rosenbrock,x0) print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) - xOpt = SteepestDescent(Rosenbrock, maxIter=20, maxIterLS=15).minimize(x0) + xOpt = SteepestDescent(maxIter=20, maxIterLS=15).minimize(Rosenbrock, x0) print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) def simplePass(x): diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py new file mode 100644 index 00000000..f0239875 --- /dev/null +++ b/SimPEG/regularization/Regularization.py @@ -0,0 +1,113 @@ +from SimPEG.utils import sdiag + +class Regularization(object): + """docstring for Regularization""" + + @property + def mref(self): + return self._mref + @mref.setter + def mref(self, value): + self._mref = value + + @property + def Wx(self): + if self._Wx is None: + self._Wx = mesh.cellGradx + return self._Wx + + @property + def Wy(self): + if self._Wy is None: + self._Wy = mesh.cellGrady + return self._Wy + + @property + def Wz(self): + if self._Wz is None: + self._Wz = mesh.cellGradz + return self._Wz + + @property + def Ws(self): + if self._Ws is None: + self._Ws = sdiag(self.mesh.vol) + return self._Ws + + + def __init__(self, mesh): + self.mesh = mesh + self._Wx = None + self._Wy = None + self._Wz = None + self.alpha_s = 1e-6 + self.alpha_x = 1 + self.alpha_y = 1 + self.alpha_z = 1 + + def pnorm(self, r): + return 0.5*r.dot(r) + + def modelObj(self, m): + mresid = m - self.mref + + mobj = self.alpha_s * self.pnorm( self.Ws * mresid ) + + mobj += self.alpha_x * self.pnorm( self.Wx * mresid ) + + if self.mesh.dim > 1: + mobj += self.alpha_y * self.pnorm( self.Wy * mresid ) + if self.mesh.dim > 2: + mobj += self.alpha_z * self.pnorm( self.Wz * mresid ) + + return mobj + + def modelObjDeriv(self, m): + """ + + In 1D: + + .. math:: + + m_{\\text{obj}} = {1 \over 2}\\alpha_s \left\| W_s (m- m_{\\text{ref}})\\right\|^2_2 + + {1 \over 2}\\alpha_x \left\| W_x (m- m_{\\text{ref}})\\right\|^2_2 + + \\frac{ \partial m_{\\text{obj}} }{\partial m} = + \\alpha_s W_s^{\\top} W_s (m - m_{\\text{ref}}) + + \\alpha_x W_x^{\\top} W_x (m - m_{\\text{ref}}) + + + \\frac{ \partial^2 m_{\\text{obj}} }{\partial m^2} = + \\alpha_s W_s^{\\top} W_s + + \\alpha_x W_x^{\\top} W_x + + """ + + mresid = m - self.mref + + mobjDeriv = self.alpha_s * self.Ws.T * ( self.Ws * mresid) + + mobjDeriv += self.alpha_x * self.Wx.T * ( self.Wx * mresid) + + if self.mesh.dim > 1: + mobjDeriv += self.alpha_y * self.Wy.T * ( self.Wy * mresid) + if self.mesh.dim > 2: + mobjDeriv += self.alpha_z * self.Wz.T * ( self.Wz * mresid) + + return mobjDeriv + + + def modelObj2Deriv(self, m): + mresid = m - self.mref + + mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws + + mobj2Deriv += self.alpha_x * self.Wx.T * self.Wx + + if self.mesh.dim > 1: + mobj2Deriv += self.alpha_y * self.Wy.T * self.Wy + if self.mesh.dim > 2: + mobj2Deriv += self.alpha_z * self.Wz.T * self.Wz + + return mobj2Deriv + diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index 9b2158c4..83738a75 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -163,13 +163,19 @@ class OrderTest(unittest.TestCase): print '' self.assertTrue(passTest) -def Rosenbrock(x): +def Rosenbrock(x, return_g=True, return_H=True): """Rosenbrock function for testing GaussNewton scheme""" f = 100*(x[1]-x[0]**2)**2+(1-x[0])**2 g = np.array([2*(200*x[0]**3-200*x[0]*x[1]+x[0]-1), 200*(x[1]-x[0]**2)]) H = np.array([[-400*x[1]+1200*x[0]**2+2, -400*x[0]], [-400*x[0], 200]]) - return f, g, H + + out = (f,) + if return_g: + out += (g,) + if return_H: + out += (H,) + return out def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): """