From 41b7dbe82b5415aa26bdda574fa20b2a046953bc Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 22 Oct 2013 14:40:52 -0700 Subject: [PATCH 01/34] bug fix --- docs/examples/mesh/plot_LogicallyOrthogonalMesh.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/examples/mesh/plot_LogicallyOrthogonalMesh.py b/docs/examples/mesh/plot_LogicallyOrthogonalMesh.py index 55350946..bb49ae9a 100644 --- a/docs/examples/mesh/plot_LogicallyOrthogonalMesh.py +++ b/docs/examples/mesh/plot_LogicallyOrthogonalMesh.py @@ -1,4 +1,5 @@ -from SimPEG import LogicallyOrthogonalMesh, utils +from SimPEG.mesh import LogicallyOrthogonalMesh +from SimPEG import utils import matplotlib.pyplot as plt X, Y = utils.exampleLomGird([3,3],'rotate') M = LogicallyOrthogonalMesh([X, Y]) From 13efdc12f90f6ab8723f3511ddd35a9118c65649 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 22 Oct 2013 14:42:30 -0700 Subject: [PATCH 02/34] updates to view, and documentation of base mesh. --- SimPEG/mesh/BaseMesh.py | 25 +++++++++++++++++++++++++ SimPEG/mesh/TensorView.py | 3 +++ 2 files changed, 28 insertions(+) diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py index 29c43dd1..6a9a8032 100644 --- a/SimPEG/mesh/BaseMesh.py +++ b/SimPEG/mesh/BaseMesh.py @@ -2,6 +2,7 @@ import numpy as np from SimPEG.utils import mkvc + class BaseMesh(object): """ BaseMesh does all the counting you don't want to do. @@ -216,6 +217,12 @@ class BaseMesh(object): :rtype: int :return: nC + + .. plot:: + + from SimPEG.mesh import TensorMesh + import numpy as np + TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(centers=True,showIt=True) """ fget = lambda self: np.prod(self.n) return locals() @@ -270,6 +277,12 @@ class BaseMesh(object): :rtype: int :return: nN + + .. plot:: + + from SimPEG.mesh import TensorMesh + import numpy as np + TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(nodes=True,showIt=True) """ fget = lambda self: np.prod(self.n + 1) return locals() @@ -324,6 +337,12 @@ class BaseMesh(object): :rtype: numpy.array (dim, ) :return: [prod(nEx), prod(nEy), prod(nEz)] + + .. plot:: + + from SimPEG.mesh import TensorMesh + import numpy as np + TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(edges=True,showIt=True) """ fget = lambda self: np.array([np.prod(x) for x in [self.nEx, self.nEy, self.nEz] if not x is None]) return locals() @@ -378,6 +397,12 @@ class BaseMesh(object): :rtype: numpy.array (dim, ) :return: [prod(nFx), prod(nFy), prod(nFz)] + + .. plot:: + + from SimPEG.mesh import TensorMesh + import numpy as np + TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(faces=True,showIt=True) """ fget = lambda self: np.array([np.prod(x) for x in [self.nFx, self.nFy, self.nFz] if not x is None]) return locals() diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py index 687a520d..0b9ff7b3 100644 --- a/SimPEG/mesh/TensorView.py +++ b/SimPEG/mesh/TensorView.py @@ -267,6 +267,9 @@ class TensorView(object): if faces: ax.plot(xs1[:, 0], xs1[:, 1], 'g>') ax.plot(xs2[:, 0], xs2[:, 1], 'g^') + if edges: + ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'c>') + ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'c^') # Plot the grid lines if lines: From 0cb9fa210acec6979ea9deb6134c4907486e5420 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 22 Oct 2013 14:43:24 -0700 Subject: [PATCH 03/34] print statements now go to problem if they are implemented there. --- SimPEG/inverse/Optimize.py | 37 ++++++++++++++++++++++++++++++++----- 1 file changed, 32 insertions(+), 5 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index eee18b16..965fb308 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -5,7 +5,12 @@ norm = np.linalg.norm class Minimize(object): - """docstring for Minimize""" + """ + + Minimize is a general class for derivative based optimization. + + + """ name = "GeneralOptimizationAlgorithm" @@ -59,12 +64,34 @@ class Minimize(object): self.xOld = x0 def printInit(self): - print "%s %s %s" % ('='*22, self.name, '='*22) - print "iter\tJc\t\tnorm(dJ)\tLS" - print "%s" % '-'*57 + """ + 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 def printIter(self): - print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS) + """ + 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) def printDone(self): print "%s STOP! %s" % ('-'*25,'-'*25) From 2a8f43aa1015cb27e242684ac51fd3efb3def9c3 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 22 Oct 2013 19:42:38 -0700 Subject: [PATCH 04/34] 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): """ From 2bb20280bfe6110c18f69dba76cfdeb031eec7ba Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 22 Oct 2013 20:11:30 -0700 Subject: [PATCH 05/34] to get working.. --- SimPEG/forward/DCProblem/DCProblem.py | 19 ++++++++++++++++--- SimPEG/forward/Problem.py | 9 +++++++++ SimPEG/inverse/Inversion.py | 11 +++++------ SimPEG/inverse/Optimize.py | 5 ++++- SimPEG/inverse/__init__.py | 1 + SimPEG/regularization/__init__.py | 1 + 6 files changed, 36 insertions(+), 10 deletions(-) create mode 100644 SimPEG/regularization/__init__.py diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py index 1df8897b..45b8a2c1 100644 --- a/SimPEG/forward/DCProblem/DCProblem.py +++ b/SimPEG/forward/DCProblem/DCProblem.py @@ -105,6 +105,10 @@ class DCProblem(Problem): if __name__ == '__main__': + + from SimPEG.regularization import Regularization + from SimPEG import inverse + # Create the mesh h1 = np.ones(100) h2 = np.ones(100) @@ -143,7 +147,7 @@ if __name__ == '__main__': dobs, Wd = synthetic.createData(mSynth, std=0.05) u = synthetic.field(mSynth) - mesh.plotImage(u[:,10], showIt=True) + mesh.plotImage(u[:,10], showIt=False) # Now set up the problem to do some minimization problem = DCProblem(mesh) @@ -153,8 +157,15 @@ if __name__ == '__main__': problem.dobs = dobs m0 = mesh.gridCC[:,0]*0+sig1 - print problem.misfit(m0) - print problem.misfit(mSynth) + # print problem.misfit(m0) + # print problem.misfit(mSynth) + + opt = inverse.InexactGaussNewton() + reg = Regularization(mesh) + + inv = inverse.Inversion(problem, reg, opt) + + inv.run(m0) # Check Derivative derChk = lambda m: [problem.misfit(m), problem.misfitDeriv(m)] @@ -166,3 +177,5 @@ if __name__ == '__main__': 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)) + + diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 04f9771a..0333b8b4 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -62,6 +62,15 @@ class Problem(object): def P(self, value): self._P = value + @property + def std(self): + """ + Estimated Standard Deviations. + """ + return self._std + @std.setter + def std(self, value): + self._std = value @property def dobs(self): diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 99173af0..8d92e15e 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,4 +1,5 @@ import numpy as np +from SimPEG.utils import sdiag class Inversion(object): """docstring for Inversion""" @@ -11,20 +12,18 @@ class Inversion(object): self.opt = opt @property - def W(self): + def Wd(self): """ Standard deviation weighting matrix. """ - return self._W - @W.setter - def W(self, value): - self._W = value + return sdiag(1/self.prob.std) def run(self, m0): + m = m0 self._iter = 0 while True: self._beta = self.getBeta() - self.opt.minimize(self.evalFunction,m) + m = self.opt.minimize(self.evalFunction,m) if self.stoppingCriteria(): break self._iter += 1 diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index c3ad8dd2..43776d33 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -111,13 +111,16 @@ class Minimize(object): self._STOP[4] = self._iter >= self.maxIter return all(self._STOP[0:3]) | any(self._STOP[3:]) + def projection(self, p): + return p + def linesearch(self, p): # Armijo linesearch descent = np.inner(self.g, p) t = 1 iterLS = 0 while iterLS < self.maxIterLS: - xt = self.xc + t*p + xt = self.projection(self.xc + t*p) ft = self.evalFunction(xt, return_g=False, return_H=False) if ft < self.f + t*self.LSreduction*descent: break diff --git a/SimPEG/inverse/__init__.py b/SimPEG/inverse/__init__.py index b2a5e506..c4ca335b 100644 --- a/SimPEG/inverse/__init__.py +++ b/SimPEG/inverse/__init__.py @@ -1 +1,2 @@ from Optimize import * +from Inversion import * diff --git a/SimPEG/regularization/__init__.py b/SimPEG/regularization/__init__.py new file mode 100644 index 00000000..0230f8c3 --- /dev/null +++ b/SimPEG/regularization/__init__.py @@ -0,0 +1 @@ +from Regularization import Regularization From 3de8ab47d925ec2fee558d9984b10fd90e3ce2ca Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 22 Oct 2013 20:25:49 -0700 Subject: [PATCH 06/34] moved DCProblem --- SimPEG/__init__.py | 2 ++ SimPEG/forward/{DCProblem => }/DCProblem.py | 33 +++++++++++++++++++-- SimPEG/forward/DCProblem/DCutils.py | 29 ------------------ SimPEG/forward/DCProblem/__init__.py | 2 -- SimPEG/forward/__init__.py | 2 +- 5 files changed, 34 insertions(+), 34 deletions(-) rename SimPEG/forward/{DCProblem => }/DCProblem.py (82%) delete mode 100644 SimPEG/forward/DCProblem/DCutils.py delete mode 100644 SimPEG/forward/DCProblem/__init__.py diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index e8946f88..a17eaa74 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1,4 +1,6 @@ import mesh import utils import inverse +import forward +import regularization from Solver import Solver diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem.py similarity index 82% rename from SimPEG/forward/DCProblem/DCProblem.py rename to SimPEG/forward/DCProblem.py index 45b8a2c1..5aad45bd 100644 --- a/SimPEG/forward/DCProblem/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -3,8 +3,8 @@ from SimPEG.forward import Problem, SyntheticProblem from SimPEG.tests import checkDerivative from SimPEG.utils import ModelBuilder, sdiag import numpy as np +import scipy as sp import scipy.sparse.linalg as linalg -import DCutils class DCProblem(Problem): """ @@ -134,7 +134,7 @@ if __name__ == '__main__': elecend = 0.5+spacelec*(nelec-1) elecLocR = np.linspace(elecini, elecend, nelec) rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5 - q, Q, rxmidloc = DCutils.genTxRxmat(nelec, spacelec, surfloc, elecini, mesh) + q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, mesh) P = Q.T # Create some data @@ -179,3 +179,32 @@ if __name__ == '__main__': print v.dot(problem.Jt(mSynth, w, u=u)) + + + +def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh): + """ Generate projection matrix (Q) and """ + elecend = 0.5+spacelec*(nelec-1) + elecLocR = np.linspace(elecini, elecend, nelec) + elecLocT = elecLocR+1 + nrx = nelec-1 + ntx = nelec-1 + q = np.zeros((mesh.nC, ntx)) + Q = np.zeros((mesh.nC, nrx)) + + for i in range(nrx): + + rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i])) + rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1])) + + txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i])) + txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1])) + + q[txind1,i] = 1 + q[txind2,i] = -1 + Q[rxind1,i] = 1 + Q[rxind2,i] = -1 + + Q = sp.csr_matrix(Q) + rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5 + return q, Q, rxmidLoc diff --git a/SimPEG/forward/DCProblem/DCutils.py b/SimPEG/forward/DCProblem/DCutils.py deleted file mode 100644 index f3445096..00000000 --- a/SimPEG/forward/DCProblem/DCutils.py +++ /dev/null @@ -1,29 +0,0 @@ -import numpy as np -import scipy.sparse as sp - -def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh): - """ Generate projection matrix (Q) and """ - elecend = 0.5+spacelec*(nelec-1) - elecLocR = np.linspace(elecini, elecend, nelec) - elecLocT = elecLocR+1 - nrx = nelec-1 - ntx = nelec-1 - q = np.zeros((mesh.nC, ntx)) - Q = np.zeros((mesh.nC, nrx)) - - for i in range(nrx): - - rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i])) - rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1])) - - txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i])) - txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1])) - - q[txind1,i] = 1 - q[txind2,i] = -1 - Q[rxind1,i] = 1 - Q[rxind2,i] = -1 - - Q = sp.csr_matrix(Q) - rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5 - return q, Q, rxmidLoc diff --git a/SimPEG/forward/DCProblem/__init__.py b/SimPEG/forward/DCProblem/__init__.py deleted file mode 100644 index a868cf80..00000000 --- a/SimPEG/forward/DCProblem/__init__.py +++ /dev/null @@ -1,2 +0,0 @@ -from DCProblem import * -from DCutils import * diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py index fe849d41..5e326a59 100644 --- a/SimPEG/forward/__init__.py +++ b/SimPEG/forward/__init__.py @@ -1,2 +1,2 @@ from Problem import * -import DCProblem +from DCProblem import DCProblem From 9cd6e7ed08782584c9c8e6bed3177853e96ff254 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 23 Oct 2013 16:44:33 -0700 Subject: [PATCH 07/34] Issue #9 Ensure that h is a float in all meshes. --- SimPEG/mesh/Cyl1DMesh.py | 24 ++++++++++++------------ SimPEG/mesh/LogicallyOrthogonalMesh.py | 2 +- SimPEG/mesh/TensorMesh.py | 2 +- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/SimPEG/mesh/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py index 5307964b..e0d6529b 100644 --- a/SimPEG/mesh/Cyl1DMesh.py +++ b/SimPEG/mesh/Cyl1DMesh.py @@ -5,8 +5,8 @@ from SimPEG.utils import mkvc, ndgrid, sdiag class Cyl1DMesh(object): """ - Cyl1DMesh is a mesh class for cylindrically symmetric 1D problems - """ + Cyl1DMesh is a mesh class for cylindrically symmetric 1D problems + """ _meshType = 'CYL1D' @@ -20,7 +20,7 @@ class Cyl1DMesh(object): assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) # Ensure h contains 1D vectors - self._h = [mkvc(x) for x in h] + self._h = [mkvc(x.astype(float)) for x in h] if z0 is None: z0 = 0 @@ -146,7 +146,7 @@ class Cyl1DMesh(object): def vectorCCz(): doc = "Cell centered grid vector (1D) in the z direction" - fget = lambda self: self.hz.cumsum() - self.hz/2 + self._z0 + fget = lambda self: self.hz.cumsum() - self.hz/2 + self._z0 return locals() vectorCCz = property(**vectorCCz()) @@ -177,7 +177,7 @@ class Cyl1DMesh(object): self._gridFr = ndgrid([self.vectorNr, self.vectorCCz]) return self._gridFr return locals() - _gridFr = None + _gridFr = None gridFr = property(**gridFr()) def gridFz(): @@ -187,7 +187,7 @@ class Cyl1DMesh(object): self._gridFz = ndgrid([self.vectorCCr, self.vectorNz]) return self._gridFz return locals() - _gridFz = None + _gridFz = None gridFz = property(**gridFz()) #################################################### @@ -350,23 +350,23 @@ class Cyl1DMesh(object): np.all(loc[:,1]<=self.vectorNz.max()), \ "Points outside of mesh" - + if locType=='fz': Q = sp.lil_matrix((loc.shape[0], self.nF), dtype=float) for i, iloc in enumerate(loc): # Point is on a z-interface - if np.any(np.abs(self.vectorNz-iloc[1])<0.001): + if np.any(np.abs(self.vectorNz-iloc[1])<0.001): dFz = self.gridFz-iloc #Distance to z faces dFz[dFz[:,0]>0,:] = np.inf #Looking for next face to the left... indL = np.argmin(np.sum(dFz**2, axis=1)) #Closest one if self.gridFz[indL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation) - zFL = self.gridFz[indL,:] - zFLL = self.gridFz[indL-1,:] + zFL = self.gridFz[indL,:] + zFLL = self.gridFz[indL-1,:] Q[i, indL+self.nFr] = (iloc[0] - zFLL[0])/(zFL[0] - zFLL[0]) Q[i, indL+self.nFr-1] = -(iloc[0] - zFL[0])/(zFL[0] - zFLL[0]) else: - zFL = self.gridFz[indL,:] + zFL = self.gridFz[indL,:] zFR = self.gridFz[indL+1,:] Q[i,indL+self.nFr] = (zFR[0] - iloc[0])/(zFR[0] - zFL[0]) Q[i,indL+self.nFr+1] = (iloc[0] - zFL[0])/(zFR[0] - zFL[0]) @@ -400,7 +400,7 @@ class Cyl1DMesh(object): Q[i, indAL+self.nFr-1] = -(dzB/DZ)*(drL/DR) Q[i, indAL+self.nFr] = (dzB/DZ)*(drLL/DR) else: - indBR = indBL+1 # Face below and to the right + indBR = indBL+1 # Face below and to the right indAR = indAL + 1 # Face above and to the right zF_BR = self.gridFz[indBR,:] diff --git a/SimPEG/mesh/LogicallyOrthogonalMesh.py b/SimPEG/mesh/LogicallyOrthogonalMesh.py index 7b1e1bca..b510a754 100644 --- a/SimPEG/mesh/LogicallyOrthogonalMesh.py +++ b/SimPEG/mesh/LogicallyOrthogonalMesh.py @@ -38,7 +38,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView): # Save nodes to private variable _gridN as vectors self._gridN = np.ones((nodes[0].size, self.dim)) for i, node_i in enumerate(nodes): - self._gridN[:, i] = mkvc(node_i) + self._gridN[:, i] = mkvc(node_i.astype(float)) def gridCC(): doc = "Cell-centered grid." diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index a3329c14..d5123c98 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -38,7 +38,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) # Ensure h contains 1D vectors - self._h = [mkvc(x) for x in h] + self._h = [mkvc(x.astype(float)) for x in h] def __str__(self): outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) From 548ba9aa0a2d1c359aba28a99a06e539e7dfb065 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 23 Oct 2013 18:19:33 -0700 Subject: [PATCH 08/34] bug fixes to framework. --- SimPEG/forward/DCProblem.py | 67 +++++++++++++------------ SimPEG/inverse/Inversion.py | 45 +++++++++-------- SimPEG/inverse/Optimize.py | 4 +- SimPEG/mesh/DiffOperators.py | 62 +++++++++++++++++++++++ SimPEG/regularization/Regularization.py | 29 ++++++----- 5 files changed, 139 insertions(+), 68 deletions(-) diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index 5aad45bd..5a011bb1 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -3,7 +3,7 @@ from SimPEG.forward import Problem, SyntheticProblem from SimPEG.tests import checkDerivative from SimPEG.utils import ModelBuilder, sdiag import numpy as np -import scipy as sp +import scipy.sparse as sp import scipy.sparse.linalg as linalg class DCProblem(Problem): @@ -104,14 +104,44 @@ class DCProblem(Problem): return Jtv + +def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh): + """ Generate projection matrix (Q) and """ + elecend = 0.5+spacelec*(nelec-1) + elecLocR = np.linspace(elecini, elecend, nelec) + elecLocT = elecLocR+1 + nrx = nelec-1 + ntx = nelec-1 + q = np.zeros((mesh.nC, ntx)) + Q = np.zeros((mesh.nC, nrx)) + + for i in range(nrx): + + rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i])) + rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1])) + + txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i])) + txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1])) + + q[txind1,i] = 1 + q[txind2,i] = -1 + Q[rxind1,i] = 1 + Q[rxind2,i] = -1 + + Q = sp.csr_matrix(Q) + rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5 + return q, Q, rxmidLoc + + + if __name__ == '__main__': from SimPEG.regularization import Regularization from SimPEG import inverse # Create the mesh - h1 = np.ones(100) - h2 = np.ones(100) + h1 = np.ones(10) + h2 = np.ones(10) mesh = TensorMesh([h1,h2]) # Create some parameters for the model @@ -153,14 +183,14 @@ if __name__ == '__main__': problem = DCProblem(mesh) problem.P = P problem.RHS = q - problem.W = Wd problem.dobs = dobs + problem.std = dobs*0 + 0.05 m0 = mesh.gridCC[:,0]*0+sig1 # print problem.misfit(m0) # print problem.misfit(mSynth) - opt = inverse.InexactGaussNewton() + opt = inverse.InexactGaussNewton(maxIterLS=20) reg = Regularization(mesh) inv = inverse.Inversion(problem, reg, opt) @@ -181,30 +211,3 @@ if __name__ == '__main__': - -def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh): - """ Generate projection matrix (Q) and """ - elecend = 0.5+spacelec*(nelec-1) - elecLocR = np.linspace(elecini, elecend, nelec) - elecLocT = elecLocR+1 - nrx = nelec-1 - ntx = nelec-1 - q = np.zeros((mesh.nC, ntx)) - Q = np.zeros((mesh.nC, nrx)) - - for i in range(nrx): - - rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i])) - rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1])) - - txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i])) - txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1])) - - q[txind1,i] = 1 - q[txind2,i] = -1 - Q[rxind1,i] = 1 - Q[rxind2,i] = -1 - - Q = sp.csr_matrix(Q) - rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5 - return q, Q, rxmidLoc diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 8d92e15e..9fd10abd 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,5 +1,6 @@ import numpy as np -from SimPEG.utils import sdiag +import scipy.sparse as sp +from SimPEG.utils import sdiag, mkvc class Inversion(object): """docstring for Inversion""" @@ -16,7 +17,10 @@ class Inversion(object): """ Standard deviation weighting matrix. """ - return sdiag(1/self.prob.std) + if getattr(self,'_Wd',None) is None: + eps = np.linalg.norm(mkvc(self.prob.dobs),2)*1e-5 + self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps) + return self._Wd def run(self, m0): m = m0 @@ -41,7 +45,7 @@ class Inversion(object): u = self.prob.field(m) phi_d = self.dataObj(m, u) - phi_m = self.modelObj(m) + phi_m = self.reg.modelObj(m) self._phi_d_last = phi_d self._phi_m_last = phi_m @@ -50,27 +54,24 @@ class Inversion(object): out = (f,) if return_g: - phi_dDeriv = self.dataObjDeriv(m, u) - phi_mDeriv = self.modelObjDeriv(m) + phi_dDeriv = self.dataObjDeriv(m, u=u) + phi_mDeriv = self.reg.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 + phi_d2Deriv = self.dataObj2Deriv(m, v, u=u) + phi_m2Deriv = self.reg.modelObj2Deriv(m)*v return phi_d2Deriv + self._beta * phi_m2Deriv - out += (H_fun,) + operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=float ) + out += (operator,) return out - def modelObj(self, m, u=None): - self.reg.misfit(m) - - def dataObj(self, m, u=None): """ :param numpy.array m: geophysical model @@ -87,7 +88,7 @@ class Inversion(object): 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 = self.Wd*self.prob.misfit(m, u=u) R = mkvc(R) return 0.5*R.dot(R) @@ -123,17 +124,17 @@ class Inversion(object): """ if u is None: - u = self.field(m) + u = self.prob.field(m) - R = self.W*(self.dpred(m, u=u) - self.dobs) + R = self.Wd*self.prob.misfit(m, u=u) 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]) + 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]) return dmisfit - def dataObj2Deriv(self, m, u=None): + def dataObj2Deriv(self, m, v, u=None): """ :param numpy.array m: geophysical model :param numpy.array u: fields @@ -167,12 +168,12 @@ class Inversion(object): """ if u is None: - u = self.field(m) + u = self.prob.field(m) - R = self.W*(self.dpred(m, u=u) - self.dobs) + R = self.Wd*self.prob.misfit(m, u=u) 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]) + 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]) return dmisfit diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 43776d33..1bfda72c 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -2,6 +2,7 @@ import numpy as np import matplotlib.pyplot as plt from SimPEG.utils import mkvc, sdiag norm = np.linalg.norm +import scipy.sparse as sp class Minimize(object): @@ -149,7 +150,8 @@ class GaussNewton(Minimize): class InexactGaussNewton(Minimize): name = 'InexactGaussNewton' def findSearchDirection(self): - return sparse.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=10) + p, info = sp.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=10) + return p class SteepestDescent(Minimize): diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 110fe9bd..380e0d19 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -161,6 +161,68 @@ class DiffOperators(object): _cellGrad = None cellGrad = property(**cellGrad()) + def cellGradx(): + doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." + + def fget(self): + if getattr(self, '_cellGradx', None) is None: + BC = ['neumann', 'neumann'] + n = self.n + if(self.dim == 1): + G1 = ddxCellGrad(n[0], BC) + elif(self.dim == 2): + G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC)) + elif(self.dim == 3): + G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC)) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F','Fx', 'V') + V = self.vol + self._cellGradx = sdiag(S)*G1*sdiag(1/V) + return self._cellGradx + return locals() + cellGradx = property(**cellGradx()) + + + def cellGrady(): + doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." + def fget(self): + if self.dim < 2: + return None + if getattr(self, '_cellGrady', None) is None: + BC = ['neumann', 'neumann'] + n = self.n + if(self.dim == 2): + G2 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC)) + elif(self.dim == 3): + G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0])) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F','Fy', 'V') + V = self.vol + self._cellGrady = sdiag(S)*G2*sdiag(1/V) + return self._cellGrady + return locals() + cellGrady = property(**cellGrady()) + + + + def cellGradz(): + doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." + def fget(self): + if self.dim < 3: + return None + if getattr(self, '_cellGradz', None) is None: + BC = ['neumann', 'neumann'] + n = self.n + G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0])) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F','Fz', 'V') + V = self.vol + self._cellGradz = sdiag(S)*G3*sdiag(1/V) + return self._cellGradz + return locals() + cellGradz = property(**cellGradz()) + + def edgeCurl(): doc = "Construct the 3D curl operator." diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index f0239875..7e211df6 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/regularization/Regularization.py @@ -1,10 +1,13 @@ from SimPEG.utils import sdiag +import numpy as np class Regularization(object): """docstring for Regularization""" @property def mref(self): + if getattr(self, '_mref', None) is None: + self._mref = np.zeros(self.mesh.nC); return self._mref @mref.setter def mref(self, value): @@ -12,25 +15,25 @@ class Regularization(object): @property def Wx(self): - if self._Wx is None: - self._Wx = mesh.cellGradx + if getattr(self, '_Wx', None) is None: + self._Wx = self.mesh.cellGradx return self._Wx @property def Wy(self): - if self._Wy is None: - self._Wy = mesh.cellGrady + if getattr(self, '_Wy', None) is None: + self._Wy = self.mesh.cellGrady return self._Wy @property def Wz(self): - if self._Wz is None: - self._Wz = mesh.cellGradz + if getattr(self, '_Wz', None) is None: + self._Wz = self.mesh.cellGradz return self._Wz @property def Ws(self): - if self._Ws is None: + if getattr(self,'_Ws', None) is None: self._Ws = sdiag(self.mesh.vol) return self._Ws @@ -87,12 +90,12 @@ class Regularization(object): mobjDeriv = self.alpha_s * self.Ws.T * ( self.Ws * mresid) - mobjDeriv += self.alpha_x * self.Wx.T * ( self.Wx * mresid) + mobjDeriv = 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) + mobjDeriv = 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) + mobjDeriv = mobjDeriv + self.alpha_z * self.Wz.T * ( self.Wz * mresid) return mobjDeriv @@ -102,12 +105,12 @@ class Regularization(object): mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws - mobj2Deriv += self.alpha_x * self.Wx.T * self.Wx + mobj2Deriv = mobj2Deriv + self.alpha_x * self.Wx.T * self.Wx if self.mesh.dim > 1: - mobj2Deriv += self.alpha_y * self.Wy.T * self.Wy + mobj2Deriv = mobj2Deriv + self.alpha_y * self.Wy.T * self.Wy if self.mesh.dim > 2: - mobj2Deriv += self.alpha_z * self.Wz.T * self.Wz + mobj2Deriv = mobj2Deriv + self.alpha_z * self.Wz.T * self.Wz return mobj2Deriv From d94d3f7f374a69faa5752165f86b3a2e50c13bbc Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 23 Oct 2013 18:35:37 -0700 Subject: [PATCH 09/34] =?UTF-8?q?fixes=20to=20inversion=20framework.=20DC?= =?UTF-8?q?=20problem=20is=20kinda=20working=E2=80=A6=3F!?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- SimPEG/forward/DCProblem.py | 24 +++++++++++++----------- SimPEG/inverse/Inversion.py | 14 ++++++++++++-- 2 files changed, 25 insertions(+), 13 deletions(-) diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index 5a011bb1..c2494daf 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -140,8 +140,8 @@ if __name__ == '__main__': from SimPEG import inverse # Create the mesh - h1 = np.ones(10) - h2 = np.ones(10) + h1 = np.ones(100) + h2 = np.ones(100) mesh = TensorMesh([h1,h2]) # Create some parameters for the model @@ -177,7 +177,7 @@ if __name__ == '__main__': dobs, Wd = synthetic.createData(mSynth, std=0.05) u = synthetic.field(mSynth) - mesh.plotImage(u[:,10], showIt=False) + # mesh.plotImage(u[:,10], showIt=False) # Now set up the problem to do some minimization problem = DCProblem(mesh) @@ -190,23 +190,25 @@ if __name__ == '__main__': # print problem.misfit(m0) # print problem.misfit(mSynth) - opt = inverse.InexactGaussNewton(maxIterLS=20) + opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=1) reg = Regularization(mesh) inv = inverse.Inversion(problem, reg, opt) - inv.run(m0) + m = inv.run(m0) + + mesh.plotImage(m,showIt=True) # Check Derivative - derChk = lambda m: [problem.misfit(m), problem.misfitDeriv(m)] + 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)) + # 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)) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 9fd10abd..f217be83 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -22,6 +22,15 @@ class Inversion(object): self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps) return self._Wd + @property + def phi_d_target(self): + if getattr(self, '_phi_d_target', None) is None: + return self.prob.dobs.size + return self._phi_d_target + @phi_d_target.setter + def phi_d_target(self, value): + self._phi_d_target = value + def run(self, m0): m = m0 self._iter = 0 @@ -30,13 +39,14 @@ class Inversion(object): m = self.opt.minimize(self.evalFunction,m) if self.stoppingCriteria(): break self._iter += 1 + return m def getBeta(self): - return 1 + return 1e3 def stoppingCriteria(self): self._STOP = np.zeros(2,dtype=bool) - self._STOP[0] = self._iter >= maxIter + self._STOP[0] = self._iter >= self.maxIter self._STOP[1] = self._phi_d_last <= self.phi_d_target return np.any(self._STOP) From f4c63d47a3ae58aae1a5a3520b16e86e8cc71de6 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 24 Oct 2013 14:12:37 -0700 Subject: [PATCH 10/34] 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 From fc4294eb4dc6fe3f4f43c1c413ea5ccfc5011bcc Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 24 Oct 2013 15:11:48 -0700 Subject: [PATCH 11/34] Added A sample Linear problem that runs! --- SimPEG/forward/LinearProblem.py | 91 +++++++++++++++++++++++++ SimPEG/forward/Problem.py | 22 +++--- SimPEG/inverse/Inversion.py | 11 ++- SimPEG/inverse/Optimize.py | 2 +- SimPEG/regularization/Regularization.py | 20 +++--- 5 files changed, 124 insertions(+), 22 deletions(-) create mode 100644 SimPEG/forward/LinearProblem.py diff --git a/SimPEG/forward/LinearProblem.py b/SimPEG/forward/LinearProblem.py new file mode 100644 index 00000000..d8f3621e --- /dev/null +++ b/SimPEG/forward/LinearProblem.py @@ -0,0 +1,91 @@ +import numpy as np +from SimPEG.mesh import TensorMesh +from SimPEG.forward import Problem +from SimPEG.regularization import Regularization +from SimPEG.inverse import * +import matplotlib.pyplot as plt + +N = 100 +h = np.ones(N)/N +M = TensorMesh([h]) + +nk = 20 +jk = np.linspace(1.,20.,nk) +p = -0.25 +q = 0.25 + + + +g = lambda k: np.exp(p*jk[k]*M.vectorCCx)*np.cos(2*np.pi*q*jk[k]*M.vectorCCx) + +G = np.empty((nk, M.nC)) + +for i in range(nk): + G[i,:] = g(i) + + + +plt.figure(1) +for i in range(nk): + plt.plot(G[i,:]) + + +m_true = np.zeros(M.nC) +m_true[M.vectorCCx > 0.3] = 1. +m_true[M.vectorCCx > 0.45] = -0.5 +m_true[M.vectorCCx > 0.6] = 0 + + +d_true = G.dot(m_true) +noise = 0.1 * np.random.rand(d_true.size) + +d_obs = d_true + noise + +# plt.figure(3) +# plt.plot(d_true,'-o') +# plt.plot(d_obs,'r-o') + + + +class LinearProblem(Problem): + """docstring for LinearProblem""" + + def dpred(self, m, u=None): + return self.G.dot(m) + + def J(self, m, v, u=None): + return G.dot(v) + + def Jt(self, m, v, u=None): + return G.T.dot(v) + + def modelTransform(self, m): + return m + + def modelTransformDeriv(self, m): + return sp.eye(m.size) + +prob = LinearProblem(M) +prob.G = G +prob.dobs = d_obs +prob.std = np.ones_like(d_obs)*0.1 + +reg = Regularization(M) + +opt = InexactGaussNewton(maxIter=10) + +inv = Inversion(prob,reg,opt) + +m0 = np.zeros_like(m_true) + +mrec = inv.run(m0) + + +plt.figure(2) + +plt.plot(M.vectorCCx, m_true, 'b-') +plt.plot(M.vectorCCx, mrec, 'r-') + + + +plt.show() diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 0333b8b4..32dad675 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -82,6 +82,17 @@ class Problem(object): def dobs(self, value): self._dobs = value + def dpred(self, m, u=None): + """ + Predicted data. + + .. math:: + d_\\text{pred} = Pu(m) + """ + if u is None: + u = self.field(m) + return self.P*u + def misfit(self, m, u=None): """ :param numpy.array m: geophysical model @@ -152,17 +163,6 @@ class Problem(object): """ pass - def dpred(self, m, u=None): - """ - Predicted data. - - .. math:: - d_\\text{pred} = Pu(m) - """ - if u is None: - u = self.field(m) - return self.P*u - def modelTransform(self, m): """ :param numpy.array m: model diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index fd166318..cf96676f 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -25,8 +25,15 @@ class Inversion(object): @property def phi_d_target(self): + """ + target for phi_d + + By default this is the number of data. + + Note that we do not set the target if it is None, but we return the default value. + """ if getattr(self, '_phi_d_target', None) is None: - return self.prob.dobs.size + return self.prob.dobs.size # return self._phi_d_target @phi_d_target.setter def phi_d_target(self, value): @@ -43,7 +50,7 @@ class Inversion(object): return m def getBeta(self): - return 1e2 + return 1e-2 def stoppingCriteria(self): self._STOP = np.zeros(2,dtype=bool) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 2545a381..4ac96d6a 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -22,7 +22,7 @@ class Minimize(object): tolF = 1e-4 tolX = 1e-4 tolG = 1e-4 - eps = 1e-16 + eps = 1e-5 printIter = [] # push to here if you want to print these on iter diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index 7e211df6..6f5970e6 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/regularization/Regularization.py @@ -13,29 +13,33 @@ class Regularization(object): def mref(self, value): self._mref = value + @property + def Ws(self): + if getattr(self,'_Ws', None) is None: + self._Ws = sdiag(self.mesh.vol) + return self._Ws + @property def Wx(self): if getattr(self, '_Wx', None) is None: - self._Wx = self.mesh.cellGradx + a = self.mesh.r(self.mesh.area,'F','Fx','V') + self._Wx = sdiag(a)*self.mesh.cellGradx return self._Wx @property def Wy(self): if getattr(self, '_Wy', None) is None: - self._Wy = self.mesh.cellGrady + a = self.mesh.r(self.mesh.area,'F','Fy','V') + self._Wy = sdiag(a)*self.mesh.cellGrady return self._Wy @property def Wz(self): if getattr(self, '_Wz', None) is None: - self._Wz = self.mesh.cellGradz + a = self.mesh.r(self.mesh.area,'F','Fz','V') + self._Wz = sdiag(a)*self.mesh.cellGradz return self._Wz - @property - def Ws(self): - if getattr(self,'_Ws', None) is None: - self._Ws = sdiag(self.mesh.vol) - return self._Ws def __init__(self, mesh): From d974843d9871b66306f03b3b06b6838c21e9646b Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 24 Oct 2013 15:33:07 -0700 Subject: [PATCH 12/34] Moved model transforms to different file. --- SimPEG/forward/DCProblem.py | 4 +-- SimPEG/forward/ModelTransforms.py | 49 +++++++++++++++++++++++++++++++ SimPEG/forward/Problem.py | 26 ++++------------ SimPEG/forward/__init__.py | 1 + 4 files changed, 57 insertions(+), 23 deletions(-) create mode 100644 SimPEG/forward/ModelTransforms.py diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index 48f25f12..a2ab3967 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -1,5 +1,5 @@ from SimPEG.mesh import TensorMesh -from SimPEG.forward import Problem, SyntheticProblem +from SimPEG.forward import Problem, SyntheticProblem, ModelTransforms from SimPEG.tests import checkDerivative from SimPEG.utils import ModelBuilder, sdiag, mkvc from SimPEG import Solver @@ -7,7 +7,7 @@ import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as linalg -class DCProblem(Problem): +class DCProblem(Problem, ModelTransforms.LogModel): """ **DCProblem** diff --git a/SimPEG/forward/ModelTransforms.py b/SimPEG/forward/ModelTransforms.py new file mode 100644 index 00000000..8b7dc65a --- /dev/null +++ b/SimPEG/forward/ModelTransforms.py @@ -0,0 +1,49 @@ +import numpy as np + + +class LogModel(object): + """docstring for LogModel""" + def modelTransform(self, m): + """ + :param numpy.array m: model + :rtype: numpy.array + :return: transformed model + + The modelTransform changes the model into the physical property. + + A common example of this is to invert for electrical conductivity + in log space. In this case, your model will be log(sigma) and to + get back to sigma, you can take the exponential: + + .. math:: + + m = \log{\sigma} + + \exp{m} = \exp{\log{\sigma}} = \sigma + """ + return np.exp(mkvc(m)) + + def modelTransformDeriv(self, m): + """ + :param numpy.array m: model + :rtype: scipy.csr_matrix + :return: derivative of transformed model + + The modelTransform changes the model into the physical property. + The modelTransformDeriv provides the derivative of the modelTransform. + + If the model transform is: + + .. math:: + + m = \log{\sigma} + + \exp{m} = \exp{\log{\sigma}} = \sigma + + Then the derivative is: + + .. math:: + + \\frac{\partial \exp{m}}{\partial m} = \\text{sdiag}(\exp{m}) + """ + return sdiag(np.exp(mkvc(m))) diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 32dad675..679e8686 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -175,13 +175,8 @@ class Problem(object): in log space. In this case, your model will be log(sigma) and to get back to sigma, you can take the exponential: - .. math:: - - m = \log{\sigma} - - \exp{m} = \exp{\log{\sigma}} = \sigma """ - return np.exp(mkvc(m)) + return m def modelTransformDeriv(self, m): """ @@ -191,22 +186,8 @@ class Problem(object): The modelTransform changes the model into the physical property. The modelTransformDeriv provides the derivative of the modelTransform. - - If the model transform is: - - .. math:: - - m = \log{\sigma} - - \exp{m} = \exp{\log{\sigma}} = \sigma - - Then the derivative is: - - .. math:: - - \\frac{\partial \exp{m}}{\partial m} = \\text{sdiag}(\exp{m}) """ - return sdiag(np.exp(mkvc(m))) + return sp.eye(m.size) @@ -239,3 +220,6 @@ class SyntheticProblem(object): eps = np.linalg.norm(mkvc(dobs),2)*1e-5 Wd = 1/(abs(dobs)*std+eps) return dobs, Wd + + + diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py index 5e326a59..5f9dafaf 100644 --- a/SimPEG/forward/__init__.py +++ b/SimPEG/forward/__init__.py @@ -1,2 +1,3 @@ from Problem import * from DCProblem import DCProblem +import ModelTransforms From 3bbe258eabab0d5f6022dd065c7841407d9a0fed Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 25 Oct 2013 11:14:21 -0700 Subject: [PATCH 13/34] beta schedualing --- SimPEG/inverse/BetaSchedual.py | 12 ++++++++++++ SimPEG/inverse/Inversion.py | 12 ++++++++++-- SimPEG/inverse/Optimize.py | 1 + SimPEG/inverse/__init__.py | 1 + 4 files changed, 24 insertions(+), 2 deletions(-) create mode 100644 SimPEG/inverse/BetaSchedual.py diff --git a/SimPEG/inverse/BetaSchedual.py b/SimPEG/inverse/BetaSchedual.py new file mode 100644 index 00000000..62fe8416 --- /dev/null +++ b/SimPEG/inverse/BetaSchedual.py @@ -0,0 +1,12 @@ + + +class Cooling(object): + """Simple Beta Schedual""" + + beta0 = 1.e6 + beta_coolingFactor = 5. + + def getBeta(self): + if self._beta is None: + return beta0 + return self._beta * beta_coolingFactor diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index cf96676f..b30b666e 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,8 +1,9 @@ import numpy as np import scipy.sparse as sp from SimPEG.utils import sdiag, mkvc +from SimPEG.inverse import BS -class Inversion(object): +class Inversion(object, BetaSchedual.Cooling): """docstring for Inversion""" maxIter = 10 @@ -42,6 +43,7 @@ class Inversion(object): def run(self, m0): m = m0 self._iter = 0 + self._beta = None while True: self._beta = self.getBeta() m = self.opt.minimize(self.evalFunction,m) @@ -49,8 +51,13 @@ class Inversion(object): self._iter += 1 return m + beta0 = 1.e6 + beta_coolingFactor = 5. + def getBeta(self): - return 1e-2 + if self._beta is None: + return beta0 + return self._beta * beta_coolingFactor def stoppingCriteria(self): self._STOP = np.zeros(2,dtype=bool) @@ -191,3 +198,4 @@ class Inversion(object): 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 4ac96d6a..ac0045c4 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -158,6 +158,7 @@ class GaussNewton(Minimize): class InexactGaussNewton(Minimize): name = 'InexactGaussNewton' def findSearchDirection(self): + # TODO: use BFGS as a preconditioner or gauss sidel of the WtW or solve WtW directly p, info = sp.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=5) return p diff --git a/SimPEG/inverse/__init__.py b/SimPEG/inverse/__init__.py index c4ca335b..2ca95fb5 100644 --- a/SimPEG/inverse/__init__.py +++ b/SimPEG/inverse/__init__.py @@ -1,2 +1,3 @@ from Optimize import * from Inversion import * +import BetaSchedual From dbd24b71c1276cb590570fdb0ed567d4c28d4259 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 26 Oct 2013 14:16:50 -0700 Subject: [PATCH 14/34] Added kwargs to the inversion. can use these to initiate your inversion object. --- SimPEG/forward/LinearProblem.py | 10 ++-------- SimPEG/inverse/Inversion.py | 20 ++++++++++++++------ SimPEG/inverse/Optimize.py | 6 +++--- 3 files changed, 19 insertions(+), 17 deletions(-) diff --git a/SimPEG/forward/LinearProblem.py b/SimPEG/forward/LinearProblem.py index d8f3621e..2260f8db 100644 --- a/SimPEG/forward/LinearProblem.py +++ b/SimPEG/forward/LinearProblem.py @@ -59,12 +59,6 @@ class LinearProblem(Problem): def Jt(self, m, v, u=None): return G.T.dot(v) - def modelTransform(self, m): - return m - - def modelTransformDeriv(self, m): - return sp.eye(m.size) - prob = LinearProblem(M) prob.G = G prob.dobs = d_obs @@ -72,9 +66,9 @@ prob.std = np.ones_like(d_obs)*0.1 reg = Regularization(M) -opt = InexactGaussNewton(maxIter=10) +opt = InexactGaussNewton(maxIter=20) -inv = Inversion(prob,reg,opt) +inv = Inversion(prob,reg,opt,beta0=1e-4) m0 = np.zeros_like(m_true) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index b30b666e..02a81daa 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,18 +1,26 @@ import numpy as np import scipy.sparse as sp from SimPEG.utils import sdiag, mkvc -from SimPEG.inverse import BS -class Inversion(object, BetaSchedual.Cooling): +class Inversion(object): """docstring for Inversion""" maxIter = 10 - def __init__(self, prob, reg, opt): + def __init__(self, prob, reg, opt, **kwargs): self.prob = prob self.reg = reg self.opt = opt self.opt.parent = self + self.setKwargs(**kwargs) + + def setKwargs(self, **kwargs): + # Set the variables, throw an error if they don't exist. + for attr in kwargs: + if hasattr(self, attr): + setattr(self, attr, kwargs[attr]) + else: + raise Exception('%s attr is not recognized' % attr) @property def Wd(self): @@ -51,13 +59,13 @@ class Inversion(object, BetaSchedual.Cooling): self._iter += 1 return m - beta0 = 1.e6 + beta0 = 1.e2 beta_coolingFactor = 5. def getBeta(self): if self._beta is None: - return beta0 - return self._beta * beta_coolingFactor + return self.beta0 + return self._beta * self.beta_coolingFactor def stoppingCriteria(self): self._STOP = np.zeros(2,dtype=bool) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index ac0045c4..9e3b5463 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -19,9 +19,9 @@ class Minimize(object): maxIterLS = 10 LSreduction = 1e-4 LSshorten = 0.5 - tolF = 1e-4 - tolX = 1e-4 - tolG = 1e-4 + tolF = 1e-1 + tolX = 1e-1 + tolG = 1e-1 eps = 1e-5 printIter = [] # push to here if you want to print these on iter From 5cbc980b9c283bf2fbaef1c29245aa1a6786fc25 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 26 Oct 2013 15:55:42 -0700 Subject: [PATCH 15/34] Mixins work differently than I thought, this fixes the DCProblem code. --- SimPEG/forward/DCProblem.py | 2 +- SimPEG/forward/ModelTransforms.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index a2ab3967..b8e8248d 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -7,7 +7,7 @@ import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as linalg -class DCProblem(Problem, ModelTransforms.LogModel): +class DCProblem(ModelTransforms.LogModel, Problem): """ **DCProblem** diff --git a/SimPEG/forward/ModelTransforms.py b/SimPEG/forward/ModelTransforms.py index 8b7dc65a..ea89b974 100644 --- a/SimPEG/forward/ModelTransforms.py +++ b/SimPEG/forward/ModelTransforms.py @@ -1,5 +1,5 @@ import numpy as np - +from SimPEG.utils import mkvc, sdiag class LogModel(object): """docstring for LogModel""" From 1555caa7758aa46d5771a3221391c171fbfd645e Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 26 Oct 2013 16:02:04 -0700 Subject: [PATCH 16/34] Changed to SimPEG solver for DC problem. --- SimPEG/forward/DCProblem.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index b8e8248d..c2a615ed 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -40,13 +40,8 @@ class DCProblem(ModelTransforms.LogModel, Problem): def field(self, m): A = self.createMatrix(m) - solve = linalg.factorized(A) - - nRHSs = self.RHS.shape[1] # Number of RHSs - phi = np.zeros((self.mesh.nC, nRHSs)) + np.nan - for ii in range(nRHSs): - phi[:,ii] = solve(self.RHS[:,ii]) - + solve = Solver(A) + phi = solve.solve(self.RHS) return phi def J(self, m, v, u=None): From 49d08f3245d3df883d18902a3cd5f15bb849bd83 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 26 Oct 2013 16:05:17 -0700 Subject: [PATCH 17/34] add factorize option to SimPEG.Solver --- SimPEG/Solver.py | 44 ++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 40 insertions(+), 4 deletions(-) diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index 2087ef5f..cbd4a3ea 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -36,21 +36,27 @@ class Solver(object): del self.dsolve self.dsolve = None - def solveDirect(self, b, backend='scipy'): + def solveDirect(self, b, factorize=False, backend='scipy'): assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch' - if self.dsolve is None: + if factorize and self.dsolve is None: self.A = self.A.tocsc() # for efficiency self.dsolve = linalg.factorized(self.A) if len(b.shape) == 1 or b.shape[1] == 1: # Just one RHS - return self.dsolve(b) + if factorize: + return self.dsolve(b) + else: + return linalg.dsolve.spsolve(self.A, b) # Multiple RHSs X = np.empty_like(b) for i in range(b.shape[1]): - X[:,i] = self.dsolve(b[:,i]) + if factorize: + X[:,i] = self.dsolve(b[:,i]) + else: + X[:,i] = linalg.dsolve.spsolve(self.A,b[:,i]) return X @@ -73,3 +79,33 @@ class Solver(object): for i in range(b.shape[1]): X[:,i] = b[:,i]/diagA return X + + +if __name__ == '__main__': + from SimPEG.mesh import TensorMesh + from time import time + h1 = np.ones(20)*100. + h2 = np.ones(20)*100. + h3 = np.ones(20)*100. + + h = [h1,h2,h3] + + M = TensorMesh(h) + + D = M.faceDiv + G = M.cellGrad + Msig = M.getFaceMass() + A = D*Msig*G + + rhs = np.random.rand(M.nC) + + + tic = time() + solve = Solver(A, options={'factorize':True}) + x = solve.solve(rhs) + print 'Factorized', time() - tic + tic = time() + solve = Solver(A, options={'factorize':False}) + x = solve.solve(rhs) + print 'spsolve', time() - tic + From 86c1080631d2213ac8402a62fd131b5ae91bf2e3 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sat, 26 Oct 2013 16:19:02 -0700 Subject: [PATCH 18/34] Documentation updates. --- SimPEG/forward/LinearProblem.py | 126 +++++++++--------- SimPEG/forward/__init__.py | 1 + .../{BetaSchedual.py => BetaSchedule.py} | 2 +- SimPEG/inverse/__init__.py | 2 +- docs/api_Optimize.rst | 15 +++ docs/api_Problem.rst | 10 +- 6 files changed, 89 insertions(+), 67 deletions(-) rename SimPEG/inverse/{BetaSchedual.py => BetaSchedule.py} (86%) diff --git a/SimPEG/forward/LinearProblem.py b/SimPEG/forward/LinearProblem.py index 2260f8db..d30a5b4d 100644 --- a/SimPEG/forward/LinearProblem.py +++ b/SimPEG/forward/LinearProblem.py @@ -5,47 +5,6 @@ from SimPEG.regularization import Regularization from SimPEG.inverse import * import matplotlib.pyplot as plt -N = 100 -h = np.ones(N)/N -M = TensorMesh([h]) - -nk = 20 -jk = np.linspace(1.,20.,nk) -p = -0.25 -q = 0.25 - - - -g = lambda k: np.exp(p*jk[k]*M.vectorCCx)*np.cos(2*np.pi*q*jk[k]*M.vectorCCx) - -G = np.empty((nk, M.nC)) - -for i in range(nk): - G[i,:] = g(i) - - - -plt.figure(1) -for i in range(nk): - plt.plot(G[i,:]) - - -m_true = np.zeros(M.nC) -m_true[M.vectorCCx > 0.3] = 1. -m_true[M.vectorCCx > 0.45] = -0.5 -m_true[M.vectorCCx > 0.6] = 0 - - -d_true = G.dot(m_true) -noise = 0.1 * np.random.rand(d_true.size) - -d_obs = d_true + noise - -# plt.figure(3) -# plt.plot(d_true,'-o') -# plt.plot(d_obs,'r-o') - - class LinearProblem(Problem): """docstring for LinearProblem""" @@ -59,27 +18,72 @@ class LinearProblem(Problem): def Jt(self, m, v, u=None): return G.T.dot(v) -prob = LinearProblem(M) -prob.G = G -prob.dobs = d_obs -prob.std = np.ones_like(d_obs)*0.1 +if __name__ == '__main__': + N = 100 + h = np.ones(N)/N + M = TensorMesh([h]) -reg = Regularization(M) - -opt = InexactGaussNewton(maxIter=20) - -inv = Inversion(prob,reg,opt,beta0=1e-4) - -m0 = np.zeros_like(m_true) - -mrec = inv.run(m0) - - -plt.figure(2) - -plt.plot(M.vectorCCx, m_true, 'b-') -plt.plot(M.vectorCCx, mrec, 'r-') + nk = 20 + jk = np.linspace(1.,20.,nk) + p = -0.25 + q = 0.25 -plt.show() + g = lambda k: np.exp(p*jk[k]*M.vectorCCx)*np.cos(2*np.pi*q*jk[k]*M.vectorCCx) + + G = np.empty((nk, M.nC)) + + for i in range(nk): + G[i,:] = g(i) + + + + plt.figure(1) + for i in range(nk): + plt.plot(G[i,:]) + + + m_true = np.zeros(M.nC) + m_true[M.vectorCCx > 0.3] = 1. + m_true[M.vectorCCx > 0.45] = -0.5 + m_true[M.vectorCCx > 0.6] = 0 + + + d_true = G.dot(m_true) + noise = 0.1 * np.random.rand(d_true.size) + + d_obs = d_true + noise + + # plt.figure(3) + # plt.plot(d_true,'-o') + # plt.plot(d_obs,'r-o') + + + + + + prob = LinearProblem(M) + prob.G = G + prob.dobs = d_obs + prob.std = np.ones_like(d_obs)*0.1 + + reg = Regularization(M) + + opt = InexactGaussNewton(maxIter=20) + + inv = Inversion(prob,reg,opt,beta0=1e-4) + + m0 = np.zeros_like(m_true) + + mrec = inv.run(m0) + + + plt.figure(2) + + plt.plot(M.vectorCCx, m_true, 'b-') + plt.plot(M.vectorCCx, mrec, 'r-') + + + + plt.show() diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py index 5f9dafaf..a7c58967 100644 --- a/SimPEG/forward/__init__.py +++ b/SimPEG/forward/__init__.py @@ -1,3 +1,4 @@ from Problem import * from DCProblem import DCProblem +from LinearProblem import LinearProblem import ModelTransforms diff --git a/SimPEG/inverse/BetaSchedual.py b/SimPEG/inverse/BetaSchedule.py similarity index 86% rename from SimPEG/inverse/BetaSchedual.py rename to SimPEG/inverse/BetaSchedule.py index 62fe8416..2c633a6a 100644 --- a/SimPEG/inverse/BetaSchedual.py +++ b/SimPEG/inverse/BetaSchedule.py @@ -1,7 +1,7 @@ class Cooling(object): - """Simple Beta Schedual""" + """Simple Beta Schedule""" beta0 = 1.e6 beta_coolingFactor = 5. diff --git a/SimPEG/inverse/__init__.py b/SimPEG/inverse/__init__.py index 2ca95fb5..14bddce7 100644 --- a/SimPEG/inverse/__init__.py +++ b/SimPEG/inverse/__init__.py @@ -1,3 +1,3 @@ from Optimize import * from Inversion import * -import BetaSchedual +import BetaSchedule diff --git a/docs/api_Optimize.rst b/docs/api_Optimize.rst index 9765bd3b..d726b5f8 100644 --- a/docs/api_Optimize.rst +++ b/docs/api_Optimize.rst @@ -6,3 +6,18 @@ Optimize .. automodule:: SimPEG.inverse.Optimize :members: :undoc-members: + + +Inversion +********* + +.. automodule:: SimPEG.inverse.Inversion + :members: + :undoc-members: + +Beta Schedule +************* + +.. automodule:: SimPEG.inverse.BetaSchedule + :members: + :undoc-members: diff --git a/docs/api_Problem.rst b/docs/api_Problem.rst index 83250f3e..068e30da 100644 --- a/docs/api_Problem.rst +++ b/docs/api_Problem.rst @@ -13,14 +13,16 @@ Problem DCProblem ********* -.. automodule:: SimPEG.forward.DCProblem.DCProblem +.. automodule:: SimPEG.forward.DCProblem :members: :undoc-members: -DCutils -******* -.. automodule:: SimPEG.forward.DCProblem.DCutils +Linear Problem +************** + +.. automodule:: SimPEG.forward.LinearProblem :members: :undoc-members: + From c33c9bee55e65485f36b8c116ff38369342bd9cb Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 30 Oct 2013 22:42:25 -0600 Subject: [PATCH 19/34] Added forward and backwards solvers implemented in python. Added tests for direct solvers. --- SimPEG/Solver.py | 39 +++++++++++-- SimPEG/tests/test_Solver.py | 111 ++++++++++++++++++++++++++++++++++++ 2 files changed, 146 insertions(+), 4 deletions(-) create mode 100644 SimPEG/tests/test_Solver.py diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index cbd4a3ea..a9105c84 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -1,4 +1,5 @@ import numpy as np +import scipy.sparse as sparse import scipy.sparse.linalg as linalg @@ -64,10 +65,36 @@ class Solver(object): pass def solveBackward(self, b): - pass + "Perform a backwards solve with upper triangular A in CSR format." + if type(self.A) is not sparse.csr.csr_matrix: + from scipy.sparse import csr_matrix + self.A = csr_matrix(self.A) + vals = self.A.data + rowptr = self.A.indptr + colind = self.A.indices + x = np.empty_like(b) # empty() is faster than zeros(). + for i in reversed(xrange(self.A.shape[0])): + ith_row = vals[rowptr[i] : rowptr[i+1]] + cols = colind[rowptr[i] : rowptr[i+1]] + x_vals = x[cols] + x[i] = (b[i] - np.dot(ith_row[1:], x_vals[1:])) / ith_row[0] + return x def solveForward(self, b): - pass + "Perform a forward solve with lower triangular A in CSR format." + if type(self.A) is not sparse.csr.csr_matrix: + from scipy.sparse import csr_matrix + self.A = csr_matrix(self.A) + vals = self.A.data + rowptr = self.A.indptr + colind = self.A.indices + x = np.empty_like(b) # empty() is faster than zeros(). + for i in xrange(self.A.shape[0]): + ith_row = vals[rowptr[i] : rowptr[i+1]] + cols = colind[rowptr[i] : rowptr[i+1]] + x_vals = x[cols] + x[i] = (b[i] - np.dot(ith_row[:-1], x_vals[:-1])) / ith_row[-1] + return x def solveDiagonal(self, b): diagA = self.A.diagonal() @@ -96,16 +123,20 @@ if __name__ == '__main__': G = M.cellGrad Msig = M.getFaceMass() A = D*Msig*G + A[0,0] *= 10 # remove the constant null space from the matrix - rhs = np.random.rand(M.nC) - + e = np.ones(M.nC) + rhs = A.dot(e) tic = time() solve = Solver(A, options={'factorize':True}) x = solve.solve(rhs) print 'Factorized', time() - tic + print np.linalg.norm(e-x,np.inf) tic = time() solve = Solver(A, options={'factorize':False}) x = solve.solve(rhs) print 'spsolve', time() - tic + print np.linalg.norm(e-x,np.inf) + diff --git a/SimPEG/tests/test_Solver.py b/SimPEG/tests/test_Solver.py new file mode 100644 index 00000000..9b5cc0e7 --- /dev/null +++ b/SimPEG/tests/test_Solver.py @@ -0,0 +1,111 @@ +import unittest +from SimPEG import Solver +from SimPEG.mesh import TensorMesh +from SimPEG.utils import sdiag +import numpy as np +import scipy.sparse as sparse + +TOL = 1e-10 +numRHS = 5 + + +class TestSolver(unittest.TestCase): + + def setUp(self): + h1 = np.ones(10)*100. + h2 = np.ones(10)*100. + h3 = np.ones(10)*100. + + h = [h1,h2,h3] + + M = TensorMesh(h) + + D = M.faceDiv + G = M.cellGrad + Msig = M.getFaceMass() + A = D*Msig*G + A[0,0] *= 10 # remove the constant null space from the matrix + + self.A = A + self.M = M + + def test_directFactored_1(self): + solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':True,'backend':'scipy'}) + e = np.ones(self.M.nC) + rhs = self.A.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + + def test_directFactored_M(self): + solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':True,'backend':'scipy'}) + e = np.ones((self.M.nC,numRHS)) + rhs = self.A.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directSpsolve_1(self): + solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':False,'backend':'scipy'}) + e = np.ones(self.M.nC) + rhs = self.A.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directSpsolve_M(self): + solve = Solver(self.A, doDirect=True, flag=None, options={'factorize':False,'backend':'scipy'}) + e = np.ones((self.M.nC, numRHS)) + rhs = self.A.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directLower_1(self): + AL = sparse.tril(self.A) + solve = Solver(AL, doDirect=True, flag='L', options={}) + e = np.ones(self.M.nC) + rhs = AL.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directLower_M(self): + AL = sparse.tril(self.A) + solve = Solver(AL, doDirect=True, flag='L', options={}) + e = np.ones((self.M.nC,numRHS)) + rhs = AL.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directUpper_1(self): + AU = sparse.triu(self.A) + solve = Solver(AU, doDirect=True, flag='U', options={}) + e = np.ones(self.M.nC) + rhs = AU.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directUpper_M(self): + AU = sparse.triu(self.A) + solve = Solver(AU, doDirect=True, flag='U', options={}) + e = np.ones((self.M.nC,numRHS)) + rhs = AU.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directDiagonal_1(self): + AD = sdiag(self.A.diagonal()) + solve = Solver(AD, doDirect=True, flag='D', options={}) + e = np.ones(self.M.nC) + rhs = AD.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + def test_directDiagonal_M(self): + AD = sdiag(self.A.diagonal()) + solve = Solver(AD, doDirect=True, flag='D', options={}) + e = np.ones((self.M.nC,numRHS)) + rhs = AD.dot(e) + x = solve.solve(rhs) + self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True) + + +if __name__ == '__main__': + unittest.main() From 3f843b28439849fe29f0318ce26fd4de53fcd2d1 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 30 Oct 2013 23:00:42 -0600 Subject: [PATCH 20/34] Added documentation for SimPEG.Solver --- SimPEG/Solver.py | 79 +++++++++++++++++++++++++++++++++++++++++---- docs/api_Solver.rst | 9 ++++++ docs/index.rst | 1 + 3 files changed, 82 insertions(+), 7 deletions(-) create mode 100644 docs/api_Solver.rst diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index a9105c84..127a8c9f 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -4,9 +4,30 @@ import scipy.sparse.linalg as linalg class Solver(object): - """docstring for Solver""" - def __init__(self, A, doDirect=True, flag=None, options={}): + """ + Solver is a light wrapper on the various types of + linear solvers available in python. + :param scipy.sparse A: Matrix + :param bool doDirect: if you want a direct solver + :param string flag: Matrix type flag for special solves: [None, 'L', 'U', 'D'] + :param dict options: options which are passed to each sub solver, see each for details. + :rtype: Solver + :return: Solver + + To use for direct solvers:: + + solve = Solver(A, doDirect=True, flag=None, options={'factorize':True,'backend':'scipy'}) + x = solve.solve(rhs) + + Or in one line:: + + x = Solver(A).solve(rhs) + + The flag can be set to None, 'L', 'U', or 'D', for general, lower, upper, and diagonal matrices, respectively. + + """ + def __init__(self, A, doDirect=True, flag=None, options={}): assert type(doDirect) is bool, 'doDirect must be a boolean' assert flag in [None, 'L', 'U', 'D'], "flag must be set to None, 'L', 'U', or 'D'" @@ -18,6 +39,17 @@ class Solver(object): self.options = options def solve(self, b): + """ + Solves the linear system. + + .. math:: + + Ax=b + + :param numpy.ndarray b: the right hand side + :rtype: numpy.ndarray + :return: x + """ if self.flag is None and self.doDirect: return self.solveDirect(b, **self.options) elif self.flag is None and not self.doDirect: @@ -38,6 +70,14 @@ class Solver(object): self.dsolve = None def solveDirect(self, b, factorize=False, backend='scipy'): + """ + Use solve instead of this interface. + + :param bool factorize: if you want to factorize and store factors + :param str backend: which backend to use. Default is scipy + :rtype: numpy.ndarray + :return: x + """ assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch' if factorize and self.dsolve is None: @@ -64,8 +104,16 @@ class Solver(object): def solveIter(self, b, M=None, iterSolver='CG'): pass - def solveBackward(self, b): - "Perform a backwards solve with upper triangular A in CSR format." + def solveBackward(self, b, backend='python'): + """ + Use solve instead of this interface. + + Perform a backwards solve with upper triangular A in CSR format (best, if not, it will be converted). + + :param str backend: which backend to use. Default is python. + :rtype: numpy.ndarray + :return: x + """ if type(self.A) is not sparse.csr.csr_matrix: from scipy.sparse import csr_matrix self.A = csr_matrix(self.A) @@ -80,8 +128,16 @@ class Solver(object): x[i] = (b[i] - np.dot(ith_row[1:], x_vals[1:])) / ith_row[0] return x - def solveForward(self, b): - "Perform a forward solve with lower triangular A in CSR format." + def solveForward(self, backend='python'): + """ + Use solve instead of this interface. + + Perform a forward solve with lower triangular A in CSR format (best, if not, it will be converted). + + :param str backend: which backend to use. Default is python. + :rtype: numpy.ndarray + :return: x + """ if type(self.A) is not sparse.csr.csr_matrix: from scipy.sparse import csr_matrix self.A = csr_matrix(self.A) @@ -96,7 +152,16 @@ class Solver(object): x[i] = (b[i] - np.dot(ith_row[:-1], x_vals[:-1])) / ith_row[-1] return x - def solveDiagonal(self, b): + def solveDiagonal(self, backend='python'): + """ + Use solve instead of this interface. + + Perform a diagonal solve with diagonal matrix A. + + :param str backend: which backend to use. Default is python. + :rtype: numpy.ndarray + :return: x + """ diagA = self.A.diagonal() if len(b.shape) == 1 or b.shape[1] == 1: # Just one RHS diff --git a/docs/api_Solver.rst b/docs/api_Solver.rst new file mode 100644 index 00000000..74a81707 --- /dev/null +++ b/docs/api_Solver.rst @@ -0,0 +1,9 @@ +.. _api_Solver: + +Solver +****** + +.. automodule:: SimPEG.Solver + :members: + :undoc-members: + diff --git a/docs/index.rst b/docs/index.rst index 692363a2..fa127864 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -62,6 +62,7 @@ Utility Codes .. toctree:: :maxdepth: 2 + api_Solver api_Utils From 11128b2e9ba40ef89ff29487211b2da8ab843ba6 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 30 Oct 2013 23:13:49 -0600 Subject: [PATCH 21/34] Bug Fixes. Test_DCproblem is still not working. --- SimPEG/Solver.py | 4 ++-- SimPEG/forward/Problem.py | 1 + SimPEG/forward/__init__.py | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index 127a8c9f..d9ac4a1d 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -128,7 +128,7 @@ class Solver(object): x[i] = (b[i] - np.dot(ith_row[1:], x_vals[1:])) / ith_row[0] return x - def solveForward(self, backend='python'): + def solveForward(self, b, backend='python'): """ Use solve instead of this interface. @@ -152,7 +152,7 @@ class Solver(object): x[i] = (b[i] - np.dot(ith_row[:-1], x_vals[:-1])) / ith_row[-1] return x - def solveDiagonal(self, backend='python'): + def solveDiagonal(self, b, backend='python'): """ Use solve instead of this interface. diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 679e8686..395156a2 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -1,5 +1,6 @@ import numpy as np from SimPEG.utils import mkvc, sdiag +import scipy.sparse as sp norm = np.linalg.norm diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py index a7c58967..33c9a6b1 100644 --- a/SimPEG/forward/__init__.py +++ b/SimPEG/forward/__init__.py @@ -1,4 +1,4 @@ from Problem import * -from DCProblem import DCProblem +import DCProblem from LinearProblem import LinearProblem import ModelTransforms From 831eca50221accde729c3cb22c6d957f0b38804b Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 30 Oct 2013 23:28:14 -0600 Subject: [PATCH 22/34] Some documentation. --- SimPEG/tests/TestUtils.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index 83738a75..ff7ee668 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -192,6 +192,16 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): :rtype: bool :return: did you pass the test?! + + .. plot:: + :include-source: + + from SimPEG.tests import checkDerivative + from SimPEG.utils import sdiag + import numpy as np + def simplePass(x): + return np.sin(x), sdiag(np.cos(x)) + checkDerivative(simplePass, np.random.randn(5)) """ print "%s checkDerivative %s" % ('='*20, '='*20) From f44a85766978b4b0101f7424dc6d30be6c74f622 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 1 Nov 2013 15:00:17 -0700 Subject: [PATCH 23/34] Moved solver to utilities. But it is still available as SimPEG.Solver (as defined in the main __init__.py file) --- SimPEG/__init__.py | 4 ++-- SimPEG/{ => utils}/Solver.py | 0 SimPEG/utils/__init__.py | 2 ++ docs/api_Solver.rst | 2 +- 4 files changed, 5 insertions(+), 3 deletions(-) rename SimPEG/{ => utils}/Solver.py (100%) diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index a17eaa74..7f059a74 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1,6 +1,6 @@ -import mesh import utils +from utils import Solver +import mesh import inverse import forward import regularization -from Solver import Solver diff --git a/SimPEG/Solver.py b/SimPEG/utils/Solver.py similarity index 100% rename from SimPEG/Solver.py rename to SimPEG/utils/Solver.py diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 7c976ce3..d5918b0d 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -2,6 +2,8 @@ import matutils import sputils import lomutils import ModelBuilder +import Solver +from Solver import Solver from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind from sputils import spzeros, kron3, speye, sdiag from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird diff --git a/docs/api_Solver.rst b/docs/api_Solver.rst index 74a81707..db75860b 100644 --- a/docs/api_Solver.rst +++ b/docs/api_Solver.rst @@ -3,7 +3,7 @@ Solver ****** -.. automodule:: SimPEG.Solver +.. automodule:: SimPEG.utils.Solver :members: :undoc-members: From 581f1b15f137e414fc3ed628ffaf56c2f082ec23 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 1 Nov 2013 15:02:37 -0700 Subject: [PATCH 24/34] renamed misfit to dataResidual in Problem class. fixed betaSchedule bug. Printing in the inversion. --- SimPEG/forward/Problem.py | 2 +- SimPEG/inverse/BetaSchedule.py | 2 +- SimPEG/inverse/Inversion.py | 22 +++++++++++++++++----- 3 files changed, 19 insertions(+), 7 deletions(-) diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 395156a2..32ab663f 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -94,7 +94,7 @@ class Problem(object): u = self.field(m) return self.P*u - def misfit(self, m, u=None): + def dataResidual(self, m, u=None): """ :param numpy.array m: geophysical model :param numpy.array u: fields diff --git a/SimPEG/inverse/BetaSchedule.py b/SimPEG/inverse/BetaSchedule.py index 2c633a6a..fe197340 100644 --- a/SimPEG/inverse/BetaSchedule.py +++ b/SimPEG/inverse/BetaSchedule.py @@ -9,4 +9,4 @@ class Cooling(object): def getBeta(self): if self._beta is None: return beta0 - return self._beta * beta_coolingFactor + return self._beta / beta_coolingFactor diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 02a81daa..b58a8318 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -6,6 +6,7 @@ class Inversion(object): """docstring for Inversion""" maxIter = 10 + name = 'SimPEG Inversion' def __init__(self, prob, reg, opt, **kwargs): self.prob = prob @@ -22,6 +23,14 @@ class Inversion(object): else: raise Exception('%s attr is not recognized' % attr) + def printInit(self): + print "%s %s %s" % ('='*22, self.name, '='*22) + print " # beta phi_d phi_m f norm(dJ) #LS" + print "%s" % '-'*62 + + def printIter(self): + print "%3d %1.2e %1.2e %1.2e %1.2e %1.2e %3d" % (self.opt._iter, self._beta, self._phi_d_last, self._phi_m_last, self.opt.f, np.linalg.norm(self.opt.g), self.opt._iterLS) + @property def Wd(self): """ @@ -65,7 +74,7 @@ class Inversion(object): def getBeta(self): if self._beta is None: return self.beta0 - return self._beta * self.beta_coolingFactor + return self._beta / self.beta_coolingFactor def stoppingCriteria(self): self._STOP = np.zeros(2,dtype=bool) @@ -121,9 +130,10 @@ class Inversion(object): 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(m, u=u) + # TODO: ensure that this is a data is vector and Wd is a matrix. + R = self.Wd*self.prob.dataResidual(m, u=u) R = mkvc(R) - return 0.5*R.dot(R) + return 0.5*np.vdot(R, R) def dataObjDeriv(self, m, u=None): """ @@ -159,7 +169,7 @@ class Inversion(object): if u is None: u = self.prob.field(m) - R = self.Wd*self.prob.misfit(m, u=u) + R = self.Wd*self.prob.dataResidual(m, u=u) dmisfit = self.prob.Jt(m, self.Wd * R, u=u) @@ -201,8 +211,10 @@ class Inversion(object): if u is None: u = self.prob.field(m) - R = self.Wd*self.prob.misfit(m, u=u) + R = self.Wd*self.prob.dataResidual(m, u=u) + # TODO: abstract to different norms a little cleaner. + # \/ it goes here. in l2 it is the identity. dmisfit = self.prob.Jt(m, self.Wd * self.Wd * self.prob.J(m, v, u=u), u=u) return dmisfit From 510ca4a54e8aa41fccb89d66a5d3285273fca312 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 1 Nov 2013 15:04:16 -0700 Subject: [PATCH 25/34] Issue #11 PubSub based communications in Minimize. --- SimPEG/inverse/Optimize.py | 50 ++++++++++++++++++++++++++++++-------- docs/requirements.txt | 1 + 2 files changed, 41 insertions(+), 10 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 9e3b5463..d75b64d7 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -3,6 +3,7 @@ import matplotlib.pyplot as plt from SimPEG.utils import mkvc, sdiag norm = np.linalg.norm import scipy.sparse as sp +from pubsub import pub class Minimize(object): @@ -24,9 +25,8 @@ class Minimize(object): tolG = 1e-1 eps = 1e-5 - printIter = [] # push to here if you want to print these on iter - def __init__(self, **kwargs): + self._id = int(np.random.rand()*1e6) # create a unique identifier to this program to be used in pubsub self.setKwargs(**kwargs) def setKwargs(self, **kwargs): @@ -51,13 +51,19 @@ class Minimize(object): while True: self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) + pub.sendMessage('Minimize.evalFunction', minimize=self, f=self.f, g=self.g, H=self.H) self.printIter() if self.stoppingCriteria(): break p = self.findSearchDirection() - xt, passLS = self.linesearch(p) + #TODO: Scale search direction + pub.sendMessage('Minimize.searchDirection', minimize=self, p=p) + xt, passLS = self.linesearch(p) ## TODO: should be called modifyStep to be inclusive of trust region stuff etc. + pub.sendMessage('Minimize.linesearch', minimize=self, xt=xt) if not passLS: xt = self.linesearchBreak(p) + return self.xc self.doEndIteration(xt) + pub.sendMessage('Minimize.endIteration', minimize=self, xt=xt) self.printDone() @@ -65,7 +71,9 @@ class Minimize(object): @property def parent(self): - """This is the parent of the optimization routine.""" + """ + This is the parent of the optimization routine. + """ return getattr(self, '_parent', None) @parent.setter def parent(self, value): @@ -85,6 +93,10 @@ class Minimize(object): printIter is called at the beginning of the optimization routine. """ + pub.sendMessage('Minimize.printInit', minimize=self) + if self.parent is not None and hasattr(self.parent, 'printInit'): + self.parent.printInit() + return print "%s %s %s" % ('='*22, self.name, '='*22) print "iter\tJc\t\tnorm(dJ)\tLS" print "%s" % '-'*57 @@ -94,12 +106,22 @@ class Minimize(object): printIter is called directly after function evaluations. """ + pub.sendMessage('Minimize.printIter', minimize=self) + if self.parent is not None and hasattr(self.parent, 'printIter'): + self.parent.printIter() + return print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS) def printDone(self): + pub.sendMessage('Minimize.printDone', minimize=self) + if self.parent is not None and hasattr(self.parent, 'printDone'): + self.parent.printDone() + return print "%s STOP! %s" % ('-'*25,'-'*25) - print "%d : |fc-fOld| = %1.4e <= tolF*(1+|fStop|) = %1.4e" % (self._STOP[0], abs(self.f-self.fOld), self.tolF*(1+abs(self.fStop))) - print "%d : |xc-xOld| = %1.4e <= tolX*(1+|x0|) = %1.4e" % (self._STOP[1], norm(self.xc-self.xOld), self.tolX*(1+norm(self.x0))) + # TODO: put controls on gradient value, min model update, and function value + if self._iter > 0: + print "%d : |fc-fOld| = %1.4e <= tolF*(1+|fStop|) = %1.4e" % (self._STOP[0], abs(self.f-self.fOld), self.tolF*(1+abs(self.fStop))) + print "%d : |xc-xOld| = %1.4e <= tolX*(1+|x0|) = %1.4e" % (self._STOP[1], norm(self.xc-self.xOld), self.tolX*(1+norm(self.x0))) print "%d : |g| = %1.4e <= tolG*(1+|fStop|) = %1.4e" % (self._STOP[2], norm(self.g), self.tolG*(1+abs(self.fStop))) print "%d : |g| = %1.4e <= 1e3*eps = %1.4e" % (self._STOP[3], norm(self.g), 1e3*self.eps) print "%d : iter = %3d\t <= maxIter\t = %3d" % (self._STOP[4], self._iter, self.maxIter) @@ -140,7 +162,7 @@ class Minimize(object): return xt, iterLS < self.maxIterLS def linesearchBreak(self, p): - raise Exception('The linesearch got broken. Boo.') + print 'The linesearch got broken. Boo.' def doEndIteration(self, xt): # store old values @@ -159,7 +181,7 @@ class InexactGaussNewton(Minimize): name = 'InexactGaussNewton' def findSearchDirection(self): # TODO: use BFGS as a preconditioner or gauss sidel of the WtW or solve WtW directly - p, info = sp.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=5) + p, info = sp.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=10) return p @@ -170,11 +192,19 @@ class SteepestDescent(Minimize): if __name__ == '__main__': from SimPEG.tests import Rosenbrock, checkDerivative + import matplotlib.pyplot as plt x0 = np.array([2.6, 3.7]) checkDerivative(Rosenbrock, x0, plotIt=False) - xOpt = GaussNewton(maxIter=20).minimize(Rosenbrock,x0) + + def listener1(minimize,p): + plt.plot(p) + plt.show() + print p + pub.subscribe(listener1, 'Minimize.searchDirection') + + xOpt = GaussNewton(maxIter=20,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock,x0) print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) - xOpt = SteepestDescent(maxIter=20, maxIterLS=15).minimize(Rosenbrock, x0) + xOpt = SteepestDescent(maxIter=30, maxIterLS=15,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock, x0) print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) def simplePass(x): diff --git a/docs/requirements.txt b/docs/requirements.txt index 24ce15ab..d7809121 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1 +1,2 @@ numpy +pubsub From 187ba97112f97c9267373536ef011f1ef0e345de Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 1 Nov 2013 15:19:07 -0700 Subject: [PATCH 26/34] Bug fix in regularization. --- SimPEG/mesh/DiffOperators.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 380e0d19..598b392a 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -192,7 +192,7 @@ class DiffOperators(object): BC = ['neumann', 'neumann'] n = self.n if(self.dim == 2): - G2 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC)) + G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0])) elif(self.dim == 3): G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0])) # Compute areas of cell faces & volumes From a1a6e679257012e5360a1e684077fb4b6f62ed83 Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Sun, 3 Nov 2013 10:36:22 -0800 Subject: [PATCH 27/34] Added max step size option to Optimize code. --- SimPEG/inverse/Optimize.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index d75b64d7..01a6f453 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -18,6 +18,7 @@ class Minimize(object): maxIter = 20 maxIterLS = 10 + maxStep = np.inf LSreduction = 1e-4 LSshorten = 0.5 tolF = 1e-1 @@ -55,7 +56,8 @@ class Minimize(object): self.printIter() if self.stoppingCriteria(): break p = self.findSearchDirection() - #TODO: Scale search direction + if self.maxStep < np.abs(p.max()): + p = self.maxStep*p/np.abs(p.max()) pub.sendMessage('Minimize.searchDirection', minimize=self, p=p) xt, passLS = self.linesearch(p) ## TODO: should be called modifyStep to be inclusive of trust region stuff etc. pub.sendMessage('Minimize.linesearch', minimize=self, xt=xt) From 224b0cea9e2fbf65aec8b06f32398d9304f6a119 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 09:43:49 -0800 Subject: [PATCH 28/34] Issue #17 Scaled gradient. moved this to an overridable function. Changed linesearch to modifySearchDirection --- SimPEG/inverse/Optimize.py | 49 ++++++++++++++++++++++++++------------ 1 file changed, 34 insertions(+), 15 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 01a6f453..d4c46c69 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -18,7 +18,7 @@ class Minimize(object): maxIter = 20 maxIterLS = 10 - maxStep = np.inf + maxStep = np.inf LSreduction = 1e-4 LSshorten = 0.5 tolF = 1e-1 @@ -56,14 +56,14 @@ class Minimize(object): self.printIter() if self.stoppingCriteria(): break p = self.findSearchDirection() - if self.maxStep < np.abs(p.max()): - p = self.maxStep*p/np.abs(p.max()) pub.sendMessage('Minimize.searchDirection', minimize=self, p=p) - xt, passLS = self.linesearch(p) ## TODO: should be called modifyStep to be inclusive of trust region stuff etc. - pub.sendMessage('Minimize.linesearch', minimize=self, xt=xt) + p = self.scaleSearchDirection(p) + pub.sendMessage('Minimize.scaleSearchDirection', minimize=self, p=p) + xt, passLS = self.modifySearchDirection(p) + pub.sendMessage('Minimize.modifySearchDirection', minimize=self, xt=xt) if not passLS: - xt = self.linesearchBreak(p) - return self.xc + xt, caught = self.modifySearchDirectionBreak(p) + if not caught: return self.xc self.doEndIteration(xt) pub.sendMessage('Minimize.endIteration', minimize=self, xt=xt) @@ -129,9 +129,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 findSearchDirection(self): - return -self.g - def stoppingCriteria(self): if self._iter == 0: self.fStop = self.f # Save this for stopping criteria @@ -147,7 +144,15 @@ class Minimize(object): def projection(self, p): return p - def linesearch(self, p): + def findSearchDirection(self): + return -self.g + + def scaleSearchDirection(self, p): + if self.maxStep < np.abs(p.max()): + p = self.maxStep*p/np.abs(p.max()) + return p + + def modifySearchDirection(self, p): # Armijo linesearch descent = np.inner(self.g, p) t = 1 @@ -163,8 +168,24 @@ class Minimize(object): self._iterLS = iterLS return xt, iterLS < self.maxIterLS - def linesearchBreak(self, p): + def modifySearchDirectionBreak(self, p): + """ + Code is called if modifySearchDirection fails + to find a descent direction. + + The search direction is passed as input and + this function must pass back both a new searchDirection, + and if the searchDirection break has been caught. + + By default, no additional work is done, and the + function returns a False indicating the break was not caught. + + :param numpy.ndarray p: searchDirection + :rtype: numpy.ndarray,bool + :return: (xt, breakCaught) + """ print 'The linesearch got broken. Boo.' + return p, False def doEndIteration(self, xt): # store old values @@ -199,9 +220,7 @@ if __name__ == '__main__': checkDerivative(Rosenbrock, x0, plotIt=False) def listener1(minimize,p): - plt.plot(p) - plt.show() - print p + print 'hi: ', p pub.subscribe(listener1, 'Minimize.searchDirection') xOpt = GaussNewton(maxIter=20,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock,x0) From 56bacc4a27b79841a6ed0a70992261f508f5ba7f Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 09:53:38 -0800 Subject: [PATCH 29/34] Issue #15 Approximate Hessians. By default these are not implemented, and are the same as J and Jt --- SimPEG/forward/Problem.py | 30 +++++++++++++++++++++++++++++- SimPEG/inverse/Inversion.py | 2 +- 2 files changed, 30 insertions(+), 2 deletions(-) diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 32ab663f..2e6831f7 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -150,10 +150,38 @@ class Problem(object): :rtype: numpy.array :return: JTv - Transpose of J + Effect of transpose of J on a vector v. """ pass + + def J_approx(self, m, v, u=None): + """ + + :param numpy.array m: model + :param numpy.array v: vector to multiply + :param numpy.array u: fields + :rtype: numpy.array + :return: Jv + + Approximate effect of J on a vector v + + """ + return self.J(m, v, u) + + def Jt_approx(self, m, v, u=None): + """ + :param numpy.array m: model + :param numpy.array v: vector to multiply + :param numpy.array u: fields + :rtype: numpy.array + :return: JTv + + Approximate transpose of J*v + + """ + return self.Jt(m, v, u) + def field(self, m): """ The field given the model. diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index b58a8318..827e12f8 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -215,7 +215,7 @@ class Inversion(object): # TODO: abstract to different norms a little cleaner. # \/ it goes here. in l2 it is the identity. - dmisfit = self.prob.Jt(m, self.Wd * self.Wd * self.prob.J(m, v, u=u), u=u) + dmisfit = self.prob.Jt_approx(m, self.Wd * self.Wd * self.prob.J_approx(m, v, u=u), u=u) return dmisfit From 61189d2e5b0b39ad61179dfbb52f95b71554b206 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 10:11:14 -0800 Subject: [PATCH 30/34] Documentation updates. --- docs/api_LOMView.rst | 8 -------- docs/api_LogicallyOrthogonalMesh.rst | 8 ++++++++ docs/api_TensorMesh.rst | 7 +++++++ docs/api_TensorView.rst | 8 -------- docs/index.rst | 7 ------- 5 files changed, 15 insertions(+), 23 deletions(-) delete mode 100644 docs/api_LOMView.rst delete mode 100644 docs/api_TensorView.rst diff --git a/docs/api_LOMView.rst b/docs/api_LOMView.rst deleted file mode 100644 index cfbfa6fa..00000000 --- a/docs/api_LOMView.rst +++ /dev/null @@ -1,8 +0,0 @@ -.. _api_LOMView: - -LOM View -******** - -.. automodule:: SimPEG.mesh.LomView - :members: - :undoc-members: diff --git a/docs/api_LogicallyOrthogonalMesh.rst b/docs/api_LogicallyOrthogonalMesh.rst index af9e73d0..9d7d516b 100644 --- a/docs/api_LogicallyOrthogonalMesh.rst +++ b/docs/api_LogicallyOrthogonalMesh.rst @@ -6,3 +6,11 @@ Logically Orthogonal Mesh .. automodule:: SimPEG.mesh.LogicallyOrthogonalMesh :members: :undoc-members: + + +LOM View +******** + +.. automodule:: SimPEG.mesh.LomView + :members: + :undoc-members: diff --git a/docs/api_TensorMesh.rst b/docs/api_TensorMesh.rst index 36fb925d..9de206f1 100644 --- a/docs/api_TensorMesh.rst +++ b/docs/api_TensorMesh.rst @@ -6,3 +6,10 @@ Tensor Mesh .. automodule:: SimPEG.mesh.TensorMesh :members: :undoc-members: + +Tensor View +*********** + +.. automodule:: SimPEG.mesh.TensorView + :members: + :undoc-members: diff --git a/docs/api_TensorView.rst b/docs/api_TensorView.rst deleted file mode 100644 index c6b5a9b8..00000000 --- a/docs/api_TensorView.rst +++ /dev/null @@ -1,8 +0,0 @@ -.. _api_TensorView: - -Tensor View -*********** - -.. automodule:: SimPEG.mesh.TensorView - :members: - :undoc-members: diff --git a/docs/index.rst b/docs/index.rst index fa127864..981e5ddc 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -1,8 +1,3 @@ -.. SimPEG documentation master file, created by - sphinx-quickstart on Fri Aug 30 18:42:44 2013. - You can adapt this file completely to your liking, but it should at least - contain the root `toctree` directive. - SimPEG ====== @@ -24,10 +19,8 @@ Meshing & Operators api_BaseMesh api_TensorMesh - api_TensorView api_LogicallyOrthogonalMesh api_Cyl1DMesh - api_LOMView api_DiffOperators api_InnerProducts From faa2be56a58f33396e2078afd99444d262dfd715 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 10:17:16 -0800 Subject: [PATCH 31/34] update requirements.txt for docs. --- docs/requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/requirements.txt b/docs/requirements.txt index d7809121..dcecbe99 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,2 +1,2 @@ numpy -pubsub +pypubsub From ac099f08cd4f7171546c0cc8524834f04eba710b Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 10:26:56 -0800 Subject: [PATCH 32/34] There are some compatibility issues with pypubsub installing on the documentation servers. and confusion when installing this on other peoples comps. I have put a warning in. Code should work regardless of if you have pypubsub installed or not. --- SimPEG/inverse/Optimize.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index d4c46c69..3e50e3df 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -3,7 +3,14 @@ import matplotlib.pyplot as plt from SimPEG.utils import mkvc, sdiag norm = np.linalg.norm import scipy.sparse as sp -from pubsub import pub + +try: + from pubsub import pub + doPub = True +except Exception, e: + print 'Warning: you may not have the required pubsub installed, use pypubsub. You will not be able to listen to events.' + doPub = False + class Minimize(object): @@ -52,20 +59,20 @@ class Minimize(object): while True: self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) - pub.sendMessage('Minimize.evalFunction', minimize=self, f=self.f, g=self.g, H=self.H) + if doPub: pub.sendMessage('Minimize.evalFunction', minimize=self, f=self.f, g=self.g, H=self.H) self.printIter() if self.stoppingCriteria(): break p = self.findSearchDirection() - pub.sendMessage('Minimize.searchDirection', minimize=self, p=p) + if doPub: pub.sendMessage('Minimize.searchDirection', minimize=self, p=p) p = self.scaleSearchDirection(p) - pub.sendMessage('Minimize.scaleSearchDirection', minimize=self, p=p) + if doPub: pub.sendMessage('Minimize.scaleSearchDirection', minimize=self, p=p) xt, passLS = self.modifySearchDirection(p) - pub.sendMessage('Minimize.modifySearchDirection', minimize=self, xt=xt) + if doPub: pub.sendMessage('Minimize.modifySearchDirection', minimize=self, xt=xt) if not passLS: xt, caught = self.modifySearchDirectionBreak(p) if not caught: return self.xc self.doEndIteration(xt) - pub.sendMessage('Minimize.endIteration', minimize=self, xt=xt) + if doPub: pub.sendMessage('Minimize.endIteration', minimize=self, xt=xt) self.printDone() @@ -95,7 +102,7 @@ class Minimize(object): printIter is called at the beginning of the optimization routine. """ - pub.sendMessage('Minimize.printInit', minimize=self) + if doPub: pub.sendMessage('Minimize.printInit', minimize=self) if self.parent is not None and hasattr(self.parent, 'printInit'): self.parent.printInit() return @@ -108,14 +115,14 @@ class Minimize(object): printIter is called directly after function evaluations. """ - pub.sendMessage('Minimize.printIter', minimize=self) + if doPub: pub.sendMessage('Minimize.printIter', minimize=self) if self.parent is not None and hasattr(self.parent, 'printIter'): self.parent.printIter() return print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS) def printDone(self): - pub.sendMessage('Minimize.printDone', minimize=self) + if doPub: pub.sendMessage('Minimize.printDone', minimize=self) if self.parent is not None and hasattr(self.parent, 'printDone'): self.parent.printDone() return @@ -221,7 +228,7 @@ if __name__ == '__main__': def listener1(minimize,p): print 'hi: ', p - pub.subscribe(listener1, 'Minimize.searchDirection') + if doPub: pub.subscribe(listener1, 'Minimize.searchDirection') xOpt = GaussNewton(maxIter=20,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock,x0) print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) From df51919d8135ab0d8a94da0ecca3e615246b2f62 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 11:49:55 -0800 Subject: [PATCH 33/34] Documentation Updates. Issue #11 --- SimPEG/inverse/Optimize.py | 166 +++++++++++++++++++++++++++++++++-- SimPEG/mesh/InnerProducts.py | 8 +- 2 files changed, 163 insertions(+), 11 deletions(-) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 3e50e3df..d4829b7b 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -3,6 +3,7 @@ import matplotlib.pyplot as plt from SimPEG.utils import mkvc, sdiag norm = np.linalg.norm import scipy.sparse as sp +from SimPEG import Solver try: from pubsub import pub @@ -16,7 +17,7 @@ except Exception, e: class Minimize(object): """ - Minimize is a general class for derivative based optimization. + Minimize is a general class for derivative based optimization. """ @@ -38,6 +39,7 @@ class Minimize(object): self.setKwargs(**kwargs) def setKwargs(self, **kwargs): + """Sets key word arguments (kwargs) that are present in the object.""" # Set the variables, throw an error if they don't exist. for attr in kwargs: if hasattr(self, attr): @@ -47,11 +49,58 @@ class Minimize(object): def minimize(self, evalFunction, x0): """ + Minimizes the function (evalFunction) starting at the location x0. + + :param def evalFunction: function handle that evaluates: f, g, H = F(x) + :param numpy.ndarray x0: starting location + :rtype: numpy.ndarray + :return: x, the last iterate of the optimization algorithm evalFunction is a function handle:: - evalFunction(x, return_g=True, return_H=True ) + (f[, g][, H]) = evalFunction(x, return_g=False, return_H=False ) + + Events are fired with the following inputs via pypubsub:: + + Minimize.printInit (minimize) + Minimize.evalFunction (minimize, f, g, H) + Minimize.printIter (minimize) + Minimize.searchDirection (minimize, p) + Minimize.scaleSearchDirection (minimize, p) + Minimize.modifySearchDirection (minimize, xt, passLS) + Minimize.endIteration (minimize, xt) + Minimize.printDone (minimize) + + To hook into one of these events (must have pypubsub installed):: + + from pubsub import pub + def listener(minimize,p): + print 'The search direction is: ', p + pub.subscribe(listener, 'Minimize.searchDirection') + + You can use pubsub communication to debug your code, it is not used internally. + + + The algorithm for general minimization is as follows:: + + startup(x0) + printInit() + + while True: + f, g, H = evalFunction(xc) + printIter() + if stoppingCriteria(): break + p = findSearchDirection() + p = scaleSearchDirection(p) + xt, passLS = modifySearchDirection(p) + if not passLS: + xt, caught = modifySearchDirectionBreak(p) + if not caught: return xc + doEndIteration(xt) + + printDone() + return xc """ self.evalFunction = evalFunction self.startup(x0) @@ -67,7 +116,7 @@ class Minimize(object): p = self.scaleSearchDirection(p) if doPub: pub.sendMessage('Minimize.scaleSearchDirection', minimize=self, p=p) xt, passLS = self.modifySearchDirection(p) - if doPub: pub.sendMessage('Minimize.modifySearchDirection', minimize=self, xt=xt) + if doPub: pub.sendMessage('Minimize.modifySearchDirection', minimize=self, xt=xt, passLS=passLS) if not passLS: xt, caught = self.modifySearchDirectionBreak(p) if not caught: return self.xc @@ -89,6 +138,19 @@ class Minimize(object): self._parent = value def startup(self, x0): + """ + **startup** is called at the start of any new minimize call. + + This will set:: + + x0 = x0 + xc = x0 + _iter = _iterLS = 0 + + :param numpy.ndarray x0: initial x + :rtype: None + :return: None + """ self._iter = 0 self._iterLS = 0 self._STOP = np.zeros((5,1),dtype=bool) @@ -99,7 +161,10 @@ class Minimize(object): def printInit(self): """ - printIter is called at the beginning of the optimization routine. + printInit is called at the beginning of the optimization routine. + + If there is a parent object, printInit will check for a + parent.printInit function and call that. """ if doPub: pub.sendMessage('Minimize.printInit', minimize=self) @@ -114,6 +179,9 @@ class Minimize(object): """ printIter is called directly after function evaluations. + If there is a parent object, printIter will check for a + parent.printIter function and call that. + """ if doPub: pub.sendMessage('Minimize.printIter', minimize=self) if self.parent is not None and hasattr(self.parent, 'printIter'): @@ -122,6 +190,13 @@ class Minimize(object): print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS) def printDone(self): + """ + printDone is called at the end of the optimization routine. + + If there is a parent object, printDone will check for a + parent.printDone function and call that. + + """ if doPub: pub.sendMessage('Minimize.printDone', minimize=self) if self.parent is not None and hasattr(self.parent, 'printDone'): self.parent.printDone() @@ -149,17 +224,79 @@ class Minimize(object): return all(self._STOP[0:3]) | any(self._STOP[3:]) def projection(self, p): + """ + projects the search direction. + + by default, no projection is applied. + + :param numpy.ndarray p: searchDirection + :rtype: numpy.ndarray + :return: p, projected search direction + """ return p def findSearchDirection(self): + """ + **findSearchDirection** should return an approximation of: + + .. math:: + + H p = - g + + Where you are solving for the search direction, p + + The default is: + + .. math:: + + H = I + + p = - g + + And corresponds to SteepestDescent. + + The latest function evaluations are present in:: + + self.f, self.g, self.H + + :rtype: numpy.ndarray + :return: p, Search Direction + """ return -self.g def scaleSearchDirection(self, p): + """ + **scaleSearchDirection** should scale the search direction if appropriate. + + Set the parameter **maxStep** in the minimize object, to scale back the gradient to a maximum size. + + :param numpy.ndarray p: searchDirection + :rtype: numpy.ndarray + :return: p, Scaled Search Direction + """ + if self.maxStep < np.abs(p.max()): p = self.maxStep*p/np.abs(p.max()) return p def modifySearchDirection(self, p): + """ + **modifySearchDirection** changes the search direction based on some sort of linesearch or trust-region criteria. + + By default, an Armijo backtracking linesearch is preformed with the following parameters: + + * maxIterLS, the maximum number of linesearch iterations + * LSreduction, the expected reduction expected, default: 1e-4 + * LSshorten, how much the step is reduced, default: 0.5 + + If the linesearch is completed, and a descent direction is found, passLS is returned as True. + + Else, a modifySearchDirectionBreak call is preformed. + + :param numpy.ndarray p: searchDirection + :rtype: numpy.ndarray,bool + :return: (xt, passLS) + """ # Armijo linesearch descent = np.inner(self.g, p) t = 1 @@ -185,7 +322,7 @@ class Minimize(object): and if the searchDirection break has been caught. By default, no additional work is done, and the - function returns a False indicating the break was not caught. + evalFunction returns a False indicating the break was not caught. :param numpy.ndarray p: searchDirection :rtype: numpy.ndarray,bool @@ -195,6 +332,17 @@ class Minimize(object): return p, False def doEndIteration(self, xt): + """ + **doEndIteration** is called at the end of each minimize iteration. + + By default, function values and x locations are shuffled to store 1 past iteration in memory. + + self.xc must be updated in this code. + + :param numpy.ndarray xt: tested new iterate that ensures a descent direction. + :rtype: None + :return: None + """ # store old values self.fOld = self.f self.xOld, self.xc = self.xc, xt @@ -204,14 +352,18 @@ class Minimize(object): class GaussNewton(Minimize): name = 'GaussNewton' def findSearchDirection(self): - return np.linalg.solve(self.H,-self.g) + return Solver(self.H).solve(-self.g) class InexactGaussNewton(Minimize): name = 'InexactGaussNewton' + + maxIterCG = 10 + tolCG = 1e-5 + def findSearchDirection(self): # TODO: use BFGS as a preconditioner or gauss sidel of the WtW or solve WtW directly - p, info = sp.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=10) + p, info = sp.linalg.cg(self.H, -self.g, tol=self.tolCG, maxiter=self.maxIterCG) return p diff --git a/SimPEG/mesh/InnerProducts.py b/SimPEG/mesh/InnerProducts.py index 9fe84ac3..f0ac4ab0 100644 --- a/SimPEG/mesh/InnerProducts.py +++ b/SimPEG/mesh/InnerProducts.py @@ -81,9 +81,9 @@ class InnerProducts(object): def getFaceInnerProduct(self, mu=None, returnP=False): """Wrapper function, - :py:func:`SimPEG.InnerProducts.getEdgeInnerProduct` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct` - :py:func:`SimPEG.InnerProducts.getEdgeInnerProduct2D` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct2D` """ if self.dim == 2: return getFaceInnerProduct2D(self, mu, returnP) @@ -93,9 +93,9 @@ class InnerProducts(object): def getEdgeInnerProduct(self, sigma=None, returnP=False): """Wrapper function, - :py:func:`SimPEG.InnerProducts.getFaceInnerProduct` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct` - :py:func:`SimPEG.InnerProducts.getFaceInnerProduct2D` + :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct2D` """ if self.dim == 2: return getEdgeInnerProduct2D(self, sigma, returnP) From 6f141ecaf299e91d973244f173e97950996cf927 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 15:11:10 -0800 Subject: [PATCH 34/34] Updates to DCProblem and testing. --- SimPEG/forward/DCProblem.py | 76 ++++++++++++++------------ SimPEG/inverse/Inversion.py | 2 +- SimPEG/inverse/Optimize.py | 18 ++---- SimPEG/tests/TestUtils.py | 35 ++++++++++-- SimPEG/tests/test_forward_DCproblem.py | 26 +++++++-- SimPEG/tests/test_forward_problem.py | 10 +++- SimPEG/tests/test_utils.py | 23 ++++++++ 7 files changed, 129 insertions(+), 61 deletions(-) diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index c2a615ed..132966a8 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -7,6 +7,7 @@ import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as linalg + class DCProblem(ModelTransforms.LogModel, Problem): """ **DCProblem** @@ -18,6 +19,11 @@ class DCProblem(ModelTransforms.LogModel, Problem): super(DCProblem, self).__init__(mesh) self.mesh.setCellGradBC('neumann') + def reshapeFields(self, u): + if len(u.shape) == 1: + u = u.reshape([-1, self.RHS.shape[1]], order='F') + return u + def createMatrix(self, m): """ Makes the matrix A(m) for the DC resistivity problem. @@ -38,11 +44,25 @@ class DCProblem(ModelTransforms.LogModel, Problem): A = D*Msig*G return A.tocsc() + def dpred(self, m, u=None): + """ + Predicted data. + + .. math:: + d_\\text{pred} = Pu(m) + """ + if u is None: + u = self.field(m) + + u = self.reshapeFields(u) + + return mkvc(self.P*u) + def field(self, m): A = self.createMatrix(m) solve = Solver(A) phi = solve.solve(self.RHS) - return phi + return mkvc(phi) def J(self, m, v, u=None): """ @@ -69,6 +89,8 @@ class DCProblem(ModelTransforms.LogModel, Problem): if u is None: u = self.field(m) + u = self.reshapeFields(u) + P = self.P D = self.mesh.faceDiv G = self.mesh.cellGrad @@ -83,13 +105,18 @@ class DCProblem(ModelTransforms.LogModel, Problem): dCdm[:, i] = D * ( sdiag( G * ui ) * ( Av_dm * ( mT_dm * v ) ) ) solve = Solver(dCdu) - # solve = linalg.factorized(dCdu) - Jv = - P * solve.solve(dCdm) - return Jv + return mkvc(Jv) def Jt(self, m, v, u=None): """Takes data, turns it into a model..ish""" + + if u is None: + u = self.field(m) + + u = self.reshapeFields(u) + v = self.reshapeFields(v) + P = self.P D = self.mesh.faceDiv G = self.mesh.cellGrad @@ -147,7 +174,7 @@ if __name__ == '__main__': import matplotlib.pyplot as plt # Create the mesh - h1 = np.ones(100) + h1 = np.ones(20) h2 = np.ones(100) mesh = TensorMesh([h1,h2]) @@ -156,12 +183,12 @@ if __name__ == '__main__': sig2 = np.log(0.01) # Create a synthetic model from a block in a half-space - p0 = [20, 20] - p1 = [50, 50] + p0 = [5, 10] + p1 = [15, 50] condVals = [sig1, sig2] mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals) plt.colorbar(mesh.plotImage(mSynth)) - # plt.show() + plt.show() # Set up the projection nelec = 50 @@ -184,7 +211,9 @@ if __name__ == '__main__': dobs, Wd = synthetic.createData(mSynth, std=0.05) u = synthetic.field(mSynth) - # mesh.plotImage(u[:,10], showIt=False) + u = synthetic.reshapeFields(u) + mesh.plotImage(u[:,10]) + # plt.show() # Now set up the problem to do some minimization problem = DCProblem(mesh) @@ -194,39 +223,16 @@ if __name__ == '__main__': problem.std = dobs*0 + 0.05 m0 = mesh.gridCC[:,0]*0+sig2 - - - # 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) + opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) reg = Regularization(mesh) - - inv = inverse.Inversion(problem, reg, opt) + inv = inverse.Inversion(problem, reg, opt, beta0=1e4) # Check Derivative derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)] checkDerivative(derChk, mSynth) + print inv.dataObj(m0) print inv.dataObj(mSynth) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 827e12f8..3aa8ce58 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -16,7 +16,7 @@ class Inversion(object): self.setKwargs(**kwargs) def setKwargs(self, **kwargs): - # Set the variables, throw an error if they don't exist. + """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist.""" for attr in kwargs: if hasattr(self, attr): setattr(self, attr, kwargs[attr]) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index d4829b7b..6221246f 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -39,8 +39,7 @@ class Minimize(object): self.setKwargs(**kwargs) def setKwargs(self, **kwargs): - """Sets key word arguments (kwargs) that are present in the object.""" - # Set the variables, throw an error if they don't exist. + """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist.""" for attr in kwargs: if hasattr(self, attr): setattr(self, attr, kwargs[attr]) @@ -161,7 +160,7 @@ class Minimize(object): def printInit(self): """ - printInit is called at the beginning of the optimization routine. + **printInit** is called at the beginning of the optimization routine. If there is a parent object, printInit will check for a parent.printInit function and call that. @@ -177,7 +176,7 @@ class Minimize(object): def printIter(self): """ - printIter is called directly after function evaluations. + **printIter** is called directly after function evaluations. If there is a parent object, printIter will check for a parent.printIter function and call that. @@ -191,7 +190,7 @@ class Minimize(object): def printDone(self): """ - printDone is called at the end of the optimization routine. + **printDone** is called at the end of the optimization routine. If there is a parent object, printDone will check for a parent.printDone function and call that. @@ -386,12 +385,3 @@ if __name__ == '__main__': print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) xOpt = SteepestDescent(maxIter=30, maxIterLS=15,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock, x0) print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) - - def simplePass(x): - return np.sin(x), sdiag(np.cos(x)) - - def simpleFail(x): - return np.sin(x), -sdiag(np.cos(x)) - - checkDerivative(simplePass, np.random.randn(5), plotIt=False) - checkDerivative(simpleFail, np.random.randn(5), plotIt=False) diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index ff7ee668..29f9086f 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -1,12 +1,15 @@ import numpy as np import matplotlib.pyplot as plt from pylab import norm -from SimPEG.utils import mkvc +from SimPEG.utils import mkvc, sdiag from SimPEG import utils from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh import numpy as np import unittest +import inspect +happiness = ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!', 'Testing is important.', 'You are awesome.', 'Go Test Go!', 'Once upon a time, a happy little test passed.', 'And then everyone was happy.'] +sadness = ['No gold star for you.','Try again soon.','Thankfully, persistence is a great substitute for talent.','It might be easier to call this a feature...','Coffee break?', 'Boooooooo :(', 'Testing is important. Do it again.'] class OrderTest(unittest.TestCase): """ @@ -157,9 +160,10 @@ class OrderTest(unittest.TestCase): print '---------------------------------------------' passTest = np.mean(np.array(order)) > self.tolerance*self._expectedOrder if passTest: - print ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!'][np.random.randint(5)] + print happiness[np.random.randint(len(happiness))] else: print 'Failed to pass test on ' + self._meshType + '.' + print sadness[np.random.randint(len(sadness))] print '' self.assertTrue(passTest) @@ -222,7 +226,11 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): for i in range(num): Jt = fctn(x0+t[i]*dx) E0[i] = l2norm(Jt[0]-Jc[0]) # 0th order Taylor - E1[i] = l2norm(Jt[0]-Jc[0]-t[i]*Jc[1].dot(dx)) # 1st order Taylor + if inspect.isfunction(Jc[1]): + E1[i] = l2norm(Jt[0]-Jc[0]-t[i]*Jc[1](dx)) # 1st order Taylor + else: + # We assume it is a numpy.ndarray + E1[i] = l2norm(Jt[0]-Jc[0]-t[i]*Jc[1].dot(dx)) # 1st order Taylor order0 = np.log10(E0[:-1]/E0[1:]) order1 = np.log10(E1[:-1]/E1[1:]) print "%d\t%1.2e\t%1.3e\t\t%1.3e\t\t%1.3f" % (i, t[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1]) @@ -238,9 +246,12 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): passTest = belowTol or correctOrder if passTest: - print "%s PASS! %s\n" % ('='*25, '='*25) + print "%s PASS! %s" % ('='*25, '='*25) + print happiness[np.random.randint(len(happiness))]+'\n' else: print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57) + print sadness[np.random.randint(len(sadness))]+'\n' + if plotIt: plt.figure() @@ -254,3 +265,19 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): plt.show() return passTest + + +if __name__ == '__main__': + + def simplePass(x): + return np.sin(x), sdiag(np.cos(x)) + + def simpleFunction(x): + return np.sin(x), lambda xi: sdiag(np.cos(x))*xi + + def simpleFail(x): + return np.sin(x), -sdiag(np.cos(x)) + + checkDerivative(simplePass, np.random.randn(5), plotIt=False) + checkDerivative(simpleFunction, np.random.randn(5), plotIt=False) + checkDerivative(simpleFail, np.random.randn(5), plotIt=False) diff --git a/SimPEG/tests/test_forward_DCproblem.py b/SimPEG/tests/test_forward_DCproblem.py index 6ffa4650..c6e6f9c2 100644 --- a/SimPEG/tests/test_forward_DCproblem.py +++ b/SimPEG/tests/test_forward_DCproblem.py @@ -3,9 +3,11 @@ import unittest from SimPEG.mesh import TensorMesh from SimPEG.utils import ModelBuilder, sdiag from SimPEG.forward import Problem, SyntheticProblem -from SimPEG.forward.DCProblem import DCProblem, DCutils +from SimPEG.forward.DCProblem import * from TestUtils import checkDerivative from scipy.sparse.linalg import dsolve +from SimPEG.regularization import Regularization +from SimPEG import inverse class DCProblemTests(unittest.TestCase): @@ -34,7 +36,7 @@ class DCProblemTests(unittest.TestCase): elecend = 0.5+spacelec*(nelec-1) elecLocR = np.linspace(elecini, elecend, nelec) rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5 - q, Q, rxmidloc = DCutils.genTxRxmat(nelec, spacelec, surfloc, elecini, mesh) + q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, mesh) P = Q.T # Create some data @@ -52,22 +54,27 @@ class DCProblemTests(unittest.TestCase): problem.RHS = q problem.W = Wd problem.dobs = dobs + problem.std = dobs*0 + 0.05 + opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + reg = Regularization(mesh) + inv = inverse.Inversion(problem, reg, opt, beta0=1e4) + + self.inv = inv + self.reg = reg self.p = problem self.mesh = mesh self.m0 = mSynth self.dobs = dobs - def test_misfit(self): - print 'SimPEG.forward.DCProblem: Testing Misfit' - derChk = lambda m: [self.p.misfit(m), self.p.misfitDeriv(m)] + derChk = lambda m: [self.p.dpred(m), lambda mx: self.p.J(self.m0, mx)] passed = checkDerivative(derChk, self.m0, plotIt=False) self.assertTrue(passed) def test_adjoint(self): # Adjoint Test - u = np.random.rand(self.mesh.nC) + u = np.random.rand(self.mesh.nC*self.p.RHS.shape[1]) v = np.random.rand(self.mesh.nC) w = np.random.rand(self.dobs.shape[0]) wtJv = w.dot(self.p.J(self.m0, v, u=u)) @@ -75,6 +82,13 @@ class DCProblemTests(unittest.TestCase): passed = (wtJv - vtJtw) < 1e-10 self.assertTrue(passed) + def test_dataObj(self): + derChk = lambda m: [self.inv.dataObj(m), self.inv.dataObjDeriv(m)] + checkDerivative(derChk, self.m0, plotIt=False) + + def test_modelObj(self): + derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)] + checkDerivative(derChk, self.m0, plotIt=False) if __name__ == '__main__': diff --git a/SimPEG/tests/test_forward_problem.py b/SimPEG/tests/test_forward_problem.py index d913a697..ba49efd8 100644 --- a/SimPEG/tests/test_forward_problem.py +++ b/SimPEG/tests/test_forward_problem.py @@ -2,6 +2,7 @@ import numpy as np import unittest from SimPEG.mesh import TensorMesh from SimPEG.forward import Problem +from SimPEG.regularization import Regularization from TestUtils import checkDerivative from scipy.sparse.linalg import dsolve @@ -15,7 +16,7 @@ class ProblemTests(unittest.TestCase): c = np.array([1, 4]) self.mesh2 = TensorMesh([a, b], np.array([3, 5])) self.p2 = Problem(self.mesh2) - + self.reg = Regularization(self.mesh2) def test_modelTransform(self): print 'SimPEG.forward.Problem: Testing Model Transform' @@ -23,6 +24,13 @@ class ProblemTests(unittest.TestCase): passed = checkDerivative(lambda m : [self.p2.modelTransform(m), self.p2.modelTransformDeriv(m)], m, plotIt=False) self.assertTrue(passed) + def test_regularization(self): + derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)] + mSynth = np.random.randn(self.mesh2.nC) + checkDerivative(derChk, mSynth, plotIt=False) + + + if __name__ == '__main__': unittest.main() diff --git a/SimPEG/tests/test_utils.py b/SimPEG/tests/test_utils.py index f0b867ec..058dba56 100644 --- a/SimPEG/tests/test_utils.py +++ b/SimPEG/tests/test_utils.py @@ -1,6 +1,28 @@ import numpy as np import unittest from SimPEG.utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal +from SimPEG.tests import checkDerivative + + +class TestCheckDerivative(unittest.TestCase): + + def test_simplePass(self): + def simplePass(x): + return np.sin(x), sdiag(np.cos(x)) + passed = checkDerivative(simplePass, np.random.randn(5), plotIt=False) + self.assertTrue(passed, True) + + def test_simpleFunction(self): + def simpleFunction(x): + return np.sin(x), lambda xi: sdiag(np.cos(x))*xi + passed = checkDerivative(simpleFunction, np.random.randn(5), plotIt=False) + self.assertTrue(passed, True) + + def test_simpleFail(self): + def simpleFail(x): + return np.sin(x), -sdiag(np.cos(x)) + passed = checkDerivative(simpleFail, np.random.randn(5), plotIt=False) + self.assertTrue(not passed, True) class TestSequenceFunctions(unittest.TestCase): @@ -85,5 +107,6 @@ class TestSequenceFunctions(unittest.TestCase): self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-12) + if __name__ == '__main__': unittest.main()