diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py deleted file mode 100644 index 2087ef5f..00000000 --- a/SimPEG/Solver.py +++ /dev/null @@ -1,75 +0,0 @@ -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 e8946f88..7f059a74 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -1,4 +1,6 @@ -import mesh import utils +from utils import Solver +import mesh import inverse -from Solver import Solver +import forward +import regularization diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py new file mode 100644 index 00000000..132966a8 --- /dev/null +++ b/SimPEG/forward/DCProblem.py @@ -0,0 +1,249 @@ +from SimPEG.mesh import TensorMesh +from SimPEG.forward import Problem, SyntheticProblem, ModelTransforms +from SimPEG.tests import checkDerivative +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 + + +class DCProblem(ModelTransforms.LogModel, Problem): + """ + **DCProblem** + + Geophysical DC resistivity problem. + + """ + def __init__(self, mesh): + 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. + + :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) + Msig = self.mesh.getFaceMass(sigma) + 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 mkvc(phi) + + def J(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 + + .. 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))) + + 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 ) ) + """ + if u is None: + u = self.field(m) + + u = self.reshapeFields(u) + + 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 + + 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 ) ) ) + + solve = Solver(dCdu) + Jv = - P * solve.solve(dCdm) + 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 + A = self.createMatrix(m) + Av_dm = self.mesh.getFaceMassDeriv() + mT_dm = self.modelTransformDeriv(m) + + dCdu = A.T + solve = Solver(dCdu) + + w = solve.solve(P.T*v) + + 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 + + + +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 + import matplotlib.pyplot as plt + + # Create the mesh + h1 = np.ones(20) + h2 = np.ones(100) + mesh = TensorMesh([h1,h2]) + + # Create some parameters for the model + sig1 = np.log(1) + sig2 = np.log(0.01) + + # Create a synthetic model from a block in a half-space + p0 = [5, 10] + p1 = [15, 50] + condVals = [sig1, sig2] + mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals) + plt.colorbar(mesh.plotImage(mSynth)) + plt.show() + + # 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 = 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) + + u = synthetic.field(mSynth) + u = synthetic.reshapeFields(u) + mesh.plotImage(u[:,10]) + # plt.show() + + # Now set up the problem to do some minimization + problem = DCProblem(mesh) + problem.P = P + problem.RHS = q + problem.dobs = dobs + problem.std = dobs*0 + 0.05 + m0 = mesh.gridCC[:,0]*0+sig2 + + 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) + + # Check Derivative + derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)] + checkDerivative(derChk, mSynth) + + + + 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/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py deleted file mode 100644 index 1df8897b..00000000 --- a/SimPEG/forward/DCProblem/DCProblem.py +++ /dev/null @@ -1,168 +0,0 @@ -from SimPEG.mesh import TensorMesh -from SimPEG.forward import Problem, SyntheticProblem -from SimPEG.tests import checkDerivative -from SimPEG.utils import ModelBuilder, sdiag -import numpy as np -import scipy.sparse.linalg as linalg -import DCutils - -class DCProblem(Problem): - """ - **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) - 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, 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))) - - 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 - G = self.mesh.cellGrad - A = self.createMatrix(m) - Av_dm = self.mesh.getFaceMassDeriv() - 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) - - 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__': - # 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) - 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) - - 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 - 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/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/LinearProblem.py b/SimPEG/forward/LinearProblem.py new file mode 100644 index 00000000..d30a5b4d --- /dev/null +++ b/SimPEG/forward/LinearProblem.py @@ -0,0 +1,89 @@ +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 + + +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) + +if __name__ == '__main__': + 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') + + + + + + 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/ModelTransforms.py b/SimPEG/forward/ModelTransforms.py new file mode 100644 index 00000000..ea89b974 --- /dev/null +++ b/SimPEG/forward/ModelTransforms.py @@ -0,0 +1,49 @@ +import numpy as np +from SimPEG.utils import mkvc, sdiag + +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 5b716f1f..2e6831f7 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 @@ -49,16 +50,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): """ @@ -72,6 +63,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): @@ -83,16 +83,35 @@ class Problem(object): def dobs(self, value): self._dobs = value - def evalFunction(self, m, doDerivative=True): + def dpred(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) + Predicted data. - return f, g, H + .. math:: + d_\\text{pred} = Pu(m) + """ + if u is None: + u = self.field(m) + return self.P*u + + def dataResidual(self, m, u=None): + """ + :param numpy.array m: geophysical model + :param numpy.array u: fields + :rtype: float + :return: data misfit + + 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): """ @@ -131,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. @@ -145,17 +192,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 @@ -168,13 +204,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): """ @@ -184,129 +215,10 @@ 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) - 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): @@ -337,3 +249,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 fe849d41..33c9a6b1 100644 --- a/SimPEG/forward/__init__.py +++ b/SimPEG/forward/__init__.py @@ -1,2 +1,4 @@ from Problem import * import DCProblem +from LinearProblem import LinearProblem +import ModelTransforms diff --git a/SimPEG/inverse/BetaSchedule.py b/SimPEG/inverse/BetaSchedule.py new file mode 100644 index 00000000..fe197340 --- /dev/null +++ b/SimPEG/inverse/BetaSchedule.py @@ -0,0 +1,12 @@ + + +class Cooling(object): + """Simple Beta Schedule""" + + 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 new file mode 100644 index 00000000..3aa8ce58 --- /dev/null +++ b/SimPEG/inverse/Inversion.py @@ -0,0 +1,221 @@ +import numpy as np +import scipy.sparse as sp +from SimPEG.utils import sdiag, mkvc + +class Inversion(object): + """docstring for Inversion""" + + maxIter = 10 + name = 'SimPEG Inversion' + + 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): + """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]) + 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): + """ + Standard deviation weighting matrix. + """ + 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 + + @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._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 + self._beta = None + while True: + self._beta = self.getBeta() + m = self.opt.minimize(self.evalFunction,m) + if self.stoppingCriteria(): break + self._iter += 1 + return m + + beta0 = 1.e2 + beta_coolingFactor = 5. + + def getBeta(self): + if self._beta is None: + return self.beta0 + return self._beta / self.beta_coolingFactor + + def stoppingCriteria(self): + self._STOP = np.zeros(2,dtype=bool) + self._STOP[0] = self._iter >= self.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.reg.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=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, v, u=u) + phi_m2Deriv = self.reg.modelObj2Deriv(m)*v + + return phi_d2Deriv + self._beta * phi_m2Deriv + + operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=float ) + out += (operator,) + return out + + + 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. + """ + # 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*np.vdot(R, 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.prob.field(m) + + R = self.Wd*self.prob.dataResidual(m, u=u) + + dmisfit = self.prob.Jt(m, self.Wd * R, u=u) + + return dmisfit + + def dataObj2Deriv(self, m, v, 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.prob.field(m) + + 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_approx(m, self.Wd * self.Wd * self.prob.J_approx(m, v, u=u), u=u) + + return dmisfit + diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index eee18b16..6221246f 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -2,54 +2,154 @@ import numpy as np 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 + 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): - """docstring for Minimize""" + """ + + Minimize is a general class for derivative based optimization. + + + """ name = "GeneralOptimizationAlgorithm" maxIter = 20 maxIterLS = 10 + maxStep = np.inf LSreduction = 1e-4 LSshorten = 0.5 - tolF = 1e-4 - tolX = 1e-4 - tolG = 1e-4 - eps = 1e-16 + tolF = 1e-1 + tolX = 1e-1 + tolG = 1e-1 + eps = 1e-5 - def __init__(self, problem, **kwargs): - self.problem = problem + 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): - # 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]) else: raise Exception('%s attr is not recognized' % attr) - def minimize(self, x0): + 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:: + + (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) 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) + 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() - xt, passLS = self.linesearch(p) + if doPub: pub.sendMessage('Minimize.searchDirection', minimize=self, p=p) + 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, passLS=passLS) if not passLS: - xt = self.linesearchBreak(p) + xt, caught = self.modifySearchDirectionBreak(p) + if not caught: return self.xc self.doEndIteration(xt) + if doPub: pub.sendMessage('Minimize.endIteration', minimize=self, xt=xt) self.printDone() 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): + """ + **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) @@ -59,29 +159,57 @@ class Minimize(object): self.xOld = x0 def printInit(self): + """ + **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) + 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 def printIter(self): + """ + **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'): + 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): + """ + **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() + 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) 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 @@ -94,14 +222,87 @@ class Minimize(object): self._STOP[4] = self._iter >= self.maxIter return all(self._STOP[0:3]) | any(self._STOP[3:]) - def linesearch(self, p): + 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 iterLS = 0 while iterLS < self.maxIterLS: - xt = self.xc + t*p - ft, temp, temp = self.evalFunction(xt, doDerivative=False) + 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 iterLS += 1 @@ -110,10 +311,37 @@ class Minimize(object): self._iterLS = iterLS return xt, iterLS < self.maxIterLS - def linesearchBreak(self, p): - raise Exception('The linesearch got broken. Boo.') + 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 + evalFunction 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): + """ + **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 @@ -123,7 +351,19 @@ 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=self.tolCG, maxiter=self.maxIterCG) + return p class SteepestDescent(Minimize): @@ -133,18 +373,15 @@ 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(Rosenbrock, maxIter=20).minimize(x0) + + def listener1(minimize,p): + print 'hi: ', p + 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]) - xOpt = SteepestDescent(Rosenbrock, maxIter=20, maxIterLS=15).minimize(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): - 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 index b2a5e506..14bddce7 100644 --- a/SimPEG/inverse/__init__.py +++ b/SimPEG/inverse/__init__.py @@ -1 +1,3 @@ from Optimize import * +from Inversion import * +import BetaSchedule 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/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py index 915bf0ef..93b82b25 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/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 110fe9bd..598b392a 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(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 + 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/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) 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 cee23f7e..ea3a0200 100644 --- a/SimPEG/mesh/TensorMesh.py +++ b/SimPEG/mesh/TensorMesh.py @@ -39,7 +39,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) 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: diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py new file mode 100644 index 00000000..6f5970e6 --- /dev/null +++ b/SimPEG/regularization/Regularization.py @@ -0,0 +1,120 @@ +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): + 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: + 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: + 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: + a = self.mesh.r(self.mesh.area,'F','Fz','V') + self._Wz = sdiag(a)*self.mesh.cellGradz + return self._Wz + + + + 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 = mobjDeriv + self.alpha_x * self.Wx.T * ( self.Wx * mresid) + + if self.mesh.dim > 1: + mobjDeriv = mobjDeriv + self.alpha_y * self.Wy.T * ( self.Wy * mresid) + if self.mesh.dim > 2: + mobjDeriv = 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 = mobj2Deriv + self.alpha_x * self.Wx.T * self.Wx + + if self.mesh.dim > 1: + mobj2Deriv = mobj2Deriv + self.alpha_y * self.Wy.T * self.Wy + if self.mesh.dim > 2: + mobj2Deriv = mobj2Deriv + self.alpha_z * self.Wz.T * self.Wz + + return mobj2Deriv + 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 diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index ee8f90e4..140b56fa 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): """ @@ -159,19 +162,26 @@ 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) -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): """ @@ -188,6 +198,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) @@ -208,7 +228,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]) @@ -224,9 +248,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() @@ -240,3 +267,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_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() 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() diff --git a/SimPEG/utils/Solver.py b/SimPEG/utils/Solver.py new file mode 100644 index 00000000..d9ac4a1d --- /dev/null +++ b/SimPEG/utils/Solver.py @@ -0,0 +1,207 @@ +import numpy as np +import scipy.sparse as sparse +import scipy.sparse.linalg as linalg + + +class Solver(object): + """ + 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'" + + self.A = A + + self.dsolve = None + self.doDirect = doDirect + self.flag = flag + 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: + 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, 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: + 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 + 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]): + if factorize: + X[:,i] = self.dsolve(b[:,i]) + else: + X[:,i] = linalg.dsolve.spsolve(self.A,b[:,i]) + + return X + + def solveIter(self, b, M=None, iterSolver='CG'): + pass + + 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) + 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, 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) + 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, 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 + return b/diagA + # Multiple RHSs + X = np.empty_like(b) + 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 + A[0,0] *= 10 # remove the constant null space from the matrix + + 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/utils/__init__.py b/SimPEG/utils/__init__.py index 9457df03..ae83222e 100644 --- a/SimPEG/utils/__init__.py +++ b/SimPEG/utils/__init__.py @@ -3,7 +3,9 @@ import sputils import lomutils import interputils 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 -from interputils import interpmat \ No newline at end of file +from interputils import interpmat 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_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: + diff --git a/docs/api_Solver.rst b/docs/api_Solver.rst new file mode 100644 index 00000000..db75860b --- /dev/null +++ b/docs/api_Solver.rst @@ -0,0 +1,9 @@ +.. _api_Solver: + +Solver +****** + +.. automodule:: SimPEG.utils.Solver + :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/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]) diff --git a/docs/index.rst b/docs/index.rst index 692363a2..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 @@ -62,6 +55,7 @@ Utility Codes .. toctree:: :maxdepth: 2 + api_Solver api_Utils diff --git a/docs/requirements.txt b/docs/requirements.txt index 24ce15ab..dcecbe99 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1 +1,2 @@ numpy +pypubsub