From a6da74c3476a07b3c8e9ac5f68dbcfa2ffe7c219 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 1 Oct 2013 20:33:57 -0700 Subject: [PATCH 01/11] Creating the inverse problem framework. Feedback welcome! --- SimPEG/GaussNewton.py | 64 +++++------ SimPEG/__init__.py | 1 + SimPEG/forward/Problem.py | 115 +++++++++++++++++++ SimPEG/forward/__init__.py | 1 + SimPEG/inverse/Optimize.py | 222 +++++++++++++++++++++++++++++++++++++ SimPEG/inverse/__init__.py | 1 + docs/api_Problem.rst | 8 ++ docs/index.rst | 3 +- 8 files changed, 382 insertions(+), 33 deletions(-) create mode 100644 SimPEG/forward/Problem.py create mode 100644 SimPEG/forward/__init__.py create mode 100644 SimPEG/inverse/Optimize.py create mode 100644 SimPEG/inverse/__init__.py create mode 100644 docs/api_Problem.rst diff --git a/SimPEG/GaussNewton.py b/SimPEG/GaussNewton.py index e0304a8d..bd6d37c3 100644 --- a/SimPEG/GaussNewton.py +++ b/SimPEG/GaussNewton.py @@ -2,25 +2,25 @@ import numpy as np import matplotlib.pyplot as plt from pylab import norm -def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3, tolX=1e-3, +def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3, tolX=1e-3, tolG=1e-3, eps=1e-16, xStop=[]): """ GaussNewton Optimization - + Input: ------ fctn - objective Function (lambda function) x0 - starting guess - + Output: ------- xOpt - numerical optimizer """ # initial output print "%s GaussNewton %s" % ('='*22,'='*22) - print "iter\tJc\t\tnorm(dJ)\tLS" + print "iter\tJc\t\tnorm(dJ)\tLS" print "%s" % '-'*57 - + # evaluate stopping criteria if xStop==[]: xStop=x0 @@ -31,26 +31,26 @@ def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3, xc = x0 STOP = np.zeros((5,1),dtype=bool) iterLS=0; iter=0 - + Jold = Jstop xOld=xc while 1: # evaluate objective function Jc,dJ,H = fctn(xc) print "%3d\t%1.2e\t%1.2e\t%d" % (iter, Jc[0],norm(dJ),iterLS) - + # check stopping rules STOP[0] = (iter>0) & (abs(Jc[0]-Jold[0]) <= tolJ*(1+abs(Jstop[0]))) STOP[1] = (iter>0) & (norm(xc-xOld) <= tolX*(1+norm(x0))) STOP[2] = norm(dJ) <= tolG*(1+abs(Jstop[0])) STOP[3] = norm(dJ) <= 1e3*eps STOP[4] = (iter >= maxIter) - if all(STOP[0:3]) | any(STOP[3:]): + if all(STOP[0:3]) | any(STOP[3:]): break - + # get search direction dx = np.linalg.solve(H,-dJ) - + # Armijo linesearch descent = np.dot(dJ.T,dx) LS =0; t = 1; iterLS=1 @@ -62,13 +62,13 @@ def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3, break iterLS = iterLS+1 t = .5*t - + # store old values Jold = Jc; xOld = xc - # update + # update xc = xt iter = iter +1 - + print "%s STOP! %s" % ('-'*25,'-'*25) print "%d : |Jc-Jold| = %1.4e <= tolJ*(1+|Jstop|) = %1.4e" % (STOP[0],abs(Jc[0]-Jold[0]),tolJ*(1+abs(Jstop[0]))) print "%d : |xc-xOld| = %1.4e <= tolX*(1+|x0|) = %1.4e" % (STOP[1],norm(xc-xOld),tolX*(1+norm(x0))) @@ -76,50 +76,50 @@ def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3, print "%d : |dJ| = %1.4e <= 1e3*eps = %1.4e" % (STOP[3],norm(dJ),1e3*eps) print "%d : iter = %3d\t <= maxIter\t = %3d" % (STOP[4],iter,maxIter) print "%s DONE! %s\n" % ('='*25,'='*25) - + return xc - + def Rosenbrock(x): """ Rosenbrock function for testing GaussNewton scheme """ - J = 100*(x[1]-x[0]**2)**2+(1-x[0])**2 + J = 100*(x[1]-x[0]**2)**2+(1-x[0])**2 dJ = 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]],dtype=float); - + return J,dJ,H - + def checkDerivative(fctn,x0): """ Basic derivative check - - Compares error decay of 0th and 1st order Taylor approximation at point + + Compares error decay of 0th and 1st order Taylor approximation at point x0 for a randomized search direction. - + Input: ------ fctn - function handle - x0 - point at which to check derivative + x0 - point at which to check derivative """ - + print "%s checkDerivative %s" % ('='*20,'='*20) print "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|" Jc,dJ,H = fctn(x0) - + dx = np.random.randn(len(x0),1) - + t = np.logspace(-1,-10,10) E0 = np.zeros(t.shape) E1 = np.zeros(t.shape) - + for i in range(0,10): Jt = fctn(x0+t[i]*dx) - E0[i] = norm(Jt[0]-Jc[0]) # 0th order Taylor - E1[i] = norm(Jt[0]-Jc[0]-t[i]*np.dot(dJ.T,dx)) # 1st order Taylor - + E0[i] = norm(Jt[0]-Jc[0]) # 0th order Taylor + E1[i] = norm(Jt[0]-Jc[0]-t[i]*np.dot(dJ.T,dx)) # 1st order Taylor + print "%d\t%1.2e\t%1.3e\t%1.3e" % (i,t[i],E0[i],E1[i]) - + print "%s DONE! %s\n" % ('='*25,'='*25) plt.figure() plt.clf() @@ -131,11 +131,11 @@ def checkDerivative(fctn,x0): plt.legend(['0th order', '1st order'],loc='upper left') plt.show() return - + if __name__ == '__main__': x0 = np.array([[2.6],[3.7]]) fctn = lambda x:Rosenbrock(x) checkDerivative(fctn,x0) xOpt = GaussNewton(fctn,x0,maxIter=20) print "xOpt=[%f,%f]" % (xOpt[0],xOpt[1]) - \ No newline at end of file + diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index ca9b9f7b..f36000a7 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1,3 +1,4 @@ from TensorMesh import TensorMesh from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh import utils +import inverse diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py new file mode 100644 index 00000000..f83ea505 --- /dev/null +++ b/SimPEG/forward/Problem.py @@ -0,0 +1,115 @@ +import numpy as np +from SimPEG.utils import mkvc, sdiag +norm = np.linalg.norm + + +class Problem(object): + """Problem is the base class for all geophysical forward problems in SimPEG""" + def __init__(self, mesh): + self.mesh = mesh + pass + + def residual(self, m): + pass + + 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))) + + def _test_modelTransformDeriv(self): + m = np.random.rand(5) + return checkDerivative(lambda m : [self.modelTransform(m), self.modelTransformDeriv(m)], m) + + def misfit(self, field): + """ + :param numpy.array field: geophysical field of interest + :rtype: float + :return: data misfit + + The data misfit using an l_2 norm is: + + .. math:: + + \mu_\\text{data} = {1\over 2}\left| \mathbf{W} (\mathbf{Pu} - 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.P*field - self.dobs) + return 0.5*mkvc(R).inner(mkvc(R)) + + def misfitDeriv(self, field): + """ + TODO: Change this documentation. + + :param numpy.array field: geophysical field of interest + :rtype: float + :return: data misfit derivative + + The data misfit using an l_2 norm is: + + .. math:: + + \mu_\\text{data} = {1\over 2}\left| \mathbf{W} (\mathbf{Pu} - 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.P*field - self.dobs) + # TODO: make in terms of the field and call Jt, e.g. if looping over RHSs using i: self.Jt(field[:,i],self.W[:,i]*R[:,i]) + return mkvc(R) + + def J(self, u): + pass + + def Jt(self, v): + pass + +if __name__ == '__main__': + from SimPEG.inverse import checkDerivative + + p = Problem(None) + m = np.random.rand(5) + checkDerivative(lambda m : [p.modelTransform(m), p.modelTransformDeriv(m)], m) diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py new file mode 100644 index 00000000..3c4d3728 --- /dev/null +++ b/SimPEG/forward/__init__.py @@ -0,0 +1 @@ +from Problem import * diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py new file mode 100644 index 00000000..53a11410 --- /dev/null +++ b/SimPEG/inverse/Optimize.py @@ -0,0 +1,222 @@ +import numpy as np +import matplotlib.pyplot as plt +from SimPEG.utils import mkvc, sdiag +norm = np.linalg.norm + + +class Minimize(object): + """docstring for Minimize""" + + name = "GeneralOptimizationAlgorithm" + + maxIter = 20 + maxIterLS = 10 + LSreduction = 1e-4 + LSshorten = 0.5 + tolF = 1e-4 + tolX = 1e-4 + tolG = 1e-4 + eps = 1e-16 + + def __init__(self, problem, **kwargs): + self.problem = problem + + # 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) + + def minimize(self, x0): + + self.startup(x0) + self.printInit() + + while True: + self.f, self.g, self.H = self.evalFunction(self.xc) + self.printIter() + if self.stoppingCriteria(): break + p = self.findSearchDirection() + xt, passLS = self.linesearch(p) + if not passLS: + xt = self.linesearchBreak(p) + self.doEndIteration(xt) + + self.printDone() + + return self.xc + + def startup(self, x0): + self._iter = 0 + self._iterLS = 0 + self._STOP = np.zeros((5,1),dtype=bool) + + self.x0 = x0 + self.xc = x0 + self.xOld = x0 + + def printInit(self): + 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) + + def printDone(self): + 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))) + 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) + 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 + + def stoppingCriteria(self): + if self._iter == 0: + self.fStop = self.f # Save this for stopping criteria + + # check stopping rules + self._STOP[0] = self._iter > 0 and (abs(self.f-self.fOld) <= self.tolF*(1+abs(self.fStop))) + self._STOP[1] = self._iter > 0 and (norm(self.xc-self.xOld) <= self.tolX*(1+norm(self.x0))) + self._STOP[2] = norm(self.g) <= self.tolG*(1+abs(self.fStop)) + self._STOP[3] = norm(self.g) <= 1e3*self.eps + self._STOP[4] = self._iter >= self.maxIter + return all(self._STOP[0:3]) | any(self._STOP[3:]) + + 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 + ft, temp, temp = self.evalFunction(xt, doDerivative=False) + if ft < self.f + t*self.LSreduction*descent: + break + iterLS += 1 + t = self.LSshorten*t + + self._iterLS = iterLS + return xt, iterLS < self.maxIterLS + + def linesearchBreak(self, p): + raise Exception('The linesearch got broken. Boo.') + + def doEndIteration(self, xt): + # store old values + self.fOld = self.f + self.xOld, self.xc = self.xc, xt + self._iter += 1 + + +class GaussNewton(Minimize): + name = 'GaussNewton' + def findSearchDirection(self): + return np.linalg.solve(self.H,-self.g) + + +class SteepestDescent(Minimize): + name = 'SteepestDescent' + def findSearchDirection(self): + return -self.g + + + +def Rosenbrock(x): + """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 + + +def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): + """ + Basic derivative check + + Compares error decay of 0th and 1st order Taylor approximation at point + x0 for a randomized search direction. + + Input: + ------ + fctn - function handle + x0 - point at which to check derivative + """ + + print "%s checkDerivative %s" % ('='*20, '='*20) + print "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|\tOrder\n%s" % ('-'*57) + + Jc = fctn(x0) + + x0 = mkvc(x0) + + if dx is None: + dx = np.random.randn(len(x0)) + + t = np.logspace(-1, -num, num) + E0 = np.ones(t.shape) + E1 = np.ones(t.shape) + + l2norm = lambda x: np.sqrt(np.inner(x, x)) # because np.norm breaks if they are scalars? + 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 + 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]) + + tolerance = 0.9 + expectedOrder = 2 + eps = 1e-10 + order0 = order0[E0[1:] > eps] + order1 = order1[E1[1:] > eps] + belowTol = order1.size == 0 and order0.size > 0 + correctOrder = order1.size > 0 and np.mean(order1) > tolerance * expectedOrder + + passTest = belowTol or correctOrder + + if passTest: + print "%s PASS! %s\n" % ('='*25, '='*25) + else: + print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57) + + if plotIt: + plt.figure() + plt.clf() + plt.loglog(t, E0, 'b') + plt.loglog(t, E1, 'g--') + plt.title('checkDerivative') + plt.xlabel('h') + plt.ylabel('error of Taylor approximation') + plt.legend(['0th order', '1st order'], loc='upper left') + plt.show() + + return passTest + +if __name__ == '__main__': + x0 = np.array([2.6, 3.7]) + checkDerivative(Rosenbrock, x0, plotIt=False) + xOpt = GaussNewton(Rosenbrock, maxIter=20).minimize(x0) + print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1]) + xOpt = SteepestDescent(Rosenbrock, maxIter=20, maxIterLS=15).minimize(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/inverse/__init__.py b/SimPEG/inverse/__init__.py new file mode 100644 index 00000000..b2a5e506 --- /dev/null +++ b/SimPEG/inverse/__init__.py @@ -0,0 +1 @@ +from Optimize import * diff --git a/docs/api_Problem.rst b/docs/api_Problem.rst new file mode 100644 index 00000000..d43616d9 --- /dev/null +++ b/docs/api_Problem.rst @@ -0,0 +1,8 @@ +.. _api_Problem: + +Problem +******* + +.. automodule:: SimPEG.forward.Problem + :members: + :undoc-members: diff --git a/docs/index.rst b/docs/index.rst index 22c9d5e9..48ebe940 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -38,12 +38,13 @@ Inversion api_GaussNewton -Example Problems +Forward Problems ================ .. toctree:: :maxdepth: 2 + api_Problem Project Index & Search From 7a4ccc9a383e0fd9373d53855d511f031b63f57d Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 1 Oct 2013 23:58:35 -0700 Subject: [PATCH 02/11] Updates to problem. Problem is defined as a general PDE that has a field and a model. --- SimPEG/forward/Problem.py | 173 +++++++++++++++++++++++++++++++++----- 1 file changed, 153 insertions(+), 20 deletions(-) diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index f83ea505..48848a07 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -4,14 +4,132 @@ norm = np.linalg.norm class Problem(object): - """Problem is the base class for all geophysical forward problems in SimPEG""" + """ + Problem is the base class for all geophysical forward problems in SimPEG. + + + The problem is a partial differential equation of the form: + + .. math:: + c(m, u) = 0 + + Here, m is the model and u is the field (or fields). + Given the model, m, we can calculate the fields u(m), + however, the data we collect is a subset of the fields, + and can be defined by a linear projection, P. + + .. math:: + d_\\text{pred} = Pu(m) + + We are interested in how changing the model transforms the data, + as such we can take write the Taylor expansion: + + .. math:: + Pu(m + hv) = Pu(m) + hP\\frac{\partial u(m)}{\partial m} v + \mathcal{O}(h^2 \left\| v \\right\| ) + + We can linearize and define the sensitivity matrix as: + + .. math:: + J = P\\frac{\partial u}{\partial m} + + The sensitivity matrix, and it's transpose will be used in the inverse problem + to (locally) find how model parameters change the data, and optimize! + """ + def __init__(self, mesh): self.mesh = mesh + + @property + def RHS(self): + """ + Source matrix. + """ + return self._RHS + @RHS.setter + 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): + """ + Projection matrix. + + .. math:: + d_\\text{pred} = Pu(m) + """ + return self._P + @P.setter + def P(self, value): + self._P = value + + + @property + def dobs(self): + """ + Observed data. + """ + return self._dobs + @dobs.setter + def dobs(self, value): + self._P = value + + + def J(self, u): + """ + Working with the general PDE, c(m, u) = 0, where m is the model and u is the field, + the sensitivity is defined as: + + .. math:: + J = P\\frac{\partial u}{\partial m} + + We can take the derivative of the PDE: + + .. math:: + \\nabla_m c(m, u) \delta m + \\nabla_u c(m, u) \delta u = 0 + + If the forward problem is invertible, then we can rearrange for du/dm: + + .. math:: + J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u) + + This can often be computed given a vector (i.e. J(v)) rather than stored, as J is a large dense matrix. + + """ pass - def residual(self, m): + def Jt(self, v): + """ + Transpose of J + """ pass + def field(self, m): + """ + The fields. + """ + 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 @@ -61,9 +179,10 @@ class Problem(object): m = np.random.rand(5) return checkDerivative(lambda m : [self.modelTransform(m), self.modelTransformDeriv(m)], m) - def misfit(self, field): + def misfit(self, m, R=None): """ - :param numpy.array field: geophysical field of interest + :param numpy.array m: geophysical model + :param numpy.array R: residual, R = W o (dpred - dobs) :rtype: float :return: data misfit @@ -71,41 +190,55 @@ class Problem(object): .. math:: - \mu_\\text{data} = {1\over 2}\left| \mathbf{W} (\mathbf{Pu} - d_\\text{obs}) \\right|_2^2 + \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.P*field - self.dobs) - return 0.5*mkvc(R).inner(mkvc(R)) + if R is None: + R = self.W*(self.dpred(m) - self.dobs) - def misfitDeriv(self, field): + R = mkvc(R) + return 0.5*R.inner(R) + + def misfitDeriv(self, m, R=None, u=None): """ - TODO: Change this documentation. - - :param numpy.array field: geophysical field of interest - :rtype: float + :param numpy.array m: geophysical model + :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} (\mathbf{Pu} - d_\\text{obs}) \\right|_2^2 + \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 + + \mathbf{R} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs} + + \mu_\\text{data} = {1\over 2}\left| \mathbf{W \circ R} \\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. + + 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.P*field - self.dobs) - # TODO: make in terms of the field and call Jt, e.g. if looping over RHSs using i: self.Jt(field[:,i],self.W[:,i]*R[:,i]) - return mkvc(R) + if R is None: + R = self.W*(self.dpred(m, u=u) - self.dobs) - def J(self, u): - pass + dmisfit = 0 + for i in range(self.RHS.shape[1]): # Loop over each right hand side + dmisfit += self.Jt(u[:,i], self.W[:,i]*R[:,i]) + + return dmisfit - def Jt(self, v): - pass if __name__ == '__main__': from SimPEG.inverse import checkDerivative From d8c676015e0b971ba825ec79865b90b9a45fb2fa Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 2 Oct 2013 15:24:51 -0700 Subject: [PATCH 03/11] Fixes to ModelBuilder. Start of the DCProblem. --- SimPEG/DiffOperators.py | 12 ++-- SimPEG/InnerProducts.py | 2 +- SimPEG/forward/DCProblem/DCProblem.py | 93 +++++++++++++++++++++++++++ SimPEG/forward/DCProblem/DCutils.py | 29 +++++++++ SimPEG/forward/Problem.py | 76 ++++++++++++++++++---- SimPEG/utils/ModelBuilder.py | 48 ++++++++------ 6 files changed, 224 insertions(+), 36 deletions(-) create mode 100644 SimPEG/forward/DCProblem/DCProblem.py create mode 100644 SimPEG/forward/DCProblem/DCutils.py diff --git a/SimPEG/DiffOperators.py b/SimPEG/DiffOperators.py index 9a8a9db9..d50cc5ba 100644 --- a/SimPEG/DiffOperators.py +++ b/SimPEG/DiffOperators.py @@ -1,6 +1,6 @@ import numpy as np from scipy import sparse as sp -from utils import mkvc, sdiag, speye, kron3, spzeros +from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros def ddx(n): @@ -287,15 +287,19 @@ class DiffOperators(object): nodalVectorAve = property(**nodalVectorAve()) def getEdgeMass(self, materialProp=None): - """mass matix for products of edge functions w'*M(materialProp)*e""" + """mass matrix for products of edge functions w'*M(materialProp)*e""" if(materialProp is None): materialProp = np.ones(self.nC) Av = self.edgeAve return sdiag(Av.T * (self.vol * mkvc(materialProp))) def getFaceMass(self, materialProp=None): - """mass matix for products of edge functions w'*M(materialProp)*e""" + """mass matrix for products of face functions w'*M(materialProp)*f""" if(materialProp is None): materialProp = np.ones(self.nC) Av = self.faceAve - return sdiag(Av.T*(self.vol*mkvc(materialProp))) + return sdiag(Av.T * (self.vol * mkvc(materialProp))) + + def getFaceMassDeriv(self): + Av = self.faceAve + return Av.T * sdiag(self.vol) diff --git a/SimPEG/InnerProducts.py b/SimPEG/InnerProducts.py index fca632f4..9fe84ac3 100644 --- a/SimPEG/InnerProducts.py +++ b/SimPEG/InnerProducts.py @@ -1,5 +1,5 @@ from scipy import sparse as sp -from utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal +from SimPEG.utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal import numpy as np diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py new file mode 100644 index 00000000..b94b35b5 --- /dev/null +++ b/SimPEG/forward/DCProblem/DCProblem.py @@ -0,0 +1,93 @@ +from SimPEG import TensorMesh +from SimPEG.forward import Problem, SyntheticProblem +from SimPEG.utils import ModelBuilder +import numpy as np +import scipy.sparse.linalg as linalg +import DCutils + +class DCProblem(Problem): + """docstring for DCProblem""" + def __init__(self, mesh): + super(DCProblem, self).__init__(mesh) + self.mesh.setCellGradBC('neumann') + + def createMatrix(self, m): + D = self.mesh.faceDiv + G = self.mesh.cellGrad + sigma = self.modelTransform(m) + Msig = self.mesh.getFaceMass(sigma) + A = D*Msig*G + return A.tocsc() + + 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]) + + return phi + + def J(self, m, v, u=None, RHSii=0, solve=None): + P = self.P + D = self.mesh.faceDiv + G = self.mesh.cellGrad + A = self.createMatrix(m) + Av_dm = self.mesh.getFaceMassDeriv() + mT_dm = self.modelTransform(m) + + dCdu = A + dCdm = - D * ( sdiag( G * u[:, RHSii] ) * ( Av_dm * ( mT_dm * v ) ) ) + + if solve is None: + solve = linalg.factorized(dCdu) + + return - P * solve(dCdm) + + + +if __name__ == '__main__': + # Create the mesh + h1 = np.ones(100) + h2 = np.ones(100) + mesh = TensorMesh([h1,h2]) + + # Create some parameters for the model + sig1 = 1 + sig2 = 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) + + + # Set up the projection + nelec = 50 + spacelec = 2 + surfloc = 0.5 + elecini = 0.5 + 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) + + + # Create some data + class syntheticDCProblem(DCProblem, SyntheticProblem): + pass + + synthetic = syntheticDCProblem(mesh); + synthetic.P = Q.T + synthetic.RHS = q + dobs, Wd = synthetic.createData(mSynth) + + # Now set up the problem to do some minimization + problem = DCProblem(mesh) + + + diff --git a/SimPEG/forward/DCProblem/DCutils.py b/SimPEG/forward/DCProblem/DCutils.py new file mode 100644 index 00000000..f3445096 --- /dev/null +++ b/SimPEG/forward/DCProblem/DCutils.py @@ -0,0 +1,29 @@ +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/Problem.py b/SimPEG/forward/Problem.py index 48848a07..ff8bdb00 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -84,8 +84,15 @@ class Problem(object): self._P = value - def J(self, u): + def J(self, m, v, u=None, RHSii=0): """ + :param numpy.array m: model + :param numpy.array v: vector to multiply + :param numpy.array u: fields + :param int RHSii: which RHS to calculate sensitivity too + :rtype: numpy.array + :return: Jv + Working with the general PDE, c(m, u) = 0, where m is the model and u is the field, the sensitivity is defined as: @@ -107,15 +114,26 @@ class Problem(object): """ pass - def Jt(self, v): + def Jt(self, m, v, u=None, RHSii=0): """ + :param numpy.array m: model + :param numpy.array v: vector to multiply + :param numpy.array u: fields + :param int RHSii: which RHS to calculate sensitivity too + :rtype: numpy.array + :return: JTv + Transpose of J """ pass def field(self, m): """ - The fields. + The field given the model. + + .. math:: + u(m) + """ pass @@ -179,10 +197,10 @@ class Problem(object): m = np.random.rand(5) return checkDerivative(lambda m : [self.modelTransform(m), self.modelTransformDeriv(m)], m) - def misfit(self, m, R=None): + def misfit(self, m, u=None): """ :param numpy.array m: geophysical model - :param numpy.array R: residual, R = W o (dpred - dobs) + :param numpy.array u: fields :rtype: float :return: data misfit @@ -195,15 +213,15 @@ class Problem(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. """ - if R is None: - R = self.W*(self.dpred(m) - self.dobs) + R = self.W*(self.dpred(m, u=u) - self.dobs) R = mkvc(R) return 0.5*R.inner(R) - def misfitDeriv(self, m, R=None, u=None): + def misfitDeriv(self, m, u=None): """ :param numpy.array m: geophysical model + :param numpy.array u: fields :rtype: numpy.array :return: data misfit derivative @@ -213,6 +231,12 @@ class Problem(object): \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{d}_\\text{pred} - \mathbf{d}_\\text{obs} \mu_\\text{data} = {1\over 2}\left| \mathbf{W \circ R} \\right|_2^2 @@ -230,8 +254,7 @@ class Problem(object): if u is None: u = self.field(m) - if R is None: - R = self.W*(self.dpred(m, u=u) - self.dobs) + 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 @@ -240,9 +263,40 @@ class Problem(object): return dmisfit +class SyntheticProblem(object): + """ + Has helpful functions when dealing with synthetic problems + + To use this class, inherit to your problem:: + + class mySyntheticExample(Problem, SyntheticProblem): + pass + """ + def createData(self, m, std=0.05): + """ + :param numpy.array m: geophysical model + :param numpy.array std: standard deviation + :rtype: numpy.array, numpy.array + :return: dobs, Wd + + Create synthetic data given a model, and a standard deviation. + + Returns the observed data with random Gaussian noise + and Wd which is the same size as data, and can be used to weight the inversion. + """ + dobs = self.dpred(m) + dobs = dobs + noise = std*abs(dobs)*np.random.randn(*dobs.shape) + dobs = dobs+noise + eps = np.linalg.norm(mkvc(dobs),2)*1e-5 + Wd = 1/(abs(dobs)*std+eps) + return dobs, Wd + + + if __name__ == '__main__': from SimPEG.inverse import checkDerivative p = Problem(None) m = np.random.rand(5) - checkDerivative(lambda m : [p.modelTransform(m), p.modelTransformDeriv(m)], m) + checkDerivative(lambda m : [p.modelTransform(m), p.modelTransformDeriv(m)], m, plotIt=False) diff --git a/SimPEG/utils/ModelBuilder.py b/SimPEG/utils/ModelBuilder.py index 3b1f977f..527d5eef 100644 --- a/SimPEG/utils/ModelBuilder.py +++ b/SimPEG/utils/ModelBuilder.py @@ -15,10 +15,6 @@ def getIndecesBlock(p0,p1,ccMesh): The points p0 and p1 must live in the the same dimensional space as the mesh. """ - # Validation of the input - assert type(p0) == np.ndarray, "Vector must be a numpy array" - assert type(p1) == np.ndarray, "Vector must be a numpy array" - # Validation: p0 and p1 live in the same dimensional space assert len(p0) == len(p1), "Dimension mismatch. len(p0) != len(p1)" @@ -47,7 +43,7 @@ def getIndecesBlock(p0,p1,ccMesh): ind = np.where(indX & indY) - else: + elif dimMesh == 3: # Define the points x1 = p0[0] y1 = p0[1] @@ -98,13 +94,16 @@ def defineTwoLayeredConductivity(depth,ccMesh,condVals): # Identify 1st cell centered reference point p0[0] = ccMesh[0,0] - p0[1] = ccMesh[0,1] - p0[2] = ccMesh[0,2] + if dim>1: p0[1] = ccMesh[0,1] + if dim>2: p0[2] = ccMesh[0,2] # Identify the last cell-centered reference point p1[0] = ccMesh[-1,0] - p1[1] = ccMesh[-1,1] - p1[2] = ccMesh[-1,2] - depth; + if dim>1: p1[1] = ccMesh[-1,1] + if dim>2: p1[2] = ccMesh[-1,2] + + # The depth is always defined on the last one. + p1[len(p1)-1] -= depth ind = getIndecesBlock(p0,p1,ccMesh) @@ -117,23 +116,24 @@ def scalarConductivity(ccMesh,pFunction): Define the distribution conductivity in the mesh according to the analytical expression given in pFunction """ - xCC = ccMesh[:,0] - yCC = ccMesh[:,1] - zCC = ccMesh[:,2] + dim = np.size(ccMesh[0,:]) + CC = [ccMesh[:,0]] + if dim>1: CC.append(ccMesh[:,1]) + if dim>2: CC.append(ccMesh[:,2]) - sigma = pFunction(xCC,yCC,zCC) + + sigma = pFunction(*CC) return sigma if __name__ == '__main__': - import sys - sys.path.append('../') - from TensorMesh import TensorMesh + from SimPEG import TensorMesh + from matplotlib import pyplot as plt # Define the mesh - testDim = 3 + testDim = 2 h1 = 0.3*np.ones(7) h1[0] = 0.5 h1[-1] = 0.6 @@ -157,8 +157,8 @@ if __name__ == '__main__': # ------------------- Test conductivities! -------------------------- print('Testing 1 block conductivity') - p0 = np.array([0.5,0.5,0.5]) - p1 = np.array([1.0,1.0,1.0]) + p0 = np.array([0.5,0.5,0.5])[:testDim] + p1 = np.array([1.0,1.0,1.0])[:testDim] condVals = np.array([100,1e-6]) sigma = defineBlockConductivity(p0,p1,ccMesh,condVals) @@ -167,6 +167,7 @@ if __name__ == '__main__': print sigma.shape M.plotImage(sigma) print 'Done with block! :)' + plt.show() # ----------------------------------------- print('Testing the two layered model') @@ -178,11 +179,17 @@ if __name__ == '__main__': M.plotImage(sigma) print sigma print 'layer model!' + plt.show() # ----------------------------------------- print('Testing scalar conductivity') - pFunction = lambda x,y,z: np.exp(x+y+z) + if testDim == 1: + pFunction = lambda x: np.exp(x) + elif testDim == 2: + pFunction = lambda x,y: np.exp(x+y) + elif testDim == 3: + pFunction = lambda x,y,z: np.exp(x+y+z) sigma = scalarConductivity(ccMesh,pFunction) @@ -190,5 +197,6 @@ if __name__ == '__main__': M.plotImage(sigma) print sigma print 'Scalar conductivity defined!' + plt.show() # ----------------------------------------- From 21501e7482773ecd6e5f5fbe9a39e730cdd68da0 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 3 Oct 2013 10:37:43 -0700 Subject: [PATCH 04/11] DC Problem tested and working. --- SimPEG/forward/DCProblem/DCProblem.py | 50 +++++++++++++++++---- SimPEG/forward/Problem.py | 62 ++++++++++++++++++++++----- 2 files changed, 93 insertions(+), 19 deletions(-) diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py index b94b35b5..fa557ea5 100644 --- a/SimPEG/forward/DCProblem/DCProblem.py +++ b/SimPEG/forward/DCProblem/DCProblem.py @@ -1,6 +1,7 @@ from SimPEG import TensorMesh from SimPEG.forward import Problem, SyntheticProblem -from SimPEG.utils import ModelBuilder +from SimPEG.inverse import checkDerivative +from SimPEG.utils import ModelBuilder, sdiag import numpy as np import scipy.sparse.linalg as linalg import DCutils @@ -30,22 +31,39 @@ class DCProblem(Problem): return phi - def J(self, m, v, u=None, RHSii=0, solve=None): + def J(self, m, v, u=None, solve=None): P = self.P D = self.mesh.faceDiv G = self.mesh.cellGrad A = self.createMatrix(m) Av_dm = self.mesh.getFaceMassDeriv() - mT_dm = self.modelTransform(m) + mT_dm = self.modelTransformDeriv(m) dCdu = A - dCdm = - D * ( sdiag( G * u[:, RHSii] ) * ( Av_dm * ( mT_dm * v ) ) ) + dCdm = D * ( sdiag( G * u ) * ( Av_dm * ( mT_dm * v ) ) ) if solve is None: solve = linalg.factorized(dCdu) - return - P * solve(dCdm) + Jv = - P * solve(dCdm) + return Jv + def Jt(self, m, v, u=None, solve=None): + P = self.P + D = self.mesh.faceDiv + G = self.mesh.cellGrad + A = self.createMatrix(m) + Av_dm = self.mesh.getFaceMassDeriv() + mT_dm = self.modelTransformDeriv(m) + + dCdu = A.T + + if solve is None: + solve = linalg.factorized(dCdu.tocsc()) + w = solve(P.T*v) + + Jtv = - mT_dm.T * ( Av_dm.T * ( sdiag( G * u ) * ( D.T * w ) ) ) + return Jtv if __name__ == '__main__': @@ -75,19 +93,35 @@ if __name__ == '__main__': 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) - + P = Q.T # Create some data class syntheticDCProblem(DCProblem, SyntheticProblem): pass synthetic = syntheticDCProblem(mesh); - synthetic.P = Q.T + synthetic.P = P synthetic.RHS = q - dobs, Wd = synthetic.createData(mSynth) + dobs, Wd = synthetic.createData(mSynth, std=0.05) # Now set up the problem to do some minimization problem = DCProblem(mesh) + problem.P = P + problem.RHS = q + problem.W = Wd + problem.dobs = dobs + m0 = mesh.gridCC[:,0]*0+sig1 + print problem.misfit(m0) + print problem.misfit(mSynth) + # Check Derivative + derChk = lambda m: [problem.misfit(m), problem.misfitDeriv(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)) diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index ff8bdb00..39ddf8a8 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -81,15 +81,14 @@ class Problem(object): return self._dobs @dobs.setter def dobs(self, value): - self._P = value + self._dobs = value - def J(self, m, v, u=None, RHSii=0): + def J(self, m, v, u=None): """ :param numpy.array m: model :param numpy.array v: vector to multiply :param numpy.array u: fields - :param int RHSii: which RHS to calculate sensitivity too :rtype: numpy.array :return: Jv @@ -114,12 +113,11 @@ class Problem(object): """ pass - def Jt(self, m, v, u=None, RHSii=0): + def Jt(self, m, v, u=None): """ :param numpy.array m: model :param numpy.array v: vector to multiply :param numpy.array u: fields - :param int RHSii: which RHS to calculate sensitivity too :rtype: numpy.array :return: JTv @@ -216,7 +214,7 @@ class Problem(object): R = self.W*(self.dpred(m, u=u) - self.dobs) R = mkvc(R) - return 0.5*R.inner(R) + return 0.5*R.dot(R) def misfitDeriv(self, m, u=None): """ @@ -237,9 +235,7 @@ class Problem(object): \mathbf{d}_\\text{pred} = \mathbf{Pu(m)} - \mathbf{R} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs} - - \mu_\\text{data} = {1\over 2}\left| \mathbf{W \circ R} \\right|_2^2 + \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. @@ -248,7 +244,7 @@ class Problem(object): .. math:: - \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top (\mathbf{W \circ R}) + \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R} """ if u is None: @@ -258,7 +254,51 @@ class Problem(object): dmisfit = 0 for i in range(self.RHS.shape[1]): # Loop over each right hand side - dmisfit += self.Jt(u[:,i], self.W[:,i]*R[:,i]) + 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 From 56aa003a3ceaffaa65ed50635869ec7223c15c4d Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 3 Oct 2013 12:16:10 -0700 Subject: [PATCH 05/11] Start of DC documentation. --- SimPEG/forward/DCProblem/DCProblem.py | 14 ++++++++++++++ SimPEG/forward/DCProblem/__init__.py | 2 ++ SimPEG/forward/__init__.py | 1 + docs/api_Problem.rst | 18 ++++++++++++++++++ 4 files changed, 35 insertions(+) create mode 100644 SimPEG/forward/DCProblem/__init__.py diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py index fa557ea5..cf33d86c 100644 --- a/SimPEG/forward/DCProblem/DCProblem.py +++ b/SimPEG/forward/DCProblem/DCProblem.py @@ -32,6 +32,20 @@ class DCProblem(Problem): return phi def J(self, m, v, u=None, solve=None): + """ + :param numpy.array m: model + :param numpy.array v: vector to multiply + :param numpy.array u: fields + :rtype: numpy.array + :return: Jv + + .. math:: + c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0 + + \\nabla_u (A(m)u - q) = A(m) + + \\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m (M(mT(m))) + """ P = self.P D = self.mesh.faceDiv G = self.mesh.cellGrad diff --git a/SimPEG/forward/DCProblem/__init__.py b/SimPEG/forward/DCProblem/__init__.py new file mode 100644 index 00000000..a868cf80 --- /dev/null +++ b/SimPEG/forward/DCProblem/__init__.py @@ -0,0 +1,2 @@ +from DCProblem import * +from DCutils import * diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py index 3c4d3728..fe849d41 100644 --- a/SimPEG/forward/__init__.py +++ b/SimPEG/forward/__init__.py @@ -1 +1,2 @@ from Problem import * +import DCProblem diff --git a/docs/api_Problem.rst b/docs/api_Problem.rst index d43616d9..83250f3e 100644 --- a/docs/api_Problem.rst +++ b/docs/api_Problem.rst @@ -1,8 +1,26 @@ .. _api_Problem: + + Problem ******* .. automodule:: SimPEG.forward.Problem :members: :undoc-members: + + +DCProblem +********* + +.. automodule:: SimPEG.forward.DCProblem.DCProblem + :members: + :undoc-members: + + +DCutils +******* + +.. automodule:: SimPEG.forward.DCProblem.DCutils + :members: + :undoc-members: From bed86b9ce4bdabc02d9a2ae90e14ff450a017e06 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 3 Oct 2013 12:28:08 -0700 Subject: [PATCH 06/11] Update documentation and show example field. --- SimPEG/forward/DCProblem/DCProblem.py | 31 +++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py index cf33d86c..491123c3 100644 --- a/SimPEG/forward/DCProblem/DCProblem.py +++ b/SimPEG/forward/DCProblem/DCProblem.py @@ -7,12 +7,29 @@ import scipy.sparse.linalg as linalg import DCutils class DCProblem(Problem): - """docstring for DCProblem""" + """ + **DCProblem** + + Geophysical DC resistivity problem. + + """ def __init__(self, mesh): super(DCProblem, self).__init__(mesh) self.mesh.setCellGradBC('neumann') def createMatrix(self, m): + """ + Makes the matrix A(m) for the DC resistivity problem. + + :param numpy.array m: model + :rtype: scipy.csc_matrix + :return: A(m) + + .. math:: + c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0 + + Where M() is the mass matrix and mT is the model transform. + """ D = self.mesh.faceDiv G = self.mesh.cellGrad sigma = self.modelTransform(m) @@ -44,7 +61,14 @@ class DCProblem(Problem): \\nabla_u (A(m)u - q) = A(m) - \\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m (M(mT(m))) + \\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) + + Where M() is the mass matrix and mT is the model transform. + + .. math:: + J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u) + + J(v) = - P ( A(m)^{-1} ( G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) v ) ) """ P = self.P D = self.mesh.faceDiv @@ -118,6 +142,9 @@ if __name__ == '__main__': synthetic.RHS = q dobs, Wd = synthetic.createData(mSynth, std=0.05) + u = synthetic.field(mSynth) + mesh.plotImage(u[:,10], showIt=True) + # Now set up the problem to do some minimization problem = DCProblem(mesh) problem.P = P From 7cf8ef231003a4fb2313b484740bb8e9702e0257 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 3 Oct 2013 14:55:22 -0700 Subject: [PATCH 07/11] Moved testing code to the test directory. added __init__.py and some documentation. --- SimPEG/GaussNewton.py | 141 -------------------- SimPEG/forward/DCProblem/DCProblem.py | 2 +- SimPEG/forward/Problem.py | 4 - SimPEG/inverse/Optimize.py | 75 ----------- SimPEG/tests/{OrderTest.py => TestUtils.py} | 132 ++++++++++++++---- SimPEG/tests/__init__.py | 2 + SimPEG/tests/test_massMatrices.py | 2 +- SimPEG/tests/test_operators.py | 2 +- SimPEG/tests/test_tensorMesh.py | 6 +- SimPEG/tests/test_utils.py | 4 +- docs/api_GaussNewton.rst | 8 -- docs/api_Optimize.rst | 8 ++ docs/api_Tests.rst | 8 ++ docs/index.rst | 24 ++-- 14 files changed, 144 insertions(+), 274 deletions(-) delete mode 100644 SimPEG/GaussNewton.py rename SimPEG/tests/{OrderTest.py => TestUtils.py} (56%) create mode 100644 SimPEG/tests/__init__.py delete mode 100644 docs/api_GaussNewton.rst create mode 100644 docs/api_Optimize.rst create mode 100644 docs/api_Tests.rst diff --git a/SimPEG/GaussNewton.py b/SimPEG/GaussNewton.py deleted file mode 100644 index bd6d37c3..00000000 --- a/SimPEG/GaussNewton.py +++ /dev/null @@ -1,141 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from pylab import norm - -def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3, tolX=1e-3, - tolG=1e-3, eps=1e-16, xStop=[]): - """ - GaussNewton Optimization - - Input: - ------ - fctn - objective Function (lambda function) - x0 - starting guess - - Output: - ------- - xOpt - numerical optimizer - """ - # initial output - print "%s GaussNewton %s" % ('='*22,'='*22) - print "iter\tJc\t\tnorm(dJ)\tLS" - print "%s" % '-'*57 - - # evaluate stopping criteria - if xStop==[]: - xStop=x0 - Jstop = fctn(xStop) - print "%3d\t%1.2e" % (-1, Jstop[0]) - - # initialize - xc = x0 - STOP = np.zeros((5,1),dtype=bool) - iterLS=0; iter=0 - - Jold = Jstop - xOld=xc - while 1: - # evaluate objective function - Jc,dJ,H = fctn(xc) - print "%3d\t%1.2e\t%1.2e\t%d" % (iter, Jc[0],norm(dJ),iterLS) - - # check stopping rules - STOP[0] = (iter>0) & (abs(Jc[0]-Jold[0]) <= tolJ*(1+abs(Jstop[0]))) - STOP[1] = (iter>0) & (norm(xc-xOld) <= tolX*(1+norm(x0))) - STOP[2] = norm(dJ) <= tolG*(1+abs(Jstop[0])) - STOP[3] = norm(dJ) <= 1e3*eps - STOP[4] = (iter >= maxIter) - if all(STOP[0:3]) | any(STOP[3:]): - break - - # get search direction - dx = np.linalg.solve(H,-dJ) - - # Armijo linesearch - descent = np.dot(dJ.T,dx) - LS =0; t = 1; iterLS=1 - while (iterLS eps] - order1 = order1[E1[1:] > eps] - belowTol = order1.size == 0 and order0.size > 0 - correctOrder = order1.size > 0 and np.mean(order1) > tolerance * expectedOrder - - passTest = belowTol or correctOrder - - if passTest: - print "%s PASS! %s\n" % ('='*25, '='*25) - else: - print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57) - - if plotIt: - plt.figure() - plt.clf() - plt.loglog(t, E0, 'b') - plt.loglog(t, E1, 'g--') - plt.title('checkDerivative') - plt.xlabel('h') - plt.ylabel('error of Taylor approximation') - plt.legend(['0th order', '1st order'], loc='upper left') - plt.show() - - return passTest - if __name__ == '__main__': x0 = np.array([2.6, 3.7]) checkDerivative(Rosenbrock, x0, plotIt=False) diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/TestUtils.py similarity index 56% rename from SimPEG/tests/OrderTest.py rename to SimPEG/tests/TestUtils.py index fb305d32..79340da7 100644 --- a/SimPEG/tests/OrderTest.py +++ b/SimPEG/tests/TestUtils.py @@ -1,5 +1,7 @@ -import sys -sys.path.append('../../') +import numpy as np +import matplotlib.pyplot as plt +from pylab import norm +from SimPEG.utils import mkvc from SimPEG import TensorMesh, utils, LogicallyOrthogonalMesh import numpy as np import unittest @@ -16,46 +18,47 @@ class OrderTest(unittest.TestCase): Given are an operator A and its discretization A[h]. For a given test function f and h --> 0 we compare: - error(h) = \| A[h](f) - A(f) \|_{\infty} + .. math:: + error(h) = \| A[h](f) - A(f) \|_{\infty} Note that you can provide any norm. Test is passed when estimated rate order of convergence is at least within the specified tolerance of the estimated rate supplied by the user. - Minimal example for a curl operator: + Minimal example for a curl operator:: - class TestCURL(OrderTest): - name = "Curl" + class TestCURL(OrderTest): + name = "Curl" - def getError(self): - # For given Mesh, generate A[h], f and A(f) and return norm of error. + def getError(self): + # For given Mesh, generate A[h], f and A(f) and return norm of error. - fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) - sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) + fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) + sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) - Ex = fun(self.M.gridEx[:, 1]) - Ey = fun(self.M.gridEy[:, 2]) - Ez = fun(self.M.gridEz[:, 0]) - f = np.concatenate((Ex, Ey, Ez)) + Ex = fun(self.M.gridEx[:, 1]) + Ey = fun(self.M.gridEy[:, 2]) + Ez = fun(self.M.gridEz[:, 0]) + f = np.concatenate((Ex, Ey, Ez)) - Fx = sol(self.M.gridFx[:, 2]) - Fy = sol(self.M.gridFy[:, 0]) - Fz = sol(self.M.gridFz[:, 1]) - Af = np.concatenate((Fx, Fy, Fz)) + Fx = sol(self.M.gridFx[:, 2]) + Fy = sol(self.M.gridFy[:, 0]) + Fz = sol(self.M.gridFz[:, 1]) + Af = np.concatenate((Fx, Fy, Fz)) - # Generate DIV matrix - Ah = self.M.edgeCurl + # Generate DIV matrix + Ah = self.M.edgeCurl - curlE = Ah*E - err = np.linalg.norm((Ah*f -Af), np.inf) - return err + curlE = Ah*E + err = np.linalg.norm((Ah*f -Af), np.inf) + return err - def test_order(self): - # runs the test - self.orderTest() + def test_order(self): + # runs the test + self.orderTest() See also: test_operatorOrder.py @@ -159,5 +162,78 @@ class OrderTest(unittest.TestCase): print '' self.assertTrue(passTest) -if __name__ == '__main__': - unittest.main() +def Rosenbrock(x): + """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 + +def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): + """ + Basic derivative check + + Compares error decay of 0th and 1st order Taylor approximation at point + x0 for a randomized search direction. + + :param lambda fctn: function handle + :param numpy.array x0: point at which to check derivative + :param int num: number of times to reduce step length, h + :param bool plotIt: if you would like to plot + :param numpy.array dx: step direction + :rtype: bool + :return: did you pass the test?! + + """ + + print "%s checkDerivative %s" % ('='*20, '='*20) + print "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|\tOrder\n%s" % ('-'*57) + + Jc = fctn(x0) + + x0 = mkvc(x0) + + if dx is None: + dx = np.random.randn(len(x0)) + + t = np.logspace(-1, -num, num) + E0 = np.ones(t.shape) + E1 = np.ones(t.shape) + + l2norm = lambda x: np.sqrt(np.inner(x, x)) # because np.norm breaks if they are scalars? + 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 + 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]) + + tolerance = 0.9 + expectedOrder = 2 + eps = 1e-10 + order0 = order0[E0[1:] > eps] + order1 = order1[E1[1:] > eps] + belowTol = order1.size == 0 and order0.size > 0 + correctOrder = order1.size > 0 and np.mean(order1) > tolerance * expectedOrder + + passTest = belowTol or correctOrder + + if passTest: + print "%s PASS! %s\n" % ('='*25, '='*25) + else: + print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57) + + if plotIt: + plt.figure() + plt.clf() + plt.loglog(t, E0, 'b') + plt.loglog(t, E1, 'g--') + plt.title('checkDerivative') + plt.xlabel('h') + plt.ylabel('error of Taylor approximation') + plt.legend(['0th order', '1st order'], loc='upper left') + plt.show() + + return passTest diff --git a/SimPEG/tests/__init__.py b/SimPEG/tests/__init__.py new file mode 100644 index 00000000..f896cfd0 --- /dev/null +++ b/SimPEG/tests/__init__.py @@ -0,0 +1,2 @@ +import TestUtils +from TestUtils import checkDerivative, Rosenbrock, OrderTest diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index 13037156..79626ca1 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -1,6 +1,6 @@ import numpy as np import unittest -from OrderTest import OrderTest +from TestUtils import OrderTest # MATLAB code: diff --git a/SimPEG/tests/test_operators.py b/SimPEG/tests/test_operators.py index ab377089..6906c580 100644 --- a/SimPEG/tests/test_operators.py +++ b/SimPEG/tests/test_operators.py @@ -1,6 +1,6 @@ import numpy as np import unittest -from OrderTest import OrderTest +from TestUtils import OrderTest MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM'] call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/tests/test_tensorMesh.py index 32d47353..300e6634 100644 --- a/SimPEG/tests/test_tensorMesh.py +++ b/SimPEG/tests/test_tensorMesh.py @@ -1,9 +1,7 @@ import numpy as np import unittest -import sys -sys.path.append('../') -from TensorMesh import TensorMesh -from OrderTest import OrderTest +from SimPEG import TensorMesh +from TestUtils import OrderTest from scipy.sparse.linalg import dsolve diff --git a/SimPEG/tests/test_utils.py b/SimPEG/tests/test_utils.py index 0dcdd362..f0b867ec 100644 --- a/SimPEG/tests/test_utils.py +++ b/SimPEG/tests/test_utils.py @@ -1,8 +1,6 @@ import numpy as np import unittest -import sys -sys.path.append('../') -from utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal +from SimPEG.utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal class TestSequenceFunctions(unittest.TestCase): diff --git a/docs/api_GaussNewton.rst b/docs/api_GaussNewton.rst deleted file mode 100644 index 65012347..00000000 --- a/docs/api_GaussNewton.rst +++ /dev/null @@ -1,8 +0,0 @@ -.. _api_GaussNewton: - -Gauss Newton -************ - -.. automodule:: SimPEG.GaussNewton - :members: - :undoc-members: diff --git a/docs/api_Optimize.rst b/docs/api_Optimize.rst new file mode 100644 index 00000000..9765bd3b --- /dev/null +++ b/docs/api_Optimize.rst @@ -0,0 +1,8 @@ +.. _api_Inverse: + +Optimize +******** + +.. automodule:: SimPEG.inverse.Optimize + :members: + :undoc-members: diff --git a/docs/api_Tests.rst b/docs/api_Tests.rst new file mode 100644 index 00000000..164f5b5d --- /dev/null +++ b/docs/api_Tests.rst @@ -0,0 +1,8 @@ +.. _api_Tests: + +Testing SimPEG +************** + +.. automodule:: SimPEG.tests.TestUtils + :members: + :undoc-members: diff --git a/docs/index.rst b/docs/index.rst index 48ebe940..b08624d5 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -30,14 +30,6 @@ Meshing & Operators api_DiffOperators api_InnerProducts -Inversion -========= - -.. toctree:: - :maxdepth: 2 - - api_GaussNewton - Forward Problems ================ @@ -46,6 +38,22 @@ Forward Problems api_Problem +Inversion +========= + +.. toctree:: + :maxdepth: 2 + + api_Optimize + +Testing SimPEG +============== + +.. toctree:: + :maxdepth: 2 + + api_Tests + Project Index & Search ====================== From fb69446cfba65552d7c20628d82f2dbb971085c5 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 3 Oct 2013 15:15:08 -0700 Subject: [PATCH 08/11] Testing problem and DCProblem. --- SimPEG/forward/Problem.py | 9 --- SimPEG/tests/test_forward_DCproblem.py | 81 ++++++++++++++++++++++++++ SimPEG/tests/test_forward_problem.py | 28 +++++++++ 3 files changed, 109 insertions(+), 9 deletions(-) create mode 100644 SimPEG/tests/test_forward_DCproblem.py create mode 100644 SimPEG/tests/test_forward_problem.py diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 7c693d17..1558ecf0 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -327,12 +327,3 @@ class SyntheticProblem(object): eps = np.linalg.norm(mkvc(dobs),2)*1e-5 Wd = 1/(abs(dobs)*std+eps) return dobs, Wd - - - -if __name__ == '__main__': - from SimPEG.inverse import checkDerivative - - p = Problem(None) - m = np.random.rand(5) - checkDerivative(lambda m : [p.modelTransform(m), p.modelTransformDeriv(m)], m, plotIt=False) diff --git a/SimPEG/tests/test_forward_DCproblem.py b/SimPEG/tests/test_forward_DCproblem.py new file mode 100644 index 00000000..b5e6d844 --- /dev/null +++ b/SimPEG/tests/test_forward_DCproblem.py @@ -0,0 +1,81 @@ +import numpy as np +import unittest +from SimPEG import TensorMesh +from SimPEG.utils import ModelBuilder, sdiag +from SimPEG.forward import Problem, SyntheticProblem +from SimPEG.forward.DCProblem import DCProblem, DCutils +from TestUtils import checkDerivative +from scipy.sparse.linalg import dsolve + + +class DCProblemTests(unittest.TestCase): + + def setUp(self): + # Create the mesh + h1 = np.ones(20) + h2 = np.ones(20) + mesh = TensorMesh([h1,h2]) + + # Create some parameters for the model + sig1 = 1 + sig2 = 0.01 + + # Create a synthetic model from a block in a half-space + p0 = [2, 2] + p1 = [5, 5] + condVals = [sig1, sig2] + mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals) + + # Set up the projection + nelec = 10 + spacelec = 2 + surfloc = 0.5 + elecini = 0.5 + 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) + P = Q.T + + # Create some data + class syntheticDCProblem(DCProblem, SyntheticProblem): + pass + + synthetic = syntheticDCProblem(mesh); + synthetic.P = P + synthetic.RHS = q + dobs, Wd = synthetic.createData(mSynth, std=0.05) + + # Now set up the problem to do some minimization + problem = DCProblem(mesh) + problem.P = P + problem.RHS = q + problem.W = Wd + problem.dobs = dobs + + 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)] + passed = checkDerivative(derChk, self.m0, plotIt=False) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC) + 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)) + vtJtw = v.dot(self.p.Jt(self.m0, w, u=u)) + passed = (wtJv - vtJtw) < 1e-10 + self.assertTrue(passed) + + + +if __name__ == '__main__': + unittest.main() diff --git a/SimPEG/tests/test_forward_problem.py b/SimPEG/tests/test_forward_problem.py new file mode 100644 index 00000000..b7180dca --- /dev/null +++ b/SimPEG/tests/test_forward_problem.py @@ -0,0 +1,28 @@ +import numpy as np +import unittest +from SimPEG import TensorMesh +from SimPEG.forward import Problem +from TestUtils import checkDerivative +from scipy.sparse.linalg import dsolve + + +class ProblemTests(unittest.TestCase): + + def setUp(self): + + a = np.array([1, 1, 1]) + b = np.array([1, 2]) + c = np.array([1, 4]) + self.mesh2 = TensorMesh([a, b], np.array([3, 5])) + self.p2 = Problem(self.mesh2) + + + def test_modelTransform(self): + print 'SimPEG.forward.Problem: Testing Model Transform' + m = np.random.rand(self.mesh2.nC) + passed = checkDerivative(lambda m : [self.p2.modelTransform(m), self.p2.modelTransformDeriv(m)], m, plotIt=False) + self.assertTrue(passed) + + +if __name__ == '__main__': + unittest.main() From 65f129e7a84ccb616697f04b865ce6dbfd7af106 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Thu, 3 Oct 2013 15:31:08 -0700 Subject: [PATCH 09/11] Documented utilities. --- SimPEG/utils/ModelBuilder.py | 25 ++++++++++++++----------- SimPEG/utils/__init__.py | 5 ++++- SimPEG/utils/lomutils.py | 27 +++++++++++---------------- SimPEG/utils/matutils.py | 4 ++-- docs/api_Utils.rst | 21 +++++++++++++++++++++ docs/index.rst | 9 +++++++++ 6 files changed, 61 insertions(+), 30 deletions(-) create mode 100644 docs/api_Utils.rst diff --git a/SimPEG/utils/ModelBuilder.py b/SimPEG/utils/ModelBuilder.py index 527d5eef..170da1bf 100644 --- a/SimPEG/utils/ModelBuilder.py +++ b/SimPEG/utils/ModelBuilder.py @@ -3,16 +3,18 @@ import numpy as np def getIndecesBlock(p0,p1,ccMesh): """ - Creates a vector containing the block indexes in the cell centerd mesh. - Returns a tuple + Creates a vector containing the block indexes in the cell centerd mesh. + Returns a tuple - The block is defined by the points - p0 : describe the position of the left upper front corner, and - p1 : describe the position of the right bottom back corner. + The block is defined by the points - ccMesh represents the cell-centered mesh + p0, describe the position of the left upper front corner, and - The points p0 and p1 must live in the the same dimensional space as the mesh. + p1, describe the position of the right bottom back corner. + + ccMesh represents the cell-centered mesh + + The points p0 and p1 must live in the the same dimensional space as the mesh. """ # Validation: p0 and p1 live in the same dimensional space @@ -64,9 +66,9 @@ def getIndecesBlock(p0,p1,ccMesh): def defineBlockConductivity(p0,p1,ccMesh,condVals): """ - Build a block with the conductivity specified by condVal. Returns an array. - condVals[0] conductivity of the block - condVals[1] conductivity of the ground + Build a block with the conductivity specified by condVal. Returns an array. + condVals[0] conductivity of the block + condVals[1] conductivity of the ground """ sigma = np.zeros(ccMesh.shape[0]) + condVals[1] ind = getIndecesBlock(p0,p1,ccMesh) @@ -80,7 +82,8 @@ def defineTwoLayeredConductivity(depth,ccMesh,condVals): Define a two layered model. Depth of the first layer must be specified. CondVals vector with the conductivity values of the layers. Eg: - Convention to number the layers: + Convention to number the layers:: + <----------------------------|------------------------------------> 0 depth zf 1st layer 2nd layer diff --git a/SimPEG/utils/__init__.py b/SimPEG/utils/__init__.py index bb46a4e9..7c976ce3 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -1,4 +1,7 @@ +import matutils +import sputils +import lomutils +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 -import ModelBuilder diff --git a/SimPEG/utils/lomutils.py b/SimPEG/utils/lomutils.py index cdd6e0de..9020ce12 100644 --- a/SimPEG/utils/lomutils.py +++ b/SimPEG/utils/lomutils.py @@ -31,19 +31,18 @@ def volTetra(xyz, A, B, C, D): """ Returns the volume for tetrahedras volume specified by the indexes A to D. + :param numpy.array xyz: X,Y,Z vertex vector + :param numpy.array A,B,C,D: vert index of the tetrahedra + :rtype: numpy.array + :return: V, volume of the tetrahedra - Input: - xyz - X,Y,Z vertex vector - A,B,C,D - vert index of the tetrahedra + Algorithm http://en.wikipedia.org/wiki/Tetrahedron#Volume - Output: - V - volume + .. math:: - Algorithm: http://en.wikipedia.org/wiki/Tetrahedron#Volume + V = {1 \over 3} A h - V = 1/3 A * h - - V = 1/6 | ( a - d ) o ( ( b - d ) X ( c - d ) ) | + V = {1 \over 6} | ( a - d ) \cdot ( ( b - d ) ( c - d ) ) | """ @@ -69,7 +68,7 @@ def indexCube(nodes, gridSize, n=None): Output: index - index in the order asked e.g. 'ABCD' --> (A,B,C,D) - TWO DIMENSIONS: + TWO DIMENSIONS:: node(i,j) node(i,j+1) A -------------- B @@ -81,7 +80,7 @@ def indexCube(nodes, gridSize, n=None): node(i+1,j) node(i+1,j+1) - THREE DIMENSIONS: + THREE DIMENSIONS:: node(i,j,k+1) node(i,j+1,k+1) E --------------- F @@ -97,10 +96,6 @@ def indexCube(nodes, gridSize, n=None): D -------------- C node(i+1,j,k) node(i+1,j+1,k) - - @author Rowan Cockett - - Last modified on: 2013/07/26 """ assert type(nodes) == str, "Nodes must be a str variable: e.g. 'ABCD'" @@ -211,7 +206,7 @@ def faceInfo(xyz, A, B, C, D, average=True, normalizeNormals=True): # # So also could be viewed as the average parallelogram. # - # WARNING: This does not compute correctly for concave quadrilaterals + # TODO: This does not compute correctly for concave quadrilaterals area = (length(nA)+length(nB)+length(nC)+length(nD))/4 return N, area diff --git a/SimPEG/utils/matutils.py b/SimPEG/utils/matutils.py index 20cc33df..24286cdd 100644 --- a/SimPEG/utils/matutils.py +++ b/SimPEG/utils/matutils.py @@ -4,7 +4,7 @@ import numpy as np def mkvc(x, numDims=1): """Creates a vector with the number of dimension specified - e.g.: + e.g.:: a = np.array([1, 2, 3]) @@ -43,7 +43,7 @@ def ndgrid(*args, **kwargs): The inputs can be a list or separate arguments. - e.g. + e.g.:: a = np.array([1, 2, 3]) b = np.array([1, 2]) diff --git a/docs/api_Utils.rst b/docs/api_Utils.rst new file mode 100644 index 00000000..e952010b --- /dev/null +++ b/docs/api_Utils.rst @@ -0,0 +1,21 @@ +.. _api_Utils: + +Utilities +********* + +.. automodule:: SimPEG.utils.matutils + :members: + :undoc-members: + +.. automodule:: SimPEG.utils.sputils + :members: + :undoc-members: + +.. automodule:: SimPEG.utils.lomutils + :members: + :undoc-members: + +.. automodule:: SimPEG.utils.ModelBuilder + :members: + :undoc-members: + diff --git a/docs/index.rst b/docs/index.rst index b08624d5..4c01c4cd 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -55,6 +55,15 @@ Testing SimPEG api_Tests +Utility Codes +============= + +.. toctree:: + :maxdepth: 2 + + api_Utils + + Project Index & Search ====================== From 0089931075422708228de96062993ee9bbbb727f Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 7 Oct 2013 15:15:18 -0700 Subject: [PATCH 10/11] bug fixes for optimization --- SimPEG/inverse/Optimize.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index df6aa9a5..eee18b16 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -20,7 +20,9 @@ class Minimize(object): def __init__(self, problem, **kwargs): self.problem = problem + 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): @@ -130,6 +132,7 @@ class SteepestDescent(Minimize): return -self.g 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) From 4851295738e9aa30b2bdfcd043ee6a79032a27ab Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Wed, 16 Oct 2013 09:18:08 -0700 Subject: [PATCH 11/11] initial solver work: import SimPEG.Solver --- SimPEG/Solver.py | 75 +++++++++++++++++++++++++++++++++++++++ SimPEG/__init__.py | 1 + SimPEG/forward/Problem.py | 10 ++++++ 3 files changed, 86 insertions(+) create mode 100644 SimPEG/Solver.py diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py new file mode 100644 index 00000000..2087ef5f --- /dev/null +++ b/SimPEG/Solver.py @@ -0,0 +1,75 @@ +import numpy as np +import scipy.sparse.linalg as linalg + + +class Solver(object): + """docstring for Solver""" + 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'" + + self.A = A + + self.dsolve = None + self.doDirect = doDirect + self.flag = flag + self.options = options + + def solve(self, b): + if self.flag is None and self.doDirect: + return self.solveDirect(b, **self.options) + elif self.flag is None and not self.doDirect: + return self.solveIter(b, **self.options) + elif self.flag == 'U': + return self.solveBackward(b) + elif self.flag == 'L': + return self.solveForward(b) + elif self.flag == 'D': + return self.solveDiagonal(b) + else: + raise Exception('Unknown flag.') + pass + + def clean(self): + """Cleans up the memory""" + del self.dsolve + self.dsolve = None + + def solveDirect(self, b, backend='scipy'): + assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch' + + if 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) + + # Multiple RHSs + X = np.empty_like(b) + for i in range(b.shape[1]): + X[:,i] = self.dsolve(b[:,i]) + + return X + + def solveIter(self, b, M=None, iterSolver='CG'): + pass + + def solveBackward(self, b): + pass + + def solveForward(self, b): + pass + + def solveDiagonal(self, b): + diagA = self.A.diagonal() + if len(b.shape) == 1 or b.shape[1] == 1: + # Just one RHS + return b/diagA + # Multiple RHSs + X = np.empty_like(b) + for i in range(b.shape[1]): + X[:,i] = b[:,i]/diagA + return X diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index f36000a7..ea9faf33 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -2,3 +2,4 @@ from TensorMesh import TensorMesh from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh import utils import inverse +from Solver import Solver diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 1558ecf0..5b716f1f 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -83,6 +83,16 @@ class Problem(object): def dobs(self, value): self._dobs = value + def evalFunction(self, m, doDerivative=True): + """ + :param numpy.array m: model + :param bool doDerivative: do you want to compute the derivative? + :rtype: numpy.array + :return: Jv + """ + f = self.misfit(m) + + return f, g, H def J(self, m, v, u=None): """