diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index 5aad45bd..5a011bb1 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -3,7 +3,7 @@ from SimPEG.forward import Problem, SyntheticProblem from SimPEG.tests import checkDerivative from SimPEG.utils import ModelBuilder, sdiag import numpy as np -import scipy as sp +import scipy.sparse as sp import scipy.sparse.linalg as linalg class DCProblem(Problem): @@ -104,14 +104,44 @@ class DCProblem(Problem): return Jtv + +def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh): + """ Generate projection matrix (Q) and """ + elecend = 0.5+spacelec*(nelec-1) + elecLocR = np.linspace(elecini, elecend, nelec) + elecLocT = elecLocR+1 + nrx = nelec-1 + ntx = nelec-1 + q = np.zeros((mesh.nC, ntx)) + Q = np.zeros((mesh.nC, nrx)) + + for i in range(nrx): + + rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i])) + rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1])) + + txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i])) + txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1])) + + q[txind1,i] = 1 + q[txind2,i] = -1 + Q[rxind1,i] = 1 + Q[rxind2,i] = -1 + + Q = sp.csr_matrix(Q) + rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5 + return q, Q, rxmidLoc + + + if __name__ == '__main__': from SimPEG.regularization import Regularization from SimPEG import inverse # Create the mesh - h1 = np.ones(100) - h2 = np.ones(100) + h1 = np.ones(10) + h2 = np.ones(10) mesh = TensorMesh([h1,h2]) # Create some parameters for the model @@ -153,14 +183,14 @@ if __name__ == '__main__': problem = DCProblem(mesh) problem.P = P problem.RHS = q - problem.W = Wd problem.dobs = dobs + problem.std = dobs*0 + 0.05 m0 = mesh.gridCC[:,0]*0+sig1 # print problem.misfit(m0) # print problem.misfit(mSynth) - opt = inverse.InexactGaussNewton() + opt = inverse.InexactGaussNewton(maxIterLS=20) reg = Regularization(mesh) inv = inverse.Inversion(problem, reg, opt) @@ -181,30 +211,3 @@ if __name__ == '__main__': - -def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh): - """ Generate projection matrix (Q) and """ - elecend = 0.5+spacelec*(nelec-1) - elecLocR = np.linspace(elecini, elecend, nelec) - elecLocT = elecLocR+1 - nrx = nelec-1 - ntx = nelec-1 - q = np.zeros((mesh.nC, ntx)) - Q = np.zeros((mesh.nC, nrx)) - - for i in range(nrx): - - rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i])) - rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1])) - - txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i])) - txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1])) - - q[txind1,i] = 1 - q[txind2,i] = -1 - Q[rxind1,i] = 1 - Q[rxind2,i] = -1 - - Q = sp.csr_matrix(Q) - rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5 - return q, Q, rxmidLoc diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 8d92e15e..9fd10abd 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,5 +1,6 @@ import numpy as np -from SimPEG.utils import sdiag +import scipy.sparse as sp +from SimPEG.utils import sdiag, mkvc class Inversion(object): """docstring for Inversion""" @@ -16,7 +17,10 @@ class Inversion(object): """ Standard deviation weighting matrix. """ - return sdiag(1/self.prob.std) + if getattr(self,'_Wd',None) is None: + eps = np.linalg.norm(mkvc(self.prob.dobs),2)*1e-5 + self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps) + return self._Wd def run(self, m0): m = m0 @@ -41,7 +45,7 @@ class Inversion(object): u = self.prob.field(m) phi_d = self.dataObj(m, u) - phi_m = self.modelObj(m) + phi_m = self.reg.modelObj(m) self._phi_d_last = phi_d self._phi_m_last = phi_m @@ -50,27 +54,24 @@ class Inversion(object): out = (f,) if return_g: - phi_dDeriv = self.dataObjDeriv(m, u) - phi_mDeriv = self.modelObjDeriv(m) + phi_dDeriv = self.dataObjDeriv(m, u=u) + phi_mDeriv = self.reg.modelObjDeriv(m) g = phi_dDeriv + self._beta * phi_mDeriv out += (g,) if return_H: def H_fun(v): - phi_d2Deriv = self.dataObj2Deriv(m, u, v) - phi_m2Deriv = self.modelObj2Deriv(m)*v + phi_d2Deriv = self.dataObj2Deriv(m, v, u=u) + phi_m2Deriv = self.reg.modelObj2Deriv(m)*v return phi_d2Deriv + self._beta * phi_m2Deriv - out += (H_fun,) + operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=float ) + out += (operator,) return out - def modelObj(self, m, u=None): - self.reg.misfit(m) - - def dataObj(self, m, u=None): """ :param numpy.array m: geophysical model @@ -87,7 +88,7 @@ class Inversion(object): Where P is a projection matrix that brings the field on the full domain to the data measurement locations; u is the field of interest; d_obs is the observed data; and W is the weighting matrix. """ - R = self.Wd*self.prob.misfit(u=u) + R = self.Wd*self.prob.misfit(m, u=u) R = mkvc(R) return 0.5*R.dot(R) @@ -123,17 +124,17 @@ class Inversion(object): """ if u is None: - u = self.field(m) + u = self.prob.field(m) - R = self.W*(self.dpred(m, u=u) - self.dobs) + R = self.Wd*self.prob.misfit(m, u=u) dmisfit = 0 - for i in range(self.RHS.shape[1]): # Loop over each right hand side - dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i]) + for i in range(self.prob.RHS.shape[1]): # Loop over each right hand side + dmisfit += self.prob.Jt(m, self.Wd[:,i]*R[:,i], u=u[:,i]) return dmisfit - def dataObj2Deriv(self, m, u=None): + def dataObj2Deriv(self, m, v, u=None): """ :param numpy.array m: geophysical model :param numpy.array u: fields @@ -167,12 +168,12 @@ class Inversion(object): """ if u is None: - u = self.field(m) + u = self.prob.field(m) - R = self.W*(self.dpred(m, u=u) - self.dobs) + R = self.Wd*self.prob.misfit(m, u=u) dmisfit = 0 - for i in range(self.RHS.shape[1]): # Loop over each right hand side - dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i]) + for i in range(self.prob.RHS.shape[1]): # Loop over each right hand side + dmisfit += self.prob.Jt(m, self.Wd[:,i] * self.Wd[:,i] * self.prob.J(m, v, u=u[:,i]), u=u[:,i]) return dmisfit diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index 43776d33..1bfda72c 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -2,6 +2,7 @@ import numpy as np import matplotlib.pyplot as plt from SimPEG.utils import mkvc, sdiag norm = np.linalg.norm +import scipy.sparse as sp class Minimize(object): @@ -149,7 +150,8 @@ class GaussNewton(Minimize): class InexactGaussNewton(Minimize): name = 'InexactGaussNewton' def findSearchDirection(self): - return sparse.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=10) + p, info = sp.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=10) + return p class SteepestDescent(Minimize): diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py index 110fe9bd..380e0d19 100644 --- a/SimPEG/mesh/DiffOperators.py +++ b/SimPEG/mesh/DiffOperators.py @@ -161,6 +161,68 @@ class DiffOperators(object): _cellGrad = None cellGrad = property(**cellGrad()) + def cellGradx(): + doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." + + def fget(self): + if getattr(self, '_cellGradx', None) is None: + BC = ['neumann', 'neumann'] + n = self.n + if(self.dim == 1): + G1 = ddxCellGrad(n[0], BC) + elif(self.dim == 2): + G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC)) + elif(self.dim == 3): + G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC)) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F','Fx', 'V') + V = self.vol + self._cellGradx = sdiag(S)*G1*sdiag(1/V) + return self._cellGradx + return locals() + cellGradx = property(**cellGradx()) + + + def cellGrady(): + doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." + def fget(self): + if self.dim < 2: + return None + if getattr(self, '_cellGrady', None) is None: + BC = ['neumann', 'neumann'] + n = self.n + if(self.dim == 2): + G2 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC)) + elif(self.dim == 3): + G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0])) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F','Fy', 'V') + V = self.vol + self._cellGrady = sdiag(S)*G2*sdiag(1/V) + return self._cellGrady + return locals() + cellGrady = property(**cellGrady()) + + + + def cellGradz(): + doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." + def fget(self): + if self.dim < 3: + return None + if getattr(self, '_cellGradz', None) is None: + BC = ['neumann', 'neumann'] + n = self.n + G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0])) + # Compute areas of cell faces & volumes + S = self.r(self.area, 'F','Fz', 'V') + V = self.vol + self._cellGradz = sdiag(S)*G3*sdiag(1/V) + return self._cellGradz + return locals() + cellGradz = property(**cellGradz()) + + def edgeCurl(): doc = "Construct the 3D curl operator." diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py index f0239875..7e211df6 100644 --- a/SimPEG/regularization/Regularization.py +++ b/SimPEG/regularization/Regularization.py @@ -1,10 +1,13 @@ from SimPEG.utils import sdiag +import numpy as np class Regularization(object): """docstring for Regularization""" @property def mref(self): + if getattr(self, '_mref', None) is None: + self._mref = np.zeros(self.mesh.nC); return self._mref @mref.setter def mref(self, value): @@ -12,25 +15,25 @@ class Regularization(object): @property def Wx(self): - if self._Wx is None: - self._Wx = mesh.cellGradx + if getattr(self, '_Wx', None) is None: + self._Wx = self.mesh.cellGradx return self._Wx @property def Wy(self): - if self._Wy is None: - self._Wy = mesh.cellGrady + if getattr(self, '_Wy', None) is None: + self._Wy = self.mesh.cellGrady return self._Wy @property def Wz(self): - if self._Wz is None: - self._Wz = mesh.cellGradz + if getattr(self, '_Wz', None) is None: + self._Wz = self.mesh.cellGradz return self._Wz @property def Ws(self): - if self._Ws is None: + if getattr(self,'_Ws', None) is None: self._Ws = sdiag(self.mesh.vol) return self._Ws @@ -87,12 +90,12 @@ class Regularization(object): mobjDeriv = self.alpha_s * self.Ws.T * ( self.Ws * mresid) - mobjDeriv += self.alpha_x * self.Wx.T * ( self.Wx * mresid) + mobjDeriv = mobjDeriv + self.alpha_x * self.Wx.T * ( self.Wx * mresid) if self.mesh.dim > 1: - mobjDeriv += self.alpha_y * self.Wy.T * ( self.Wy * mresid) + mobjDeriv = mobjDeriv + self.alpha_y * self.Wy.T * ( self.Wy * mresid) if self.mesh.dim > 2: - mobjDeriv += self.alpha_z * self.Wz.T * ( self.Wz * mresid) + mobjDeriv = mobjDeriv + self.alpha_z * self.Wz.T * ( self.Wz * mresid) return mobjDeriv @@ -102,12 +105,12 @@ class Regularization(object): mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws - mobj2Deriv += self.alpha_x * self.Wx.T * self.Wx + mobj2Deriv = mobj2Deriv + self.alpha_x * self.Wx.T * self.Wx if self.mesh.dim > 1: - mobj2Deriv += self.alpha_y * self.Wy.T * self.Wy + mobj2Deriv = mobj2Deriv + self.alpha_y * self.Wy.T * self.Wy if self.mesh.dim > 2: - mobj2Deriv += self.alpha_z * self.Wz.T * self.Wz + mobj2Deriv = mobj2Deriv + self.alpha_z * self.Wz.T * self.Wz return mobj2Deriv