From 41b7dbe82b5415aa26bdda574fa20b2a046953bc Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 22 Oct 2013 14:40:52 -0700 Subject: [PATCH 01/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] =?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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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 06cd641db3482022bafb2ad18171925e24d165ae Mon Sep 17 00:00:00 2001 From: Dave Marchant Date: Mon, 4 Nov 2013 14:07:16 -0800 Subject: [PATCH 34/45] Initial implementation of interpolation matrix generation for TensorMesh (3D ONLY) --- SimPEG/mesh/TensorMesh.py | 110 +++++++++++++++++++++++++++++++++++- SimPEG/utils/__init__.py | 2 + SimPEG/utils/interputils.py | 77 +++++++++++++++++++++++++ 3 files changed, 188 insertions(+), 1 deletion(-) create mode 100644 SimPEG/utils/interputils.py diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index a3329c14..1dbd62a3 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -1,9 +1,10 @@ import numpy as np +import scipy.sparse as sp from BaseMesh import BaseMesh from TensorView import TensorView from DiffOperators import DiffOperators from InnerProducts import InnerProducts -from SimPEG.utils import ndgrid, mkvc +from SimPEG.utils import ndgrid, mkvc, spzeros, interpmat class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): @@ -312,6 +313,113 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): _edge = None edge = property(**edge()) + # --------------- Methods --------------------- + + def isInside(self, pts): + """ + Determines if a set of points are inside a mesh. + + :param numpy.ndarray pts: Location of points to test + :rtype numpy.ndarray + :return inside, numpy array of booleans + """ + + pts = np.atleast_2d(pts) + inside = (pts[:,0] >= self.vectorNx.min()) & (pts[:,0] <= self.vectorNx.max()) + if self.dim > 1: + inside = inside & ((pts[:,1] >= self.vectorNy.min()) & (pts[:,1] <= self.vectorNy.max())) + if self.dim > 2: + inside = inside & ((pts[:,2] >= self.vectorNz.min()) & (pts[:,2] <= self.vectorNz.max())) + return inside + + def getInterpolationMat(self, loc, locType): + """ Produces interpolation matrix + + :param numpy.ndarray loc: Location of points to interpolate to + :param str locType: What to interpolate (see below) + :rtype: scipy.sparse.csr.csr_matrix + :return: M, the interpolation matrix + + locType can be:: + + 'ex' -> x-component of field defined on edges + 'ey' -> y-component of field defined on edges + 'ez' -> z-component of field defined on edges + 'fx' -> x-component of field defined on edges + 'fy' -> y-component of field defined on edges + 'fz' -> z-component of field defined on edges + 'n' -> scalar field defined on nodes + 'cc' -> scalar field defined on cell centres + """ + + loc = np.atleast_2d(loc) + assert np.all(self.isInside(loc)), "Points outside of mesh" + + if self.dim == 3: + + if locType == 'fx': + Qx = interpmat(self.vectorNx, + self.vectorCCy, + self.vectorCCz, + loc[:,0], loc[:,1], loc[:,2]) + Qy = spzeros(loc.shape[0], self.nF[1]) + Qz = spzeros(loc.shape[0], self.nF[2]) + Q = sp.hstack([Qx, Qy, Qz]) + elif locType == 'fy': + Qx = spzeros(loc.shape[0], self.nF[0]) + Qy = interpmat(self.vectorCCx, + self.vectorNy, + self.vectorCCz, + loc[:,0], loc[:,1], loc[:,2]) + Qz = spzeros(loc.shape[0], self.nF[2]) + Q = sp.hstack([Qx, Qy, Qz]) + elif locType == 'fz': + Qx = spzeros(loc.shape[0], self.nF[0]) + Qy = spzeros(loc.shape[0], self.nF[1]) + Qz = interpmat(self.vectorCCx, + self.vectorCCy, + self.vectorNz, + loc[:,0], loc[:,1], loc[:,2]) + Q = sp.hstack([Qx, Qy, Qz]) + elif locType == 'ex': + Qx = interpmat(self.vectorCCx, + self.vectorNy, + self.vectorNz, + loc[:,0], loc[:,1], loc[:,2]) + Qy = spzeros(loc.shape[0], self.nF[1]) + Qz = spzeros(loc.shape[0], self.nF[2]) + Q = sp.hstack([Qx, Qy, Qz]) + elif locType == 'ey': + Qx = spzeros(loc.shape[0], self.nF[0]) + Qy = interpmat(self.vectorNx, + self.vectorCCy, + self.vectorNz, + loc[:,0], loc[:,1], loc[:,2]) + Qz = spzeros(loc.shape[0], self.nF[2]) + Q = sp.hstack([Qx, Qy, Qz]) + elif locType == 'ez': + Qx = spzeros(loc.shape[0], self.nF[0]) + Qy = spzeros(loc.shape[0], self.nF[1]) + Qz = interpmat(self.vectorNx, + self.vectorNy, + self.vectorCCz, + loc[:,0], loc[:,1], loc[:,2]) + Q = sp.hstack([Qx, Qy, Qz]) + elif locType == 'n': + Q = interpmat(self.vectorNx, + self.vectorNy, + self.vectorNz, + loc[:,0], loc[:,1], loc[:,2]) + elif locType == 'cc': + Q = interpmat(self.vectorCCx, + self.vectorCCy, + self.vectorCCz, + loc[:,0], loc[:,1], loc[:,2]) + else: + raise NotImplementedError('getInterpolationMat: locType=='+locType) + else: + raise NotImplementedError('getInterpolationMat: dim=='+str(m.dim)) + return Q if __name__ == '__main__': print('Welcome to tensor mesh!') diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index 7c976ce3..9457df03 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -1,7 +1,9 @@ import matutils import sputils import lomutils +import interputils import ModelBuilder from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind from sputils import spzeros, kron3, speye, sdiag from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird +from interputils import interpmat \ No newline at end of file diff --git a/SimPEG/utils/interputils.py b/SimPEG/utils/interputils.py new file mode 100644 index 00000000..0d528e93 --- /dev/null +++ b/SimPEG/utils/interputils.py @@ -0,0 +1,77 @@ +import numpy as np +import scipy.sparse as sp +from sputils import spzeros +from matutils import mkvc, sub2ind + +def interpmat(x,y,z,xr,yr,zr): + + """ Local interpolation computed for each receiver point in turn """ + + nx = max(x.shape) + ny = max(y.shape) + nz = max(z.shape) + npts = max(xr.shape) + + Q = sp.lil_matrix((npts, nx*ny*nz)) + + for i in range(npts): + # in x-direction + im = np.argmin(abs(x-xr[i])) + if xr[i] - x[im] >= 0: # Point on the left + ind_x1 = im + ind_x2 = im+1 + elif xr[i] - x[im] < 0: # Point on the right + ind_x1 = im-1 + ind_x2 = im + dx1 = xr[i] - x[ind_x1] + dx2 = x[ind_x2] - xr[i] + # in y-direction + im = np.argmin(abs(y-yr[i])) + if yr[i] - y[im] >= 0: # Point on the left + ind_y1 = im + ind_y2 = im+1 + elif yr[i] - y[im] < 0: # Point on the right + ind_y1 = im-1 + ind_y2 = im + dy1 = yr[i] - y[ind_y1] + dy2 = y[ind_y2] - yr[i] + # in z-direction + im = np.argmin(abs(z-zr[i])) + if zr[i] - z[im] >= 0: # Point on the left + ind_z1 = im + ind_z2 = im+1 + elif zr[i] - z[im] < 0: # Point on the right + ind_z1 = im-1 + ind_z2 = im + dz1 = zr[i] - z[ind_z1] + dz2 = z[ind_z2] - zr[i] + dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1]) + + Dx = x[ind_x2] - x[ind_x1] + Dy = y[ind_y2] - y[ind_y1] + Dz = z[ind_z2] - z[ind_z1] + + # Get the row in the matrix + + inds = sub2ind((nx,ny,nz),[ + ( ind_x1, ind_y2, ind_z1), + ( ind_x1, ind_y1, ind_z1), + ( ind_x2, ind_y1, ind_z1), + ( ind_x2, ind_y2, ind_z1), + ( ind_x1, ind_y1, ind_z2), + ( ind_x1, ind_y2, ind_z2), + ( ind_x2, ind_y1, ind_z2), + ( ind_x2, ind_y2, ind_z2)]) + + vals = [(1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz), + (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz), + (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz), + (1-dx2/Dx)*(1-dy2/Dy)*(1-dz1/Dz), + (1-dx1/Dx)*(1-dy1/Dy)*(1-dz2/Dz), + (1-dx1/Dx)*(1-dy2/Dy)*(1-dz2/Dz), + (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz), + (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz)] + + Q[i, mkvc(inds)] = vals + Q = Q.tocsr() + return Q \ No newline at end of file From 6f141ecaf299e91d973244f173e97950996cf927 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 15:11:10 -0800 Subject: [PATCH 35/45] 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() From 4dbfc77aba4018dfb2652cbc67807b9eda898208 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 15:24:35 -0800 Subject: [PATCH 36/45] changed these to be consistent with other places in the code. --- SimPEG/mesh/TensorMesh.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 1dbd62a3..4a37d7de 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -342,22 +342,22 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): locType can be:: - 'ex' -> x-component of field defined on edges - 'ey' -> y-component of field defined on edges - 'ez' -> z-component of field defined on edges - 'fx' -> x-component of field defined on edges - 'fy' -> y-component of field defined on edges - 'fz' -> z-component of field defined on edges - 'n' -> scalar field defined on nodes - 'cc' -> scalar field defined on cell centres + 'Ex' -> x-component of field defined on edges + 'Ey' -> y-component of field defined on edges + 'Ez' -> z-component of field defined on edges + 'Fx' -> x-component of field defined on edges + 'Fy' -> y-component of field defined on edges + 'Fz' -> z-component of field defined on edges + 'N' -> scalar field defined on nodes + 'CC' -> scalar field defined on cell centers """ loc = np.atleast_2d(loc) assert np.all(self.isInside(loc)), "Points outside of mesh" if self.dim == 3: - - if locType == 'fx': + + if locType == 'Fx': Qx = interpmat(self.vectorNx, self.vectorCCy, self.vectorCCz, @@ -365,7 +365,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): Qy = spzeros(loc.shape[0], self.nF[1]) Qz = spzeros(loc.shape[0], self.nF[2]) Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'fy': + elif locType == 'Fy': Qx = spzeros(loc.shape[0], self.nF[0]) Qy = interpmat(self.vectorCCx, self.vectorNy, @@ -373,7 +373,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): loc[:,0], loc[:,1], loc[:,2]) Qz = spzeros(loc.shape[0], self.nF[2]) Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'fz': + elif locType == 'Fz': Qx = spzeros(loc.shape[0], self.nF[0]) Qy = spzeros(loc.shape[0], self.nF[1]) Qz = interpmat(self.vectorCCx, @@ -381,7 +381,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): self.vectorNz, loc[:,0], loc[:,1], loc[:,2]) Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'ex': + elif locType == 'Ex': Qx = interpmat(self.vectorCCx, self.vectorNy, self.vectorNz, @@ -389,15 +389,15 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): Qy = spzeros(loc.shape[0], self.nF[1]) Qz = spzeros(loc.shape[0], self.nF[2]) Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'ey': + elif locType == 'Ey': Qx = spzeros(loc.shape[0], self.nF[0]) Qy = interpmat(self.vectorNx, self.vectorCCy, self.vectorNz, loc[:,0], loc[:,1], loc[:,2]) - Qz = spzeros(loc.shape[0], self.nF[2]) + Qz = spzeros(loc.shape[0], self.nF[2]) Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'ez': + elif locType == 'Ez': Qx = spzeros(loc.shape[0], self.nF[0]) Qy = spzeros(loc.shape[0], self.nF[1]) Qz = interpmat(self.vectorNx, @@ -405,12 +405,12 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): self.vectorCCz, loc[:,0], loc[:,1], loc[:,2]) Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'n': + elif locType == 'N': Q = interpmat(self.vectorNx, self.vectorNy, self.vectorNz, loc[:,0], loc[:,1], loc[:,2]) - elif locType == 'cc': + elif locType == 'CC': Q = interpmat(self.vectorCCx, self.vectorCCy, self.vectorCCz, From df329c49550c47ae12fdc508b0edc5764651398a Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 15:33:52 -0800 Subject: [PATCH 37/45] Bug fix for the Edge vars. --- SimPEG/mesh/TensorMesh.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 4a37d7de..cee23f7e 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -386,20 +386,20 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): self.vectorNy, self.vectorNz, loc[:,0], loc[:,1], loc[:,2]) - Qy = spzeros(loc.shape[0], self.nF[1]) - Qz = spzeros(loc.shape[0], self.nF[2]) + Qy = spzeros(loc.shape[0], self.nE[1]) + Qz = spzeros(loc.shape[0], self.nE[2]) Q = sp.hstack([Qx, Qy, Qz]) elif locType == 'Ey': - Qx = spzeros(loc.shape[0], self.nF[0]) + Qx = spzeros(loc.shape[0], self.nE[0]) Qy = interpmat(self.vectorNx, self.vectorCCy, self.vectorNz, loc[:,0], loc[:,1], loc[:,2]) - Qz = spzeros(loc.shape[0], self.nF[2]) + Qz = spzeros(loc.shape[0], self.nE[2]) Q = sp.hstack([Qx, Qy, Qz]) elif locType == 'Ez': - Qx = spzeros(loc.shape[0], self.nF[0]) - Qy = spzeros(loc.shape[0], self.nF[1]) + Qx = spzeros(loc.shape[0], self.nE[0]) + Qy = spzeros(loc.shape[0], self.nE[1]) Qz = interpmat(self.vectorNx, self.vectorNy, self.vectorCCz, From 592490844d1b4344a7587028e12cfc777c808bd4 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 16:05:16 -0800 Subject: [PATCH 38/45] Testing Functions for interpolation. --- SimPEG/tests/test_interpolation.py | 81 ++++++++++++++++++++++++++++++ SimPEG/utils/interputils.py | 26 +++++----- 2 files changed, 94 insertions(+), 13 deletions(-) create mode 100644 SimPEG/tests/test_interpolation.py diff --git a/SimPEG/tests/test_interpolation.py b/SimPEG/tests/test_interpolation.py new file mode 100644 index 00000000..d98d58d6 --- /dev/null +++ b/SimPEG/tests/test_interpolation.py @@ -0,0 +1,81 @@ +import numpy as np +import unittest +from TestUtils import OrderTest + +MESHTYPES = ['uniformTensorMesh'] +call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) +call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) +cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] +cart_row3 = lambda g, xfun, yfun, zfun: np.c_[call3(xfun, g), call3(yfun, g), call3(zfun, g)] +cartF2 = lambda M, fx, fy: np.vstack((cart_row2(M.gridFx, fx, fy), cart_row2(M.gridFy, fx, fy))) +cartE2 = lambda M, ex, ey: np.vstack((cart_row2(M.gridEx, ex, ey), cart_row2(M.gridEy, ex, ey))) +cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_row3(M.gridFy, fx, fy, fz), cart_row3(M.gridFz, fx, fy, fz))) +cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez))) + + +LOCS = np.random.rand(50,3)*0.6+0.2 + +class TestInterpolation(OrderTest): + name = "Interpolation" + meshTypes = MESHTYPES + meshDimension = 3 + + def getError(self): + funX = lambda x, y, z: np.cos(2*np.pi*y) + funY = lambda x, y, z: np.cos(2*np.pi*z) + funZ = lambda x, y, z: np.cos(2*np.pi*x) + + if 'x' in self.type: + anal = call3(funX, LOCS) + elif 'y' in self.type: + anal = call3(funY, LOCS) + elif 'z' in self.type: + anal = call3(funZ, LOCS) + + if 'F' in self.type: + Fc = cartF3(self.M, funX, funY, funZ) + grid = self.M.projectFaceVector(Fc) + elif 'E' in self.type: + Ec = cartE3(self.M, funX, funY, funZ) + grid = self.M.projectEdgeVector(Ec) + + comp = self.M.getInterpolationMat(LOCS, self.type)*grid + + err = np.linalg.norm((comp - anal), np.inf) + return err + + def test_orderFx(self): + self.type = 'Fx' + self.name = 'Interpolation Fx' + self.orderTest() + + def test_orderFy(self): + self.type = 'Fy' + self.name = 'Interpolation Fy' + self.orderTest() + + def test_orderFz(self): + self.type = 'Fz' + self.name = 'Interpolation Fz' + self.orderTest() + + def test_orderEx(self): + self.type = 'Ex' + self.name = 'Interpolation Ex' + self.orderTest() + + def test_orderEy(self): + self.type = 'Ey' + self.name = 'Interpolation Ey' + self.orderTest() + + def test_orderEz(self): + self.type = 'Ez' + self.name = 'Interpolation Ez' + self.orderTest() + + + + +if __name__ == '__main__': + unittest.main() diff --git a/SimPEG/utils/interputils.py b/SimPEG/utils/interputils.py index 0d528e93..12eb8c47 100644 --- a/SimPEG/utils/interputils.py +++ b/SimPEG/utils/interputils.py @@ -7,22 +7,22 @@ def interpmat(x,y,z,xr,yr,zr): """ Local interpolation computed for each receiver point in turn """ - nx = max(x.shape) - ny = max(y.shape) - nz = max(z.shape) - npts = max(xr.shape) + nx = x.size + ny = y.size + nz = z.size + npts = xr.shape[0] Q = sp.lil_matrix((npts, nx*ny*nz)) for i in range(npts): - # in x-direction + # in x-direction im = np.argmin(abs(x-xr[i])) if xr[i] - x[im] >= 0: # Point on the left ind_x1 = im ind_x2 = im+1 elif xr[i] - x[im] < 0: # Point on the right ind_x1 = im-1 - ind_x2 = im + ind_x2 = im dx1 = xr[i] - x[ind_x1] dx2 = x[ind_x2] - xr[i] # in y-direction @@ -32,9 +32,9 @@ def interpmat(x,y,z,xr,yr,zr): ind_y2 = im+1 elif yr[i] - y[im] < 0: # Point on the right ind_y1 = im-1 - ind_y2 = im + ind_y2 = im dy1 = yr[i] - y[ind_y1] - dy2 = y[ind_y2] - yr[i] + dy2 = y[ind_y2] - yr[i] # in z-direction im = np.argmin(abs(z-zr[i])) if zr[i] - z[im] >= 0: # Point on the left @@ -42,9 +42,9 @@ def interpmat(x,y,z,xr,yr,zr): ind_z2 = im+1 elif zr[i] - z[im] < 0: # Point on the right ind_z1 = im-1 - ind_z2 = im + ind_z2 = im dz1 = zr[i] - z[ind_z1] - dz2 = z[ind_z2] - zr[i] + dz2 = z[ind_z2] - zr[i] dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1]) Dx = x[ind_x2] - x[ind_x1] @@ -53,7 +53,7 @@ def interpmat(x,y,z,xr,yr,zr): # Get the row in the matrix - inds = sub2ind((nx,ny,nz),[ + inds = sub2ind((nx,ny,nz),[ ( ind_x1, ind_y2, ind_z1), ( ind_x1, ind_y1, ind_z1), ( ind_x2, ind_y1, ind_z1), @@ -62,7 +62,7 @@ def interpmat(x,y,z,xr,yr,zr): ( ind_x1, ind_y2, ind_z2), ( ind_x2, ind_y1, ind_z2), ( ind_x2, ind_y2, ind_z2)]) - + vals = [(1-dx1/Dx)*(1-dy2/Dy)*(1-dz1/Dz), (1-dx1/Dx)*(1-dy1/Dy)*(1-dz1/Dz), (1-dx2/Dx)*(1-dy1/Dy)*(1-dz1/Dz), @@ -74,4 +74,4 @@ def interpmat(x,y,z,xr,yr,zr): Q[i, mkvc(inds)] = vals Q = Q.tocsr() - return Q \ No newline at end of file + return Q From 967eb1216ffa63fad3df8740e83ab2a8ed1c2149 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 16:25:57 -0800 Subject: [PATCH 39/45] Testing N and CC as well. Updates to TestOrder --- SimPEG/tests/TestUtils.py | 10 ++++++---- SimPEG/tests/test_interpolation.py | 21 ++++++++++++++++++++- 2 files changed, 26 insertions(+), 5 deletions(-) diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index 9b2158c4..ee8f90e4 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -67,8 +67,7 @@ class OrderTest(unittest.TestCase): name = "Order Test" expectedOrders = 2. # This can be a list of orders, must be the same length as meshTypes - _expectedOrder = 2. - tolerance = 0.85 + tolerance = 0.85 # This can also be a list, must be the same length as meshTypes meshSizes = [4, 8, 16, 32] meshTypes = ['uniformTensorMesh'] _meshType = meshTypes[0] @@ -124,6 +123,8 @@ class OrderTest(unittest.TestCase): """ assert type(self.meshTypes) == list, 'meshTypes must be a list' + if type(self.tolerance) is not list: + self.tolerance = np.ones(len(self.meshTypes))*self.tolerance # if we just provide one expected order, repeat it for each mesh type if type(self.expectedOrders) == float or type(self.expectedOrders) == int: @@ -134,6 +135,7 @@ class OrderTest(unittest.TestCase): for ii_meshType, meshType in enumerate(self.meshTypes): self._meshType = meshType + self._tolerance = self.tolerance[ii_meshType] self._expectedOrder = self.expectedOrders[ii_meshType] order = [] @@ -144,7 +146,7 @@ class OrderTest(unittest.TestCase): err = self.getError() if ii == 0: print '' - print 'Testing convergence on ' + self.M._meshType + ' of: ' + self.name + print self._meshType + ': ' + self.name print '_____________________________________________' print ' h | error | e(i-1)/e(i) | order' print '~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~~~~|~~~~~~~~~~' @@ -155,7 +157,7 @@ class OrderTest(unittest.TestCase): err_old = err max_h_old = max_h print '---------------------------------------------' - passTest = np.mean(np.array(order)) > self.tolerance*self._expectedOrder + 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)] else: diff --git a/SimPEG/tests/test_interpolation.py b/SimPEG/tests/test_interpolation.py index d98d58d6..a068c441 100644 --- a/SimPEG/tests/test_interpolation.py +++ b/SimPEG/tests/test_interpolation.py @@ -2,7 +2,8 @@ import numpy as np import unittest from TestUtils import OrderTest -MESHTYPES = ['uniformTensorMesh'] +MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh'] +TOLERANCES = [0.9, 0.6] call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] @@ -18,7 +19,9 @@ LOCS = np.random.rand(50,3)*0.6+0.2 class TestInterpolation(OrderTest): name = "Interpolation" meshTypes = MESHTYPES + tolerance = TOLERANCES meshDimension = 3 + meshSizes = [8, 16, 32] def getError(self): funX = lambda x, y, z: np.cos(2*np.pi*y) @@ -31,6 +34,8 @@ class TestInterpolation(OrderTest): anal = call3(funY, LOCS) elif 'z' in self.type: anal = call3(funZ, LOCS) + else: + anal = call3(funX, LOCS) if 'F' in self.type: Fc = cartF3(self.M, funX, funY, funZ) @@ -38,12 +43,26 @@ class TestInterpolation(OrderTest): elif 'E' in self.type: Ec = cartE3(self.M, funX, funY, funZ) grid = self.M.projectEdgeVector(Ec) + elif 'CC' == self.type: + grid = call3(funX, self.M.gridCC) + elif 'N' == self.type: + grid = call3(funX, self.M.gridN) comp = self.M.getInterpolationMat(LOCS, self.type)*grid err = np.linalg.norm((comp - anal), np.inf) return err + def test_orderCC(self): + self.type = 'CC' + self.name = 'Interpolation CC' + self.orderTest() + + def test_orderN(self): + self.type = 'N' + self.name = 'Interpolation N' + self.orderTest() + def test_orderFx(self): self.type = 'Fx' self.name = 'Interpolation Fx' From 4784224a84838d85ea7d91f9c3330f16a05c785a Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 16:35:11 -0800 Subject: [PATCH 40/45] Put repetitive code in a function. --- SimPEG/utils/interputils.py | 42 +++++++++++++------------------------ 1 file changed, 15 insertions(+), 27 deletions(-) diff --git a/SimPEG/utils/interputils.py b/SimPEG/utils/interputils.py index 12eb8c47..0a544532 100644 --- a/SimPEG/utils/interputils.py +++ b/SimPEG/utils/interputils.py @@ -14,37 +14,25 @@ def interpmat(x,y,z,xr,yr,zr): Q = sp.lil_matrix((npts, nx*ny*nz)) - for i in range(npts): - # in x-direction - im = np.argmin(abs(x-xr[i])) - if xr[i] - x[im] >= 0: # Point on the left + + def inter1D(x, xr_i): + im = np.argmin(abs(x-xr_i)) + if xr_i - x[im] >= 0: # Point on the left ind_x1 = im ind_x2 = im+1 - elif xr[i] - x[im] < 0: # Point on the right + elif xr_i - x[im] < 0: # Point on the right ind_x1 = im-1 ind_x2 = im - dx1 = xr[i] - x[ind_x1] - dx2 = x[ind_x2] - xr[i] - # in y-direction - im = np.argmin(abs(y-yr[i])) - if yr[i] - y[im] >= 0: # Point on the left - ind_y1 = im - ind_y2 = im+1 - elif yr[i] - y[im] < 0: # Point on the right - ind_y1 = im-1 - ind_y2 = im - dy1 = yr[i] - y[ind_y1] - dy2 = y[ind_y2] - yr[i] - # in z-direction - im = np.argmin(abs(z-zr[i])) - if zr[i] - z[im] >= 0: # Point on the left - ind_z1 = im - ind_z2 = im+1 - elif zr[i] - z[im] < 0: # Point on the right - ind_z1 = im-1 - ind_z2 = im - dz1 = zr[i] - z[ind_z1] - dz2 = z[ind_z2] - zr[i] + dx1 = xr_i - x[ind_x1] + dx2 = x[ind_x2] - xr_i + return ind_x1, ind_x2, dx1, dx2 + + for i in range(npts): + # in x-direction + ind_x1, ind_x2, dx1, dx2 = inter1D(x, xr[i]) + ind_y1, ind_y2, dy1, dy2 = inter1D(y, yr[i]) + ind_z1, ind_z2, dz1, dz2 = inter1D(z, zr[i]) + dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1]) Dx = x[ind_x2] - x[ind_x1] From 51a539a291181a78e5d7cf48a1e2efdaa08e3695 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 4 Nov 2013 18:17:01 -0800 Subject: [PATCH 41/45] Generalized to any dimension. Tested. --- SimPEG/mesh/TensorMesh.py | 117 ++++++++++++--------------- SimPEG/tests/test_interpolation.py | 122 ++++++++++++++++++++++++++--- SimPEG/utils/interputils.py | 102 +++++++++++++++++++----- 3 files changed, 244 insertions(+), 97 deletions(-) diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index ea3a0200..87b81303 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -315,6 +315,45 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): # --------------- Methods --------------------- + def getTensor(self, locType): + """ Returns a tensor list. + + :param str locType: What tensor (see below) + :rtype: list + :return: list of the tensors that make up the mesh. + + locType can be:: + + 'Ex' -> x-component of field defined on edges + 'Ey' -> y-component of field defined on edges + 'Ez' -> z-component of field defined on edges + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces + 'N' -> scalar field defined on nodes + 'CC' -> scalar field defined on cell centers + """ + + if locType is 'Fx': + ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] + elif locType is 'Fy': + ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] + elif locType is 'Fz': + ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] + elif locType is 'Ex': + ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] + elif locType is 'Ey': + ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] + elif locType is 'Ez': + ten = [self.vectorNx , self.vectorNy , self.vectorCCz] + elif locType is 'CC': + ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] + elif locType is 'N': + ten = [self.vectorNx , self.vectorNy , self.vectorNz ] + + return [t for t in ten if t is not None] + + def isInside(self, pts): """ Determines if a set of points are inside a mesh. @@ -345,9 +384,9 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): 'Ex' -> x-component of field defined on edges 'Ey' -> y-component of field defined on edges 'Ez' -> z-component of field defined on edges - 'Fx' -> x-component of field defined on edges - 'Fy' -> y-component of field defined on edges - 'Fz' -> z-component of field defined on edges + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces 'N' -> scalar field defined on nodes 'CC' -> scalar field defined on cell centers """ @@ -355,70 +394,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): loc = np.atleast_2d(loc) assert np.all(self.isInside(loc)), "Points outside of mesh" - if self.dim == 3: - - if locType == 'Fx': - Qx = interpmat(self.vectorNx, - self.vectorCCy, - self.vectorCCz, - loc[:,0], loc[:,1], loc[:,2]) - Qy = spzeros(loc.shape[0], self.nF[1]) - Qz = spzeros(loc.shape[0], self.nF[2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'Fy': - Qx = spzeros(loc.shape[0], self.nF[0]) - Qy = interpmat(self.vectorCCx, - self.vectorNy, - self.vectorCCz, - loc[:,0], loc[:,1], loc[:,2]) - Qz = spzeros(loc.shape[0], self.nF[2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'Fz': - Qx = spzeros(loc.shape[0], self.nF[0]) - Qy = spzeros(loc.shape[0], self.nF[1]) - Qz = interpmat(self.vectorCCx, - self.vectorCCy, - self.vectorNz, - loc[:,0], loc[:,1], loc[:,2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'Ex': - Qx = interpmat(self.vectorCCx, - self.vectorNy, - self.vectorNz, - loc[:,0], loc[:,1], loc[:,2]) - Qy = spzeros(loc.shape[0], self.nE[1]) - Qz = spzeros(loc.shape[0], self.nE[2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'Ey': - Qx = spzeros(loc.shape[0], self.nE[0]) - Qy = interpmat(self.vectorNx, - self.vectorCCy, - self.vectorNz, - loc[:,0], loc[:,1], loc[:,2]) - Qz = spzeros(loc.shape[0], self.nE[2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'Ez': - Qx = spzeros(loc.shape[0], self.nE[0]) - Qy = spzeros(loc.shape[0], self.nE[1]) - Qz = interpmat(self.vectorNx, - self.vectorNy, - self.vectorCCz, - loc[:,0], loc[:,1], loc[:,2]) - Q = sp.hstack([Qx, Qy, Qz]) - elif locType == 'N': - Q = interpmat(self.vectorNx, - self.vectorNy, - self.vectorNz, - loc[:,0], loc[:,1], loc[:,2]) - elif locType == 'CC': - Q = interpmat(self.vectorCCx, - self.vectorCCy, - self.vectorCCz, - loc[:,0], loc[:,1], loc[:,2]) - else: - raise NotImplementedError('getInterpolationMat: locType=='+locType) + ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1 + if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind: + nF_nE = self.nF if 'F' in locType else self.nE + components = [spzeros(loc.shape[0], n) for n in nF_nE] + components[ind] = interpmat(loc, *self.getTensor(locType)) + Q = sp.hstack(components) + elif locType in ['CC', 'N']: + Q = interpmat(loc, *self.getTensor(locType)) else: - raise NotImplementedError('getInterpolationMat: dim=='+str(m.dim)) + raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) return Q if __name__ == '__main__': diff --git a/SimPEG/tests/test_interpolation.py b/SimPEG/tests/test_interpolation.py index a068c441..6931aa42 100644 --- a/SimPEG/tests/test_interpolation.py +++ b/SimPEG/tests/test_interpolation.py @@ -1,9 +1,11 @@ import numpy as np import unittest from TestUtils import OrderTest +from SimPEG.utils import mkvc MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh'] -TOLERANCES = [0.9, 0.6] +TOLERANCES = [0.9, 0.55] +call1 = lambda fun, xyz: fun(xyz) call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2]) cart_row2 = lambda g, xfun, yfun: np.c_[call2(xfun, g), call2(yfun, g)] @@ -14,14 +16,114 @@ cartF3 = lambda M, fx, fy, fz: np.vstack((cart_row3(M.gridFx, fx, fy, fz), cart_ cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_row3(M.gridEy, ex, ey, ez), cart_row3(M.gridEz, ex, ey, ez))) -LOCS = np.random.rand(50,3)*0.6+0.2 -class TestInterpolation(OrderTest): +class TestInterpolation1D(OrderTest): + LOCS = np.random.rand(50,1)*0.6+0.2 + name = "Interpolation 1D" + meshTypes = MESHTYPES + tolerance = TOLERANCES + meshDimension = 1 + meshSizes = [8, 16, 32] + + def getError(self): + funX = lambda x: np.cos(2*np.pi*x) + + anal = mkvc(call1(funX, self.LOCS)) + + if 'CC' == self.type: + grid = call1(funX, self.M.gridCC) + elif 'N' == self.type: + grid = call1(funX, self.M.gridN) + + comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid + + err = np.linalg.norm((comp - anal), 2) + return err + + def test_orderCC(self): + self.type = 'CC' + self.name = 'Interpolation 1D: CC' + self.orderTest() + + def test_orderN(self): + self.type = 'N' + self.name = 'Interpolation 1D: N' + self.orderTest() + +class TestInterpolation2d(OrderTest): + name = "Interpolation 2D" + LOCS = np.random.rand(50,2)*0.6+0.2 + meshTypes = MESHTYPES + tolerance = TOLERANCES + meshDimension = 2 + meshSizes = [8, 16, 32, 64] + + def getError(self): + funX = lambda x, y: np.cos(2*np.pi*y) + funY = lambda x, y: np.cos(2*np.pi*x) + + if 'x' in self.type: + anal = call2(funX, self.LOCS) + elif 'y' in self.type: + anal = call2(funY, self.LOCS) + else: + anal = call2(funX, self.LOCS) + + if 'F' in self.type: + Fc = cartF2(self.M, funX, funY) + grid = self.M.projectFaceVector(Fc) + elif 'E' in self.type: + Ec = cartE2(self.M, funX, funY) + grid = self.M.projectEdgeVector(Ec) + elif 'CC' == self.type: + grid = call2(funX, self.M.gridCC) + elif 'N' == self.type: + grid = call2(funX, self.M.gridN) + + comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid + + err = np.linalg.norm((comp - anal), np.inf) + return err + + def test_orderCC(self): + self.type = 'CC' + self.name = 'Interpolation 2D: CC' + self.orderTest() + + def test_orderN(self): + self.type = 'N' + self.name = 'Interpolation 2D: N' + self.orderTest() + + def test_orderFx(self): + self.type = 'Fx' + self.name = 'Interpolation 2D: Fx' + self.orderTest() + + def test_orderFy(self): + self.type = 'Fy' + self.name = 'Interpolation 2D: Fy' + self.orderTest() + + def test_orderEx(self): + self.type = 'Ex' + self.name = 'Interpolation 2D: Ex' + self.orderTest() + + def test_orderEy(self): + self.type = 'Ey' + self.name = 'Interpolation 2D: Ey' + self.orderTest() + + + +class TestInterpolation3D(OrderTest): name = "Interpolation" + LOCS = np.random.rand(50,3)*0.6+0.2 meshTypes = MESHTYPES tolerance = TOLERANCES meshDimension = 3 - meshSizes = [8, 16, 32] + meshSizes = [8, 16, 32, 64] def getError(self): funX = lambda x, y, z: np.cos(2*np.pi*y) @@ -29,13 +131,13 @@ class TestInterpolation(OrderTest): funZ = lambda x, y, z: np.cos(2*np.pi*x) if 'x' in self.type: - anal = call3(funX, LOCS) + anal = call3(funX, self.LOCS) elif 'y' in self.type: - anal = call3(funY, LOCS) + anal = call3(funY, self.LOCS) elif 'z' in self.type: - anal = call3(funZ, LOCS) + anal = call3(funZ, self.LOCS) else: - anal = call3(funX, LOCS) + anal = call3(funX, self.LOCS) if 'F' in self.type: Fc = cartF3(self.M, funX, funY, funZ) @@ -48,7 +150,7 @@ class TestInterpolation(OrderTest): elif 'N' == self.type: grid = call3(funX, self.M.gridN) - comp = self.M.getInterpolationMat(LOCS, self.type)*grid + comp = self.M.getInterpolationMat(self.LOCS, self.type)*grid err = np.linalg.norm((comp - anal), np.inf) return err @@ -94,7 +196,5 @@ class TestInterpolation(OrderTest): self.orderTest() - - if __name__ == '__main__': unittest.main() diff --git a/SimPEG/utils/interputils.py b/SimPEG/utils/interputils.py index 0a544532..0ece29c7 100644 --- a/SimPEG/utils/interputils.py +++ b/SimPEG/utils/interputils.py @@ -3,35 +3,97 @@ import scipy.sparse as sp from sputils import spzeros from matutils import mkvc, sub2ind -def interpmat(x,y,z,xr,yr,zr): +def _interp_point_1D(x, xr_i): + im = np.argmin(abs(x-xr_i)) + if xr_i - x[im] >= 0: # Point on the left + ind_x1 = im + ind_x2 = im+1 + elif xr_i - x[im] < 0: # Point on the right + ind_x1 = im-1 + ind_x2 = im + dx1 = xr_i - x[ind_x1] + dx2 = x[ind_x2] - xr_i + return ind_x1, ind_x2, dx1, dx2 + +def interpmat(locs, x, y=None, z=None): """ Local interpolation computed for each receiver point in turn """ + if y is None and z is None: + return interpmat1D(locs, x) + elif z is None: + return interpmat2D(locs, x, y) + else: + return interpmat3D(locs, x, y, z) + +def interpmat1D(locs, x): + nx = x.size + locs = mkvc(locs) + npts = locs.shape[0] + + Q = sp.lil_matrix((npts, nx)) + + for i in range(npts): + ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i]) + dv = (x[ind_x2] - x[ind_x1]) + Dx = x[ind_x2] - x[ind_x1] + # Get the row in the matrix + inds = [ind_x1, ind_x2] + vals = [(1-dx1/Dx),(1-dx2/Dx)] + Q[i, inds] = vals + return Q.tocsr() + + + +def interpmat2D(locs, x, y): + nx = x.size + ny = y.size + npts = locs.shape[0] + + Q = sp.lil_matrix((npts, nx*ny)) + + + for i in range(npts): + ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0]) + ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1]) + + dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) + + Dx = x[ind_x2] - x[ind_x1] + Dy = y[ind_y2] - y[ind_y1] + + # Get the row in the matrix + + inds = sub2ind((nx,ny),[ + ( ind_x1, ind_y2), + ( ind_x1, ind_y1), + ( ind_x2, ind_y1), + ( ind_x2, ind_y2)]) + + vals = [(1-dx1/Dx)*(1-dy2/Dy), + (1-dx1/Dx)*(1-dy1/Dy), + (1-dx2/Dx)*(1-dy1/Dy), + (1-dx2/Dx)*(1-dy2/Dy)] + + Q[i, mkvc(inds)] = vals + + return Q.tocsr() + + + +def interpmat3D(locs, x, y, z): nx = x.size ny = y.size nz = z.size - npts = xr.shape[0] + npts = locs.shape[0] Q = sp.lil_matrix((npts, nx*ny*nz)) - def inter1D(x, xr_i): - im = np.argmin(abs(x-xr_i)) - if xr_i - x[im] >= 0: # Point on the left - ind_x1 = im - ind_x2 = im+1 - elif xr_i - x[im] < 0: # Point on the right - ind_x1 = im-1 - ind_x2 = im - dx1 = xr_i - x[ind_x1] - dx2 = x[ind_x2] - xr_i - return ind_x1, ind_x2, dx1, dx2 - for i in range(npts): - # in x-direction - ind_x1, ind_x2, dx1, dx2 = inter1D(x, xr[i]) - ind_y1, ind_y2, dy1, dy2 = inter1D(y, yr[i]) - ind_z1, ind_z2, dz1, dz2 = inter1D(z, zr[i]) + ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0]) + ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1]) + ind_z1, ind_z2, dz1, dz2 = _interp_point_1D(z, locs[i, 2]) dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1]) @@ -61,5 +123,5 @@ def interpmat(x,y,z,xr,yr,zr): (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz)] Q[i, mkvc(inds)] = vals - Q = Q.tocsr() - return Q + + return Q.tocsr() From 9b00617e3fcb9e9401f03b6a6c9b899a271e546f Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 5 Nov 2013 10:32:43 -0800 Subject: [PATCH 42/45] Documentation updates. --- SimPEG/tests/test_interpolation.py | 4 +- SimPEG/utils/interputils.py | 68 +++++++++++++++++++++++++++--- docs/api_Utils.rst | 4 ++ 3 files changed, 67 insertions(+), 9 deletions(-) diff --git a/SimPEG/tests/test_interpolation.py b/SimPEG/tests/test_interpolation.py index 6931aa42..3d751dde 100644 --- a/SimPEG/tests/test_interpolation.py +++ b/SimPEG/tests/test_interpolation.py @@ -18,7 +18,7 @@ cartE3 = lambda M, ex, ey, ez: np.vstack((cart_row3(M.gridEx, ex, ey, ez), cart_ class TestInterpolation1D(OrderTest): - LOCS = np.random.rand(50,1)*0.6+0.2 + LOCS = np.random.rand(50)*0.6+0.2 name = "Interpolation 1D" meshTypes = MESHTYPES tolerance = TOLERANCES @@ -28,7 +28,7 @@ class TestInterpolation1D(OrderTest): def getError(self): funX = lambda x: np.cos(2*np.pi*x) - anal = mkvc(call1(funX, self.LOCS)) + anal = call1(funX, self.LOCS) if 'CC' == self.type: grid = call1(funX, self.M.gridCC) diff --git a/SimPEG/utils/interputils.py b/SimPEG/utils/interputils.py index 0ece29c7..c8bcb4ec 100644 --- a/SimPEG/utils/interputils.py +++ b/SimPEG/utils/interputils.py @@ -4,6 +4,15 @@ from sputils import spzeros from matutils import mkvc, sub2ind def _interp_point_1D(x, xr_i): + """ + given a point, xr_i, this will find which two integers it lies between. + + :param numpy.ndarray x: Tensor vector of 1st dimension of grid. + :param float xr_i: Location of a point + :rtype: int,int,float,float + :return: index1, index2, portion1, portion2 + """ + # TODO: This fails if the point is on the outside of the mesh. We may want to replace this by extrapolation? im = np.argmin(abs(x-xr_i)) if xr_i - x[im] >= 0: # Point on the left ind_x1 = im @@ -17,16 +26,43 @@ def _interp_point_1D(x, xr_i): def interpmat(locs, x, y=None, z=None): - """ Local interpolation computed for each receiver point in turn """ + """ + Local interpolation computed for each receiver point in turn + + :param numpy.ndarray loc: Location of points to interpolate to + :param numpy.ndarray x: Tensor vector of 1st dimension of grid. + :param numpy.ndarray y: Tensor vector of 2nd dimension of grid. None by default. + :param numpy.ndarray z: Tensor vector of 3rd dimension of grid. None by default. + :rtype: scipy.sparse.csr.csr_matrix + :return: Interpolation matrix + + .. plot:: + + import SimPEG + import numpy as np + import matplotlib.pyplot as plt + locs = np.random.rand(50)*0.8+0.1 + x = np.linspace(0,1,7) + dense = np.linspace(0,1,200) + fun = lambda x: np.cos(2*np.pi*x) + Q = SimPEG.utils.interpmat(locs, x) + plt.plot(x, fun(x), 'bs-') + plt.plot(dense, fun(dense), 'y:') + plt.plot(locs, Q*fun(x), 'mo') + plt.plot(locs, fun(locs), 'rx') + plt.show() + + """ if y is None and z is None: - return interpmat1D(locs, x) + return _interpmat1D(locs, x) elif z is None: - return interpmat2D(locs, x, y) + return _interpmat2D(locs, x, y) else: - return interpmat3D(locs, x, y, z) + return _interpmat3D(locs, x, y, z) -def interpmat1D(locs, x): +def _interpmat1D(locs, x): + """Use interpmat with only x component provided.""" nx = x.size locs = mkvc(locs) npts = locs.shape[0] @@ -45,7 +81,8 @@ def interpmat1D(locs, x): -def interpmat2D(locs, x, y): +def _interpmat2D(locs, x, y): + """Use interpmat with only x and y components provided.""" nx = x.size ny = y.size npts = locs.shape[0] @@ -81,7 +118,8 @@ def interpmat2D(locs, x, y): -def interpmat3D(locs, x, y, z): +def _interpmat3D(locs, x, y, z): + """Use interpmat.""" nx = x.size ny = y.size nz = z.size @@ -125,3 +163,19 @@ def interpmat3D(locs, x, y, z): Q[i, mkvc(inds)] = vals return Q.tocsr() + + +if __name__ == '__main__': + import SimPEG + import numpy as np + import matplotlib.pyplot as plt + locs = np.random.rand(50)*0.8+0.1 + x = np.linspace(0,1,7) + dense = np.linspace(0,1,200) + fun = lambda x: np.cos(2*np.pi*x) + Q = SimPEG.utils.interpmat(locs, x) + plt.plot(x, fun(x), 'bs-') + plt.plot(dense, fun(dense), 'y:') + plt.plot(locs, Q*fun(x), 'mo') + plt.plot(locs, fun(locs), 'rx') + plt.show() diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst index e952010b..91d87d29 100644 --- a/docs/api_Utils.rst +++ b/docs/api_Utils.rst @@ -19,3 +19,7 @@ Utilities :members: :undoc-members: +.. automodule:: SimPEG.utils.interputils + :members: + :undoc-members: + From 37e10cabaa00a0b90e8ab7b2138c3ab3f8a5e857 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 5 Nov 2013 11:12:00 -0800 Subject: [PATCH 43/45] Minor changes to how grids are calculated. Duplicate code has been removed. --- SimPEG/mesh/TensorMesh.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py index 87b81303..6cc5d8a6 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -157,7 +157,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridCC is None: - self._gridCC = ndgrid([x for x in [self.vectorCCx, self.vectorCCy, self.vectorCCz] if not x is None]) + self._gridCC = ndgrid(self.getTensor('CC')) return self._gridCC return locals() _gridCC = None # Store grid by default @@ -168,7 +168,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridN is None: - self._gridN = ndgrid([x for x in [self.vectorNx, self.vectorNy, self.vectorNz] if not x is None]) + self._gridN = ndgrid(self.getTensor('N')) return self._gridN return locals() _gridN = None # Store grid by default @@ -179,7 +179,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridFx is None: - self._gridFx = ndgrid([x for x in [self.vectorNx, self.vectorCCy, self.vectorCCz] if not x is None]) + self._gridFx = ndgrid(self.getTensor('Fx')) return self._gridFx return locals() _gridFx = None # Store grid by default @@ -190,7 +190,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridFy is None and self.dim > 1: - self._gridFy = ndgrid([x for x in [self.vectorCCx, self.vectorNy, self.vectorCCz] if not x is None]) + self._gridFy = ndgrid(self.getTensor('Fy')) return self._gridFy return locals() _gridFy = None # Store grid by default @@ -201,7 +201,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridFz is None and self.dim > 2: - self._gridFz = ndgrid([x for x in [self.vectorCCx, self.vectorCCy, self.vectorNz] if not x is None]) + self._gridFz = ndgrid(self.getTensor('Fz')) return self._gridFz return locals() _gridFz = None # Store grid by default @@ -212,7 +212,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridEx is None: - self._gridEx = ndgrid([x for x in [self.vectorCCx, self.vectorNy, self.vectorNz] if not x is None]) + self._gridEx = ndgrid(self.getTensor('Ex')) return self._gridEx return locals() _gridEx = None # Store grid by default @@ -223,7 +223,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridEy is None and self.dim > 1: - self._gridEy = ndgrid([x for x in [self.vectorNx, self.vectorCCy, self.vectorNz] if not x is None]) + self._gridEy = ndgrid(self.getTensor('Ey')) return self._gridEy return locals() _gridEy = None # Store grid by default @@ -234,7 +234,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts): def fget(self): if self._gridEz is None and self.dim > 2: - self._gridEz = ndgrid([x for x in [self.vectorNx, self.vectorNy, self.vectorCCz] if not x is None]) + self._gridEz = ndgrid(self.getTensor('Ez')) return self._gridEz return locals() _gridEz = None # Store grid by default @@ -420,5 +420,6 @@ if __name__ == '__main__': h = h[:testDim] M = TensorMesh(h) + print M xn = M.plotGrid() From ebf8f23abcd16b9a671dd5bfa6c54ec511f1b82d Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 5 Nov 2013 15:43:16 -0800 Subject: [PATCH 44/45] Document the Test results --- SimPEG/tests/HTMLTestRunner.py | 824 +++++++++ SimPEG/tests/api_TestResults.rst | 355 ++++ SimPEG/tests/report.html | 2870 ++++++++++++++++++++++++++++++ SimPEG/tests/runTests.py | 38 +- docs/api_TestResults.rst | 2862 +++++++++++++++++++++++++++++ docs/index.rst | 1 + 6 files changed, 6949 insertions(+), 1 deletion(-) create mode 100644 SimPEG/tests/HTMLTestRunner.py create mode 100644 SimPEG/tests/api_TestResults.rst create mode 100644 SimPEG/tests/report.html create mode 100644 docs/api_TestResults.rst diff --git a/SimPEG/tests/HTMLTestRunner.py b/SimPEG/tests/HTMLTestRunner.py new file mode 100644 index 00000000..af384971 --- /dev/null +++ b/SimPEG/tests/HTMLTestRunner.py @@ -0,0 +1,824 @@ +""" +A TestRunner for use with the Python unit testing framework. It +generates a HTML report to show the result at a glance. + +The simplest way to use this is to invoke its main method. E.g. + + import unittest + import HTMLTestRunner + + ... define your tests ... + + if __name__ == '__main__': + HTMLTestRunner.main() + + +For more customization options, instantiates a HTMLTestRunner object. +HTMLTestRunner is a counterpart to unittest's TextTestRunner. E.g. + + # output to a file + fp = file('my_report.html', 'wb') + runner = HTMLTestRunner.HTMLTestRunner( + stream=fp, + title='My unit test', + description='This demonstrates the report output by HTMLTestRunner.' + ) + + # Use an external stylesheet. + # See the Template_mixin class for more customizable options + runner.STYLESHEET_TMPL = '' + + # run the test + runner.run(my_test_suite) + + +------------------------------------------------------------------------ +Copyright (c) 2004-2007, Wai Yip Tung +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + +* Redistributions of source code must retain the above copyright notice, + this list of conditions and the following disclaimer. +* Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. +* Neither the name Wai Yip Tung nor the names of its contributors may be + used to endorse or promote products derived from this software without + specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +""" + +# URL: http://tungwaiyip.info/software/HTMLTestRunner.html + +__author__ = "Wai Yip Tung" +__version__ = "0.8.2" + + +""" +Change History + +Version 0.8.2 +* Show output inline instead of popup window (Viorel Lupu). + +Version in 0.8.1 +* Validated XHTML (Wolfgang Borgert). +* Added description of test classes and test cases. + +Version in 0.8.0 +* Define Template_mixin class for customization. +* Workaround a IE 6 bug that it does not treat + +%(heading)s +%(report)s +%(ending)s + + + +""" + # variables: (title, generator, stylesheet, heading, report, ending) + + + # ------------------------------------------------------------------------ + # Stylesheet + # + # alternatively use a for external style sheet, e.g. + # + + STYLESHEET_TMPL = """ + +""" + + + + # ------------------------------------------------------------------------ + # Heading + # + + HEADING_TMPL = """
+

%(title)s

+%(parameters)s +

%(description)s

+
+ +""" # variables: (title, parameters, description) + + HEADING_ATTRIBUTE_TMPL = """

%(name)s: %(value)s

+""" # variables: (name, value) + + + + # ------------------------------------------------------------------------ + # Report + # + + REPORT_TMPL = """ +

Show +Summary +Failed +All +

+ ++++++++ + + + + + + + + +%(test_list)s + + + + + + + + +
Test Group/Test caseCountPassFailErrorView
Total%(count)s%(Pass)s%(fail)s%(error)s 
+""" # variables: (test_list, count, Pass, fail, error) + + REPORT_CLASS_TMPL = r""" + + %(desc)s + %(count)s + %(Pass)s + %(fail)s + %(error)s + Detail + +""" # variables: (style, desc, count, Pass, fail, error, cid) + + + REPORT_TEST_WITH_OUTPUT_TMPL = r""" + +
%(desc)s
+ + + + + %(status)s + + + + + + +""" # variables: (tid, Class, style, desc, status) + + + REPORT_TEST_NO_OUTPUT_TMPL = r""" + +
%(desc)s
+ %(status)s + +""" # variables: (tid, Class, style, desc, status) + + + REPORT_TEST_OUTPUT_TMPL = r""" +%(id)s: %(output)s +""" # variables: (id, output) + + + + # ------------------------------------------------------------------------ + # ENDING + # + + ENDING_TMPL = """
 
""" + +# -------------------- The end of the Template class ------------------- + + +TestResult = unittest.TestResult + +class _TestResult(TestResult): + # note: _TestResult is a pure representation of results. + # It lacks the output and reporting ability compares to unittest._TextTestResult. + + def __init__(self, verbosity=1): + TestResult.__init__(self) + self.stdout0 = None + self.stderr0 = None + self.success_count = 0 + self.failure_count = 0 + self.error_count = 0 + self.verbosity = verbosity + + # result is a list of result in 4 tuple + # ( + # result code (0: success; 1: fail; 2: error), + # TestCase object, + # Test output (byte string), + # stack trace, + # ) + self.result = [] + + + def startTest(self, test): + TestResult.startTest(self, test) + # just one buffer for both stdout and stderr + self.outputBuffer = StringIO.StringIO() + stdout_redirector.fp = self.outputBuffer + stderr_redirector.fp = self.outputBuffer + self.stdout0 = sys.stdout + self.stderr0 = sys.stderr + sys.stdout = stdout_redirector + sys.stderr = stderr_redirector + + + def complete_output(self): + """ + Disconnect output redirection and return buffer. + Safe to call multiple times. + """ + if self.stdout0: + sys.stdout = self.stdout0 + sys.stderr = self.stderr0 + self.stdout0 = None + self.stderr0 = None + return self.outputBuffer.getvalue() + + + def stopTest(self, test): + # Usually one of addSuccess, addError or addFailure would have been called. + # But there are some path in unittest that would bypass this. + # We must disconnect stdout in stopTest(), which is guaranteed to be called. + self.complete_output() + + + def addSuccess(self, test): + self.success_count += 1 + TestResult.addSuccess(self, test) + output = self.complete_output() + self.result.append((0, test, output, '')) + if self.verbosity > 1: + sys.stderr.write('ok ') + sys.stderr.write(str(test)) + sys.stderr.write('\n') + else: + sys.stderr.write('.') + + def addError(self, test, err): + self.error_count += 1 + TestResult.addError(self, test, err) + _, _exc_str = self.errors[-1] + output = self.complete_output() + self.result.append((2, test, output, _exc_str)) + if self.verbosity > 1: + sys.stderr.write('E ') + sys.stderr.write(str(test)) + sys.stderr.write('\n') + else: + sys.stderr.write('E') + + def addFailure(self, test, err): + self.failure_count += 1 + TestResult.addFailure(self, test, err) + _, _exc_str = self.failures[-1] + output = self.complete_output() + self.result.append((1, test, output, _exc_str)) + if self.verbosity > 1: + sys.stderr.write('F ') + sys.stderr.write(str(test)) + sys.stderr.write('\n') + else: + sys.stderr.write('F') + + +class HTMLTestRunner(Template_mixin): + """ + """ + def __init__(self, stream=sys.stdout, verbosity=1, title=None, description=None): + self.stream = stream + self.verbosity = verbosity + if title is None: + self.title = self.DEFAULT_TITLE + else: + self.title = title + if description is None: + self.description = self.DEFAULT_DESCRIPTION + else: + self.description = description + + self.startTime = datetime.datetime.now() + + + def run(self, test): + "Run the given test case or test suite." + result = _TestResult(self.verbosity) + test(result) + self.stopTime = datetime.datetime.now() + self.generateReport(test, result) + print >>sys.stderr, '\nTime Elapsed: %s' % (self.stopTime-self.startTime) + return result + + + def sortResult(self, result_list): + # unittest does not seems to run in any particular order. + # Here at least we want to group them together by class. + rmap = {} + classes = [] + for n,t,o,e in result_list: + cls = t.__class__ + if not rmap.has_key(cls): + rmap[cls] = [] + classes.append(cls) + rmap[cls].append((n,t,o,e)) + r = [(cls, rmap[cls]) for cls in classes] + return r + + + def getReportAttributes(self, result): + """ + Return report attributes as a list of (name, value). + Override this to add custom attributes. + """ + startTime = str(self.startTime)[:19] + duration = str(self.stopTime - self.startTime) + status = [] + if result.success_count: status.append('Pass %s' % result.success_count) + if result.failure_count: status.append('Failure %s' % result.failure_count) + if result.error_count: status.append('Error %s' % result.error_count ) + if status: + status = ' '.join(status) + else: + status = 'none' + return [ + ('Start Time', startTime), + ('Duration', duration), + ('Status', status), + ] + + + def generateReport(self, test, result): + report_attrs = self.getReportAttributes(result) + generator = 'HTMLTestRunner %s' % __version__ + stylesheet = self._generate_stylesheet() + heading = self._generate_heading(report_attrs) + report = self._generate_report(result) + ending = self._generate_ending() + output = self.HTML_TMPL % dict( + title = saxutils.escape(self.title), + generator = generator, + stylesheet = stylesheet, + heading = heading, + report = report, + ending = ending, + ) + self.stream.write(output.encode('utf8')) + + + def _generate_stylesheet(self): + return self.STYLESHEET_TMPL + + + def _generate_heading(self, report_attrs): + a_lines = [] + for name, value in report_attrs: + line = self.HEADING_ATTRIBUTE_TMPL % dict( + name = saxutils.escape(name), + value = saxutils.escape(value), + ) + a_lines.append(line) + heading = self.HEADING_TMPL % dict( + title = saxutils.escape(self.title), + parameters = ''.join(a_lines), + description = saxutils.escape(self.description), + ) + return heading + + + def _generate_report(self, result): + rows = [] + sortedResult = self.sortResult(result.result) + for cid, (cls, cls_results) in enumerate(sortedResult): + # subtotal for a class + np = nf = ne = 0 + for n,t,o,e in cls_results: + if n == 0: np += 1 + elif n == 1: nf += 1 + else: ne += 1 + + # format class description + if cls.__module__ == "__main__": + name = cls.__name__ + else: + name = "%s.%s" % (cls.__module__, cls.__name__) + doc = cls.__doc__ and cls.__doc__.split("\n")[0] or "" + desc = doc and '%s: %s' % (name, doc) or name + + row = self.REPORT_CLASS_TMPL % dict( + style = ne > 0 and 'errorClass' or nf > 0 and 'failClass' or 'passClass', + desc = desc, + count = np+nf+ne, + Pass = np, + fail = nf, + error = ne, + cid = 'c%s' % (cid+1), + ) + rows.append(row) + + for tid, (n,t,o,e) in enumerate(cls_results): + self._generate_report_test(rows, cid, tid, n, t, o, e) + + report = self.REPORT_TMPL % dict( + test_list = ''.join(rows), + count = str(result.success_count+result.failure_count+result.error_count), + Pass = str(result.success_count), + fail = str(result.failure_count), + error = str(result.error_count), + ) + return report + + + def _generate_report_test(self, rows, cid, tid, n, t, o, e): + # e.g. 'pt1.1', 'ft1.1', etc + has_output = bool(o or e) + tid = (n == 0 and 'p' or 'f') + 't%s.%s' % (cid+1,tid+1) + name = t.id().split('.')[-1] + doc = t.shortDescription() or "" + desc = doc and ('%s: %s' % (name, doc)) or name + tmpl = has_output and self.REPORT_TEST_WITH_OUTPUT_TMPL or self.REPORT_TEST_NO_OUTPUT_TMPL + + # o and e should be byte string because they are collected from stdout and stderr? + if isinstance(o,str): + # TODO: some problem with 'string_escape': it escape \n and mess up formating + # uo = unicode(o.encode('string_escape')) + uo = o.decode('latin-1') + else: + uo = o + if isinstance(e,str): + # TODO: some problem with 'string_escape': it escape \n and mess up formating + # ue = unicode(e.encode('string_escape')) + ue = e.decode('latin-1') + else: + ue = e + + script = self.REPORT_TEST_OUTPUT_TMPL % dict( + id = tid, + output = saxutils.escape(uo+ue), + ) + + row = tmpl % dict( + tid = tid, + Class = (n == 0 and 'hiddenRow' or 'none'), + style = n == 2 and 'errorCase' or (n == 1 and 'failCase' or 'none'), + desc = desc, + script = script, + status = self.STATUS[n], + ) + rows.append(row) + if not has_output: + return + + def _generate_ending(self): + return self.ENDING_TMPL + + +############################################################################## +# Facilities for running tests from the command line +############################################################################## + +# Note: Reuse unittest.TestProgram to launch test. In the future we may +# build our own launcher to support more specific command line +# parameters like test title, CSS, etc. +class TestProgram(unittest.TestProgram): + """ + A variation of the unittest.TestProgram. Please refer to the base + class for command line parameters. + """ + def runTests(self): + # Pick HTMLTestRunner as the default test runner. + # base class's testRunner parameter is not useful because it means + # we have to instantiate HTMLTestRunner before we know self.verbosity. + if self.testRunner is None: + self.testRunner = HTMLTestRunner(verbosity=self.verbosity) + unittest.TestProgram.runTests(self) + +main = TestProgram + +############################################################################## +# Executing this module from the command line +############################################################################## + +if __name__ == "__main__": + main(module=None) diff --git a/SimPEG/tests/api_TestResults.rst b/SimPEG/tests/api_TestResults.rst new file mode 100644 index 00000000..0d63a328 --- /dev/null +++ b/SimPEG/tests/api_TestResults.rst @@ -0,0 +1,355 @@ +.. _api_TestResults: + +.. raw:: html + + + + + +
+

Test Report

+

Start Time: 2013-11-05 15:24:44

+

Duration: 0:00:00.007500

+

Status: Pass 22

+ +

This demonstrates the report output by Prasanna.Yelsangikar.

+
+ + + +

Show + Summary + Failed + All +

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Test Group/Test caseCountPassFailErrorView
test_basemesh.TestBaseMesh111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
test_basemesh.TestMeshNumbers2D111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
Total222200 
+ diff --git a/SimPEG/tests/report.html b/SimPEG/tests/report.html new file mode 100644 index 00000000..be98c935 --- /dev/null +++ b/SimPEG/tests/report.html @@ -0,0 +1,2870 @@ + + + + + Test Results + + + + + + + + + +
+

Test Results

+

Start Time: 2013-11-05 15:38:27

+

Duration: 0:00:32.112639

+

Status: Pass 96 Error 3

+ +

SimPEG Test Report was automatically generated.

+
+ + + +

Show +Summary +Failed +All +

+ ++++++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Test Group/Test caseCountPassFailErrorView
test_basemesh.TestBaseMesh111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
test_basemesh.TestMeshNumbers2D111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
test_forward_DCproblem.DCProblemTests4103Detail
test_adjoint
 + + + + error + + + + +
test_dataObj
 + + + + error + + + + +
test_misfit
 + + + + error + + + + +
test_modelObj
 + + + + pass + + + + +
test_forward_problem.ProblemTests2200Detail
test_modelTransform
 + + + + pass + + + + +
test_regularization
 + + + + pass + + + + +
test_interpolation.TestInterpolation1D2200Detail
test_orderCC
 + + + + pass + + + + +
test_orderN
 + + + + pass + + + + +
test_interpolation.TestInterpolation2d6600Detail
test_orderCC
 + + + + pass + + + + +
test_orderEx
 + + + + pass + + + + +
test_orderEy
 + + + + pass + + + + +
test_orderFx
 + + + + pass + + + + +
test_orderFy
 + + + + pass + + + + +
test_orderN
 + + + + pass + + + + +
test_interpolation.TestInterpolation3D8800Detail
test_orderCC
 + + + + pass + + + + +
test_orderEx
 + + + + pass + + + + +
test_orderEy
 + + + + pass + + + + +
test_orderEz
 + + + + pass + + + + +
test_orderFx
 + + + + pass + + + + +
test_orderFy
 + + + + pass + + + + +
test_orderFz
 + + + + pass + + + + +
test_orderN
 + + + + pass + + + + +
test_LogicallyOrthogonalMesh.BasicLOMTests8800Detail
test_area_3D
pass
test_edge_2D
pass
test_edge_3D
pass
test_grid
pass
test_normals
pass
test_tangents
pass
test_vol_2D
pass
test_vol_3D
pass
test_massMatrices.TestInnerProducts: Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.6600Detail
test_order1_edges
 + + + + pass + + + + +
test_order1_faces
 + + + + pass + + + + +
test_order3_edges
 + + + + pass + + + + +
test_order3_faces
 + + + + pass + + + + +
test_order6_edges
 + + + + pass + + + + +
test_order6_faces
 + + + + pass + + + + +
test_massMatrices.TestInnerProducts2D: Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.6600Detail
test_order1_edges
 + + + + pass + + + + +
test_order1_faces
 + + + + pass + + + + +
test_order2_faces
 + + + + pass + + + + +
test_order3_edges
 + + + + pass + + + + +
test_order3_faces
 + + + + pass + + + + +
test_order6_edges
 + + + + pass + + + + +
test_operators.TestCurl1100Detail
test_order
 + + + + pass + + + + +
test_operators.TestFaceDiv1100Detail
test_order
 + + + + pass + + + + +
test_operators.TestFaceDiv2D1100Detail
test_order
 + + + + pass + + + + +
test_operators.TestNodalGrad1100Detail
test_order
 + + + + pass + + + + +
test_operators.TestNodalGrad2D1100Detail
test_order
 + + + + pass + + + + +
test_Solver.TestSolver101000Detail
test_directDiagonal_1
pass
test_directDiagonal_M
pass
test_directFactored_1
pass
test_directFactored_M
pass
test_directLower_1
pass
test_directLower_M
pass
test_directSpsolve_1
pass
test_directSpsolve_M
pass
test_directUpper_1
pass
test_directUpper_M
pass
test_tensorMesh.BasicTensorMeshTests7700Detail
test_area_3D
pass
test_edge_2D
pass
test_edge_3D
pass
test_vectorCC_2D
pass
test_vectorN_2D
pass
test_vol_2D
pass
test_vol_3D
pass
test_tensorMesh.TestPoissonEqn2200Detail
test_orderBackward
 + + + + pass + + + + +
test_orderForward
 + + + + pass + + + + +
test_utils.TestCheckDerivative3300Detail
test_simpleFail
 + + + + pass + + + + +
test_simpleFunction
 + + + + pass + + + + +
test_simplePass
 + + + + pass + + + + +
test_utils.TestSequenceFunctions8800Detail
test_indexCube_2D
pass
test_indexCube_3D
pass
test_invXXXBlockDiagonal
pass
test_mkvc1
pass
test_mkvc2
pass
test_mkvc3
pass
test_ndgrid_2D
pass
test_ndgrid_3D
pass
Total999603 
+ +
 
+ + + diff --git a/SimPEG/tests/runTests.py b/SimPEG/tests/runTests.py index b6662f00..d671324e 100644 --- a/SimPEG/tests/runTests.py +++ b/SimPEG/tests/runTests.py @@ -1,11 +1,47 @@ import glob import unittest +import HTMLTestRunner # This code will run all tests in directory named test_*.py +TITLE = 'Test Results' test_file_strings = glob.glob('test_*.py') module_strings = [str[0:len(str)-3] for str in test_file_strings] suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str in module_strings] testSuite = unittest.TestSuite(suites) -text_runner = unittest.TextTestRunner().run(testSuite) +unittest.TextTestRunner(verbosity=2).run(testSuite) + + +outfile = open("report.html", "w") +runner = HTMLTestRunner.HTMLTestRunner( + stream=outfile, + title=TITLE, + description='SimPEG Test Report was automatically generated.' + ) + +runner.run(testSuite) +outfile.close() + +reader = open("report.html", "r") +writer = open("../../docs/api_TestResults.rst", "w") + +writer.write('.. _api_TestResults:\n\nTest Results\n============\n\n.. raw:: html\n\n') + +go = False +for line in reader: + skip = False + if line == ' + + + +
+

Start Time: 2013-11-05 15:38:27

+

Duration: 0:00:32.112639

+

Status: Pass 96 Error 3

+ +

SimPEG Test Report was automatically generated.

+
+ + + +

Show + Summary + Failed + All +

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
Test Group/Test caseCountPassFailErrorView
test_basemesh.TestBaseMesh111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
test_basemesh.TestMeshNumbers2D111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
test_forward_DCproblem.DCProblemTests4103Detail
test_adjoint
 + + + + error + + + + +
test_dataObj
 + + + + error + + + + +
test_misfit
 + + + + error + + + + +
test_modelObj
 + + + + pass + + + + +
test_forward_problem.ProblemTests2200Detail
test_modelTransform
 + + + + pass + + + + +
test_regularization
 + + + + pass + + + + +
test_interpolation.TestInterpolation1D2200Detail
test_orderCC
 + + + + pass + + + + +
test_orderN
 + + + + pass + + + + +
test_interpolation.TestInterpolation2d6600Detail
test_orderCC
 + + + + pass + + + + +
test_orderEx
 + + + + pass + + + + +
test_orderEy
 + + + + pass + + + + +
test_orderFx
 + + + + pass + + + + +
test_orderFy
 + + + + pass + + + + +
test_orderN
 + + + + pass + + + + +
test_interpolation.TestInterpolation3D8800Detail
test_orderCC
 + + + + pass + + + + +
test_orderEx
 + + + + pass + + + + +
test_orderEy
 + + + + pass + + + + +
test_orderEz
 + + + + pass + + + + +
test_orderFx
 + + + + pass + + + + +
test_orderFy
 + + + + pass + + + + +
test_orderFz
 + + + + pass + + + + +
test_orderN
 + + + + pass + + + + +
test_LogicallyOrthogonalMesh.BasicLOMTests8800Detail
test_area_3D
pass
test_edge_2D
pass
test_edge_3D
pass
test_grid
pass
test_normals
pass
test_tangents
pass
test_vol_2D
pass
test_vol_3D
pass
test_massMatrices.TestInnerProducts: Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.6600Detail
test_order1_edges
 + + + + pass + + + + +
test_order1_faces
 + + + + pass + + + + +
test_order3_edges
 + + + + pass + + + + +
test_order3_faces
 + + + + pass + + + + +
test_order6_edges
 + + + + pass + + + + +
test_order6_faces
 + + + + pass + + + + +
test_massMatrices.TestInnerProducts2D: Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.6600Detail
test_order1_edges
 + + + + pass + + + + +
test_order1_faces
 + + + + pass + + + + +
test_order2_faces
 + + + + pass + + + + +
test_order3_edges
 + + + + pass + + + + +
test_order3_faces
 + + + + pass + + + + +
test_order6_edges
 + + + + pass + + + + +
test_operators.TestCurl1100Detail
test_order
 + + + + pass + + + + +
test_operators.TestFaceDiv1100Detail
test_order
 + + + + pass + + + + +
test_operators.TestFaceDiv2D1100Detail
test_order
 + + + + pass + + + + +
test_operators.TestNodalGrad1100Detail
test_order
 + + + + pass + + + + +
test_operators.TestNodalGrad2D1100Detail
test_order
 + + + + pass + + + + +
test_Solver.TestSolver101000Detail
test_directDiagonal_1
pass
test_directDiagonal_M
pass
test_directFactored_1
pass
test_directFactored_M
pass
test_directLower_1
pass
test_directLower_M
pass
test_directSpsolve_1
pass
test_directSpsolve_M
pass
test_directUpper_1
pass
test_directUpper_M
pass
test_tensorMesh.BasicTensorMeshTests7700Detail
test_area_3D
pass
test_edge_2D
pass
test_edge_3D
pass
test_vectorCC_2D
pass
test_vectorN_2D
pass
test_vol_2D
pass
test_vol_3D
pass
test_tensorMesh.TestPoissonEqn2200Detail
test_orderBackward
 + + + + pass + + + + +
test_orderForward
 + + + + pass + + + + +
test_utils.TestCheckDerivative3300Detail
test_simpleFail
 + + + + pass + + + + +
test_simpleFunction
 + + + + pass + + + + +
test_simplePass
 + + + + pass + + + + +
test_utils.TestSequenceFunctions8800Detail
test_indexCube_2D
pass
test_indexCube_3D
pass
test_invXXXBlockDiagonal
pass
test_mkvc1
pass
test_mkvc2
pass
test_mkvc3
pass
test_ndgrid_2D
pass
test_ndgrid_3D
pass
Total999603 
+ diff --git a/docs/index.rst b/docs/index.rst index 981e5ddc..12db6f97 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -47,6 +47,7 @@ Testing SimPEG :maxdepth: 2 api_Tests + api_TestResults Utility Codes From a241f14e2af556894331074d7b7d0f03da945e54 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 5 Nov 2013 16:11:49 -0800 Subject: [PATCH 45/45] updates to runTests.py --- SimPEG/tests/report.html | 2870 -------------------------------------- SimPEG/tests/runTests.py | 3 + docs/api_TestResults.rst | 611 ++++---- 3 files changed, 297 insertions(+), 3187 deletions(-) delete mode 100644 SimPEG/tests/report.html diff --git a/SimPEG/tests/report.html b/SimPEG/tests/report.html deleted file mode 100644 index be98c935..00000000 --- a/SimPEG/tests/report.html +++ /dev/null @@ -1,2870 +0,0 @@ - - - - - Test Results - - - - - - - - - -
-

Test Results

-

Start Time: 2013-11-05 15:38:27

-

Duration: 0:00:32.112639

-

Status: Pass 96 Error 3

- -

SimPEG Test Report was automatically generated.

-
- - - -

Show -Summary -Failed -All -

- -------- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Test Group/Test caseCountPassFailErrorView
test_basemesh.TestBaseMesh111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
test_basemesh.TestMeshNumbers2D111100Detail
test_meshDimensions
pass
test_mesh_nc
pass
test_mesh_nc_xyz
pass
test_mesh_ne
pass
test_mesh_nf
pass
test_mesh_numbers
pass
test_mesh_r_CC_M
pass
test_mesh_r_E_M
pass
test_mesh_r_E_V
pass
test_mesh_r_F_M
pass
test_mesh_r_F_V
pass
test_forward_DCproblem.DCProblemTests4103Detail
test_adjoint
 - - - - error - - - - -
test_dataObj
 - - - - error - - - - -
test_misfit
 - - - - error - - - - -
test_modelObj
 - - - - pass - - - - -
test_forward_problem.ProblemTests2200Detail
test_modelTransform
 - - - - pass - - - - -
test_regularization
 - - - - pass - - - - -
test_interpolation.TestInterpolation1D2200Detail
test_orderCC
 - - - - pass - - - - -
test_orderN
 - - - - pass - - - - -
test_interpolation.TestInterpolation2d6600Detail
test_orderCC
 - - - - pass - - - - -
test_orderEx
 - - - - pass - - - - -
test_orderEy
 - - - - pass - - - - -
test_orderFx
 - - - - pass - - - - -
test_orderFy
 - - - - pass - - - - -
test_orderN
 - - - - pass - - - - -
test_interpolation.TestInterpolation3D8800Detail
test_orderCC
 - - - - pass - - - - -
test_orderEx
 - - - - pass - - - - -
test_orderEy
 - - - - pass - - - - -
test_orderEz
 - - - - pass - - - - -
test_orderFx
 - - - - pass - - - - -
test_orderFy
 - - - - pass - - - - -
test_orderFz
 - - - - pass - - - - -
test_orderN
 - - - - pass - - - - -
test_LogicallyOrthogonalMesh.BasicLOMTests8800Detail
test_area_3D
pass
test_edge_2D
pass
test_edge_3D
pass
test_grid
pass
test_normals
pass
test_tangents
pass
test_vol_2D
pass
test_vol_3D
pass
test_massMatrices.TestInnerProducts: Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.6600Detail
test_order1_edges
 - - - - pass - - - - -
test_order1_faces
 - - - - pass - - - - -
test_order3_edges
 - - - - pass - - - - -
test_order3_faces
 - - - - pass - - - - -
test_order6_edges
 - - - - pass - - - - -
test_order6_faces
 - - - - pass - - - - -
test_massMatrices.TestInnerProducts2D: Integrate an function over a unit cube domain using edgeInnerProducts and faceInnerProducts.6600Detail
test_order1_edges
 - - - - pass - - - - -
test_order1_faces
 - - - - pass - - - - -
test_order2_faces
 - - - - pass - - - - -
test_order3_edges
 - - - - pass - - - - -
test_order3_faces
 - - - - pass - - - - -
test_order6_edges
 - - - - pass - - - - -
test_operators.TestCurl1100Detail
test_order
 - - - - pass - - - - -
test_operators.TestFaceDiv1100Detail
test_order
 - - - - pass - - - - -
test_operators.TestFaceDiv2D1100Detail
test_order
 - - - - pass - - - - -
test_operators.TestNodalGrad1100Detail
test_order
 - - - - pass - - - - -
test_operators.TestNodalGrad2D1100Detail
test_order
 - - - - pass - - - - -
test_Solver.TestSolver101000Detail
test_directDiagonal_1
pass
test_directDiagonal_M
pass
test_directFactored_1
pass
test_directFactored_M
pass
test_directLower_1
pass
test_directLower_M
pass
test_directSpsolve_1
pass
test_directSpsolve_M
pass
test_directUpper_1
pass
test_directUpper_M
pass
test_tensorMesh.BasicTensorMeshTests7700Detail
test_area_3D
pass
test_edge_2D
pass
test_edge_3D
pass
test_vectorCC_2D
pass
test_vectorN_2D
pass
test_vol_2D
pass
test_vol_3D
pass
test_tensorMesh.TestPoissonEqn2200Detail
test_orderBackward
 - - - - pass - - - - -
test_orderForward
 - - - - pass - - - - -
test_utils.TestCheckDerivative3300Detail
test_simpleFail
 - - - - pass - - - - -
test_simpleFunction
 - - - - pass - - - - -
test_simplePass
 - - - - pass - - - - -
test_utils.TestSequenceFunctions8800Detail
test_indexCube_2D
pass
test_indexCube_3D
pass
test_invXXXBlockDiagonal
pass
test_mkvc1
pass
test_mkvc2
pass
test_mkvc3
pass
test_ndgrid_2D
pass
test_ndgrid_3D
pass
Total999603 
- -
 
- - - diff --git a/SimPEG/tests/runTests.py b/SimPEG/tests/runTests.py index d671324e..c44f9a1e 100644 --- a/SimPEG/tests/runTests.py +++ b/SimPEG/tests/runTests.py @@ -1,3 +1,4 @@ +import os import glob import unittest import HTMLTestRunner @@ -45,3 +46,5 @@ for line in reader: writer.write(' '+line) writer.close() +reader.close() +os.remove("report.html") diff --git a/docs/api_TestResults.rst b/docs/api_TestResults.rst index 667ad91d..7b276006 100644 --- a/docs/api_TestResults.rst +++ b/docs/api_TestResults.rst @@ -185,9 +185,9 @@ Test Results -->
-

Start Time: 2013-11-05 15:38:27

-

Duration: 0:00:32.112639

-

Status: Pass 96 Error 3

+

Start Time: 2013-11-05 16:10:07

+

Duration: 0:00:31.280082

+

Status: Pass 99

SimPEG Test Report was automatically generated.

@@ -345,37 +345,48 @@ Test Results pass - + test_forward_DCproblem.DCProblemTests 4 - 1 + 4 + 0 0 - 3 Detail - -
test_adjoint
+ +
test_adjoint
+ pass + + + +
test_dataObj
- - error + + pass -