From 8803956d8301747f182fb2f23c70b8c324fc76c6 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Mon, 9 May 2016 19:58:56 +0900 Subject: [PATCH 1/6] Working on SIP --- SimPEG/EM/Static/SIP/ProblemSIP.py | 422 +++++++++++++++++++++++++++++ SimPEG/EM/Static/SIP/RxSIP.py | 82 ++++++ SimPEG/EM/Static/SIP/SurveySIP.py | 32 +++ SimPEG/EM/Static/SIP/__init__.py | 2 + SimPEG/EM/Static/__init__.py | 1 + 5 files changed, 539 insertions(+) create mode 100644 SimPEG/EM/Static/SIP/ProblemSIP.py create mode 100644 SimPEG/EM/Static/SIP/RxSIP.py create mode 100644 SimPEG/EM/Static/SIP/SurveySIP.py create mode 100644 SimPEG/EM/Static/SIP/__init__.py diff --git a/SimPEG/EM/Static/SIP/ProblemSIP.py b/SimPEG/EM/Static/SIP/ProblemSIP.py new file mode 100644 index 00000000..d01ae706 --- /dev/null +++ b/SimPEG/EM/Static/SIP/ProblemSIP.py @@ -0,0 +1,422 @@ +from SimPEG import Problem, Utils, Maps, Mesh +from SimPEG.EM.Base import BaseEMProblem +from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N +from SimPEG.Utils import sdiag +import numpy as np +from SimPEG.Utils import Zero +from SimPEG.EM.Static.DC import getxBCyBC_CC +from SurveySIP import Survey + +class ColeColePropMap(Maps.PropMap): + """ + Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m) + """ + + eta = Maps.Property("Electrical Conductivity", defaultInvProp=True) + tau = Maps.Property("Electrical Conductivity", defaultVal=0.1, propertyLink=('taui', Maps.ReciprocalMap)) + taui = Maps.Property("Electrical Conductivity", defaultVal=1., propertyLink=('tau', Maps.ReciprocalMap)) + c = Maps.Property("Electrical Conductivity", defaultVal=1.) + + +class BaseSIPProblem(BaseEMProblem): + + surveyPair = Survey + fieldsPair = Fields + PropMap = ColeColePropMap + Ainv = None + sigma = None + rho = None + f = None + Ainv = None + + def DebyeTime(t): + peta = self.curModel.eta*np.exp(-self.curModel.taui*t) + return peta + + def EtaDeriv(t): + return np.exp(-self.curModel.taui*t) + + def TauiDeriv(t): + return -self.curModel.eta*t*np.exp(-self.curModel.taui*t) + + def fields(self, m): + self.curModel = m + if self.f is None: + self.f = self.fieldsPair(self.mesh, self.survey) + if self.Ainv == None: + A = self.getA() + self.Ainv = self.Solver(A, **self.solverOpts) + RHS = self.getRHS() + u = self.Ainv * RHS + Srcs = self.survey.srcList + self.f[Srcs, self._solutionType] = u + return self.f + + def forward(self, m, f=None): + + if f is None: + f = self.fields(m) + + self.curModel = m + Jv = self.dataPair(self.survey) #same size as the data + # A = self.getA() + JvAll = [] + for tind in range(len(self.survey.times)): + #Pseudo-chareability + v = DebyeTime(self.survey.times[tind]) + for src in self.survey.srcList: + u_src = f[src, self._solutionType] # solution vector + dA_dm_v = self.getADeriv(u_src, v) + dRHS_dm_v = self.getRHSDeriv(src, v) + du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v ) + for rx in src.rxList: + timeindex = rx.getTimeP(self.survey.times) + if timeindex[tind]: + df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) + Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) + JvAll.append(Utils.mkvc(Jv)) + # Conductivity (d u / d log sigma) + if self._formulation is 'EB': + return -np.hstack(JvAll) + # Conductivity (d u / d log rho) + if self._formulation is 'HJ': + return np.hstack(JvAll) + + # def Jvec(self, m, v, f=None): + + # if f is None: + # f = self.fields(m) + + # self.curModel = m + + # Jv = self.dataPair(self.survey) #same size as the data + # x1 = + # x2 = + # # A = self.getA() + # for src in self.survey.srcList: + # u_src = f[src, self._solutionType] # solution vector + + # dA_dm_v = self.getADeriv(u_src, v) + # dRHS_dm_v = self.getRHSDeriv(src, v) + # du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v ) + + # for rx in src.rxList: + + # for tind in range(len(self.survey.times)): + # df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + # df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) + # Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) + # # Conductivity (d u / d log sigma) + # if self._formulation is 'EB': + # return -Utils.mkvc(Jv) + # # Conductivity (d u / d log rho) + # if self._formulation is 'HJ': + # return Utils.mkvc(Jv) + + # def Jtvec(self, m, v, f=None): + # if f is None: + # f = self.fields(m) + + # self.curModel = m + + # # Ensure v is a data object. + # if not isinstance(v, self.dataPair): + # v = self.dataPair(self.survey, v) + + # Jtv = np.zeros(m.size) + # # AT = self.getA() + + # for src in self.survey.srcList: + # u_src = f[src, self._solutionType] + # for rx in src.rxList: + # PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m + # df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + # df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) + # ATinvdf_duT = self.Ainv * df_duT + # dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) + # dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) + # du_dmT = -dA_dmT + dRHS_dmT + # Jtv += df_dmT + du_dmT + # # Conductivity ((d u / d log sigma).T) + # if self._formulation is 'EB': + # return -Utils.mkvc(Jtv) + # # Conductivity ((d u / d log rho).T) + # if self._formulation is 'HJ': + # return Utils.mkvc(Jtv) + + def getSourceTerm(self): + """ + takes concept of source and turns it into a matrix + """ + """ + Evaluates the sources, and puts them in matrix form + + :rtype: (numpy.ndarray, numpy.ndarray) + :return: q (nC or nN, nSrc) + """ + + Srcs = self.survey.srcList + + if self._formulation is 'EB': + n = self.mesh.nN + # return NotImplementedError + + elif self._formulation is 'HJ': + n = self.mesh.nC + + q = np.zeros((n, len(Srcs))) + + for i, src in enumerate(Srcs): + q[:,i] = src.eval(self) + return q + + @property + def deleteTheseOnModelUpdate(self): + toDelete = [] + return toDelete + + # assume log rho or log cond + @property + def MeSigma(self): + """ + Edge inner product matrix for \\(\\sigma\\). Used in the E-B formulation + """ + if getattr(self, '_MeSigma', None) is None: + self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma) + return self._MeSigma + + @property + def MfRhoI(self): + """ + Inverse of :code:`MfRho` + """ + if getattr(self, '_MfRhoI', None) is None: + self._MfRhoI = self.mesh.getFaceInnerProduct(self.rho, invMat=True) + return self._MfRhoI + + def MfRhoIDeriv(self,u): + """ + Derivative of :code:`MfRhoI` with respect to the model. + """ + + dMfRhoI_dI = -self.MfRhoI**2 + dMf_drho = self.mesh.getFaceInnerProductDeriv(self.rho)(u) + drho_dlogrho = Utils.sdiag(self.rho) + return dMfRhoI_dI * ( dMf_drho * ( drho_dlogrho)) + + # TODO: This should take a vector + def MeSigmaDeriv(self, u): + """ + Derivative of MeSigma with respect to the model + """ + dsigma_dlogsigma = Utils.sdiag(self.sigma) + return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) * dsigma_dlogsigma + +class Problem3D_CC(BaseSIPProblem): + + _solutionType = 'phiSolution' + _formulation = 'HJ' # CC potentials means J is on faces + fieldsPair = Fields_CC + + def __init__(self, mesh, **kwargs): + BaseSIPProblem.__init__(self, mesh, **kwargs) + self.setBC() + + def getA(self): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = D MfRhoI D^\\top V + + """ + + D = self.Div + G = self.Grad + # TODO: this won't work for full anisotropy + MfRhoI = self.MfRhoI + A = D * MfRhoI * G + + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return V.T * A + return A + + def getADeriv(self, u, v, adjoint= False): + + D = self.Div + G = self.Grad + MfRhoIDeriv = self.MfRhoIDeriv + + if adjoint: + # if self._makeASymmetric is True: + # v = V * v + return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) ) + return D * (MfRhoIDeriv( G * u ) * v) + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return self.Vol.T * RHS + + return RHS + + def getRHSDeriv(self, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() + + def setBC(self): + if self.mesh.dim==3: + fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + gBFzm = self.mesh.gridFz[fzm,:] + gBFzp = self.mesh.gridFz[fzp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + alpha_zm, alpha_zp = temp_zm*0., temp_zp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + beta_zm, beta_zp = temp_zm, temp_zp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + gamma_zm, gamma_zp = temp_zm*0., temp_zp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp] + + elif self.mesh.dim==2: + + fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd + gBFxm = self.mesh.gridFx[fxm,:] + gBFxp = self.mesh.gridFx[fxp,:] + gBFym = self.mesh.gridFy[fym,:] + gBFyp = self.mesh.gridFy[fyp,:] + + # Setup Mixed B.C (alpha, beta, gamma) + temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0]) + temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1]) + + alpha_xm, alpha_xp = temp_xm*0., temp_xp*0. + alpha_ym, alpha_yp = temp_ym*0., temp_yp*0. + + beta_xm, beta_xp = temp_xm, temp_xp + beta_ym, beta_yp = temp_ym, temp_yp + + gamma_xm, gamma_xp = temp_xm*0., temp_xp*0. + gamma_ym, gamma_yp = temp_ym*0., temp_yp*0. + + alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp] + beta = [beta_xm, beta_xp, beta_ym, beta_yp] + gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp] + + x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma) + V = self.Vol + self.Div = V * self.mesh.faceDiv + P_BC, B = self.mesh.getBCProjWF_simple() + M = B*self.mesh.aveCC2F + self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M + + +class Problem3D_N(BaseSIPProblem): + + _solutionType = 'phiSolution' + _formulation = 'EB' # N potentials means B is on faces + fieldsPair = Fields_N + + def __init__(self, mesh, **kwargs): + BaseSIPProblem.__init__(self, mesh, **kwargs) + + def getA(self): + """ + + Make the A matrix for the cell centered DC resistivity problem + + A = D MfRhoI D^\\top V + + """ + + # TODO: this won't work for full anisotropy + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + A = Grad.T * MeSigma * Grad + + # Handling Null space of A + A[0,0] = A[0,0] + 1. + + return A + + def getADeriv(self, u, v, adjoint=False): + """ + + Product of the derivative of our system matrix with respect to the model and a vector + + """ + MeSigma = self.MeSigma + Grad = self.mesh.nodalGrad + if not adjoint: + return Grad.T*(self.MeSigmaDeriv(Grad*u)*v) + elif adjoint: + return self.MeSigmaDeriv(Grad*u).T * (Grad*v) + + + def getRHS(self): + """ + RHS for the DC problem + + q + """ + + RHS = self.getSourceTerm() + return RHS + + def getRHSDeriv(self, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + """ + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() + +if __name__ == '__main__': + + + cs = 12.5 + hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)] + hz = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + sigma = np.ones(mesh.nC) + prob = BaseSIPProblem(mesh, sigma=sigma) + + diff --git a/SimPEG/EM/Static/SIP/RxSIP.py b/SimPEG/EM/Static/SIP/RxSIP.py new file mode 100644 index 00000000..f4ff1a2f --- /dev/null +++ b/SimPEG/EM/Static/SIP/RxSIP.py @@ -0,0 +1,82 @@ +import SimPEG +import numpy as np +from SimPEG.Utils import Zero, closestPoints + +class BaseRx(SimPEG.Survey.BaseTimeRx): + locs = None + rxType = None + + knownRxTypes = { + 'phi':['phi',None], + 'ex':['e','x'], + 'ey':['e','y'], + 'ez':['e','z'], + 'jx':['j','x'], + 'jy':['j','y'], + 'jz':['j','z'], + } + + def __init__(self, locs, times, rxType, **kwargs): + SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType, **kwargs) + + @property + def projField(self): + """Field Type projection (e.g. e b ...)""" + return self.knownRxTypes[self.rxType][0] + + def projGLoc(self, f): + """Grid Location projection (e.g. Ex Fy ...)""" + comp = self.knownRxTypes[self.rxType][1] + if comp is not None: + return f._GLoc(self.rxType) + comp + return f._GLoc(self.rxType) + + def getTimeP(self, timesall): + """ + Returns the time projection matrix. + + .. note:: + + This is not stored in memory, but is created on demand. + """ + time_inds = np.in1d(timesall, self.times) + return time_inds + + def evalDeriv(self, src, mesh, f, v, adjoint=False): + P = self.getP(mesh, self.projGLoc(f)) + if not adjoint: + return P*v + elif adjoint: + return P.T*v + + +# DC.Rx.Dipole(locs) +class Dipole(BaseRx): + + def __init__(self, locsM, locsN, times, rxType = 'phi', **kwargs): + assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size' + locs = [locsM, locsN] + # We may not need this ... + BaseRx.__init__(self, locs, times, rxType) + + @property + def nD(self): + """Number of data in the receiver.""" + return self.locs[0].shape[0] * len(self.times) + + # Not sure why ... + # return int(self.locs[0].size / 2) + + + def getP(self, mesh, Gloc): + if mesh in self._Ps: + return self._Ps[mesh] + + P0 = mesh.getInterpolationMat(self.locs[0], Gloc) + P1 = mesh.getInterpolationMat(self.locs[1], Gloc) + P = P0 - P1 + + if self.storeProjections: + self._Ps[mesh] = P + + return P diff --git a/SimPEG/EM/Static/SIP/SurveySIP.py b/SimPEG/EM/Static/SIP/SurveySIP.py new file mode 100644 index 00000000..de375794 --- /dev/null +++ b/SimPEG/EM/Static/SIP/SurveySIP.py @@ -0,0 +1,32 @@ +import SimPEG +from SimPEG.EM.Base import BaseEMSurvey +from SimPEG import np, sp, Survey +from SimPEG.Utils import Zero, Identity +from SimPEG.EM.Static.DC.SrcDC import BaseSrc +from SimPEG.EM.Static.SIP.RxSIP import BaseRx + +class Survey(BaseEMSurvey): + rxPair = BaseRx + srcPair = BaseSrc + times = None + + def __init__(self, srcList, **kwargs): + self.srcList = srcList + BaseEMSurvey.__init__(self, srcList, **kwargs) + self.getUniqueTimes() + + def getUniqueTimes(self): + time_rx = [] + for src in self.srcList: + for rx in src.rxList: + time_rx.append(rx.times) + self.times = np.unique(np.hstack(time_rx)) + + def dpred(self, m, f=None): + """ + Predicted data. + + .. math:: + d_\\text{pred} = Pf(m) + """ + return self.prob.Jvec(m, m, f=f) diff --git a/SimPEG/EM/Static/SIP/__init__.py b/SimPEG/EM/Static/SIP/__init__.py new file mode 100644 index 00000000..f81ceb3c --- /dev/null +++ b/SimPEG/EM/Static/SIP/__init__.py @@ -0,0 +1,2 @@ +from ProblemSIP import Problem3D_CC, Problem3D_N +from SurveySIP import Survey diff --git a/SimPEG/EM/Static/__init__.py b/SimPEG/EM/Static/__init__.py index c9b4dc0d..f6e0e757 100644 --- a/SimPEG/EM/Static/__init__.py +++ b/SimPEG/EM/Static/__init__.py @@ -1,2 +1,3 @@ import DC import IP +import SIP From e8bd78f63d8ac84e749e15471e69e40a8174ccd8 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 19 May 2016 02:09:48 +0900 Subject: [PATCH 2/6] Working Spectral IP: - Fwd - Jvec - Jtvec --- SimPEG/EM/Static/SIP/ProblemSIP.py | 146 +++++++++++---------- SimPEG/EM/Static/SIP/RxSIP.py | 8 +- SimPEG/EM/Static/SIP/SrcSIP.py | 64 ++++++++++ SimPEG/EM/Static/SIP/SurveySIP.py | 76 ++++++++++- SimPEG/EM/Static/SIP/__init__.py | 4 +- tests/em/static/test_SIP_jvecjtvecadj.py | 154 +++++++++++++++++++++++ 6 files changed, 382 insertions(+), 70 deletions(-) create mode 100644 SimPEG/EM/Static/SIP/SrcSIP.py create mode 100644 tests/em/static/test_SIP_jvecjtvecadj.py diff --git a/SimPEG/EM/Static/SIP/ProblemSIP.py b/SimPEG/EM/Static/SIP/ProblemSIP.py index d01ae706..e4ab33d8 100644 --- a/SimPEG/EM/Static/SIP/ProblemSIP.py +++ b/SimPEG/EM/Static/SIP/ProblemSIP.py @@ -5,7 +5,7 @@ from SimPEG.Utils import sdiag import numpy as np from SimPEG.Utils import Zero from SimPEG.EM.Static.DC import getxBCyBC_CC -from SurveySIP import Survey +from SurveySIP import Survey, Data class ColeColePropMap(Maps.PropMap): """ @@ -22,6 +22,7 @@ class BaseSIPProblem(BaseEMProblem): surveyPair = Survey fieldsPair = Fields + dataPair = Data PropMap = ColeColePropMap Ainv = None sigma = None @@ -29,15 +30,17 @@ class BaseSIPProblem(BaseEMProblem): f = None Ainv = None - def DebyeTime(t): + def DebyeTime(self, t): peta = self.curModel.eta*np.exp(-self.curModel.taui*t) return peta - def EtaDeriv(t): - return np.exp(-self.curModel.taui*t) + def EtaDeriv(self, t, v): + v = np.array(v, dtype=float) + return np.exp(-self.curModel.taui*t) * (self.curModel.etaDeriv*v) - def TauiDeriv(t): - return -self.curModel.eta*t*np.exp(-self.curModel.taui*t) + def TauiDeriv(self, t, v): + v = np.array(v, dtype=float) + return -self.curModel.eta*t*np.exp(-self.curModel.taui*t) * (self.curModel.tauiDeriv*v) def fields(self, m): self.curModel = m @@ -63,7 +66,8 @@ class BaseSIPProblem(BaseEMProblem): JvAll = [] for tind in range(len(self.survey.times)): #Pseudo-chareability - v = DebyeTime(self.survey.times[tind]) + t = self.survey.times[tind] + v = self.DebyeTime(t) for src in self.survey.srcList: u_src = f[src, self._solutionType] # solution vector dA_dm_v = self.getADeriv(u_src, v) @@ -74,76 +78,88 @@ class BaseSIPProblem(BaseEMProblem): if timeindex[tind]: df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) - Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) - JvAll.append(Utils.mkvc(Jv)) + Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v) + # Conductivity (d u / d log sigma) if self._formulation is 'EB': - return -np.hstack(JvAll) - # Conductivity (d u / d log rho) + return -Utils.mkvc(Jv) + # Resistivity (d u / d log rho) if self._formulation is 'HJ': - return np.hstack(JvAll) + return Utils.mkvc(Jv) - # def Jvec(self, m, v, f=None): + def Jvec(self, m, v, f=None): - # if f is None: - # f = self.fields(m) + if f is None: + f = self.fields(m) - # self.curModel = m + self.curModel = m + Jv = self.dataPair(self.survey) #same size as the data + # A = self.getA() + JvAll = [] + #Assume only eta and tau (eta first then tau) + # v = [2*Mx1] + v = v.reshape((int(v.size/2), 2), order='F') - # Jv = self.dataPair(self.survey) #same size as the data - # x1 = - # x2 = - # # A = self.getA() - # for src in self.survey.srcList: - # u_src = f[src, self._solutionType] # solution vector + for tind in range(len(self.survey.times)): + t = self.survey.times[tind] + v0 = self.EtaDeriv(t, v[:,0]) + v1 = self.TauiDeriv(t, v[:,1]) + for src in self.survey.srcList: + u_src = f[src, self._solutionType] # solution vector + dA_dm_v0 = self.getADeriv(u_src, v0) + dRHS_dm_v0 = self.getRHSDeriv(src, v0) + du_dm_v0 = self.Ainv * ( - dA_dm_v0 + dRHS_dm_v0 ) + dA_dm_v1 = self.getADeriv(u_src, v1) + dRHS_dm_v1 = self.getRHSDeriv(src, v1) + du_dm_v1 = self.Ainv * ( - dA_dm_v1 + dRHS_dm_v1 ) + for rx in src.rxList: + timeindex = rx.getTimeP(self.survey.times) + if timeindex[tind]: + df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dm_v0 = df_dmFun(src, du_dm_v0, v0, adjoint=False) + df_dm_v1 = df_dmFun(src, du_dm_v1, v1, adjoint=False) + Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v0) + Jv[src, rx, t] += rx.evalDeriv(src, self.mesh, f, df_dm_v1) + # Conductivity (d u / d log sigma) + if self._formulation is 'EB': + return -Jv.tovec() + # Resistivity (d u / d log rho) + if self._formulation is 'HJ': + return Jv.tovec() - # dA_dm_v = self.getADeriv(u_src, v) - # dRHS_dm_v = self.getRHSDeriv(src, v) - # du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v ) + def Jtvec(self, m, v, f=None): + if f is None: + f = self.fields(m) - # for rx in src.rxList: + self.curModel = m - # for tind in range(len(self.survey.times)): - # df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) - # df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) - # Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) - # # Conductivity (d u / d log sigma) - # if self._formulation is 'EB': - # return -Utils.mkvc(Jv) - # # Conductivity (d u / d log rho) - # if self._formulation is 'HJ': - # return Utils.mkvc(Jv) + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) - # def Jtvec(self, m, v, f=None): - # if f is None: - # f = self.fields(m) + Jtv= np.zeros(m.size) + for tind in range(len(self.survey.times)): + t = self.survey.times[tind] + for src in self.survey.srcList: + u_src = f[src, self._solutionType] + for rx in src.rxList: + timeindex = rx.getTimeP(self.survey.times) + if timeindex[tind]: + PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx, t], adjoint=True) # wrt f, need possibility wrt m + df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) + ATinvdf_duT = self.Ainv * df_duT + dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) + du_dmT = -dA_dmT + dRHS_dmT + Jtv += np.r_[self.EtaDeriv(self.survey.times[tind], du_dmT), self.TauiDeriv(self.survey.times[tind], du_dmT)] - # self.curModel = m - - # # Ensure v is a data object. - # if not isinstance(v, self.dataPair): - # v = self.dataPair(self.survey, v) - - # Jtv = np.zeros(m.size) - # # AT = self.getA() - - # for src in self.survey.srcList: - # u_src = f[src, self._solutionType] - # for rx in src.rxList: - # PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m - # df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) - # df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) - # ATinvdf_duT = self.Ainv * df_duT - # dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) - # dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) - # du_dmT = -dA_dmT + dRHS_dmT - # Jtv += df_dmT + du_dmT - # # Conductivity ((d u / d log sigma).T) - # if self._formulation is 'EB': - # return -Utils.mkvc(Jtv) - # # Conductivity ((d u / d log rho).T) - # if self._formulation is 'HJ': - # return Utils.mkvc(Jtv) + # Conductivity ((d u / d log sigma).T) + if self._formulation is 'EB': + return -Jtv + # Conductivity ((d u / d log rho).T) + if self._formulation is 'HJ': + return Jtv def getSourceTerm(self): """ diff --git a/SimPEG/EM/Static/SIP/RxSIP.py b/SimPEG/EM/Static/SIP/RxSIP.py index f4ff1a2f..b8144d0f 100644 --- a/SimPEG/EM/Static/SIP/RxSIP.py +++ b/SimPEG/EM/Static/SIP/RxSIP.py @@ -62,7 +62,13 @@ class Dipole(BaseRx): @property def nD(self): """Number of data in the receiver.""" - return self.locs[0].shape[0] * len(self.times) + # return self.locs[0].shape[0] * len(self.times) + return self.locs[0].shape[0] + + @property + def nRx(self): + """Number of data in the receiver.""" + return self.locs[0].shape[0] # Not sure why ... # return int(self.locs[0].size / 2) diff --git a/SimPEG/EM/Static/SIP/SrcSIP.py b/SimPEG/EM/Static/SIP/SrcSIP.py new file mode 100644 index 00000000..b1f0a452 --- /dev/null +++ b/SimPEG/EM/Static/SIP/SrcSIP.py @@ -0,0 +1,64 @@ +import SimPEG +# from SimPEG.EM.Base import BaseEMSurvey +from SimPEG.Utils import Zero, closestPoints, mkvc +import numpy as np + +class BaseSrc(SimPEG.Survey.BaseSrc): + + current = 1.0 + loc = None + + def __init__(self, rxList, **kwargs): + SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs) + + def eval(self, prob): + raise NotImplementedError + + def evalDeriv(self, prob): + return Zero() + + @property + def nD(self): + """Number of data""" + return self.vnD.sum() + + @property + def vnD(self): + """Vector number of data""" + return np.array([rx.nD*len(rx.times) for rx in self.rxList]) + + + +class Dipole(BaseSrc): + + def __init__(self, rxList, locA, locB, **kwargs): + assert locA.shape == locB.shape, 'Shape of locA and locB should be the same' + self.loc = [locA, locB] + BaseSrc.__init__(self, rxList, **kwargs) + + def eval(self, prob): + if prob._formulation == 'HJ': + inds = closestPoints(prob.mesh, self.loc, gridLoc='CC') + q = np.zeros(prob.mesh.nC) + q[inds] = self.current * np.r_[1., -1.] + elif prob._formulation == 'EB': + qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense() + qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense() + q = self.current * mkvc(qa+qb) + return q + +class Pole(BaseSrc): + + def __init__(self, rxList, loc, **kwargs): + BaseSrc.__init__(self, rxList, loc=loc, **kwargs) + + def eval(self, prob): + if prob._formulation == 'HJ': + inds = closestPoints(prob.mesh, self.loc) + q = np.zeros(prob.mesh.nC) + q[inds] = self.current * np.r_[1.] + elif prob._formulation == 'EB': + q = prob.mesh.getInterpolationMat(self.loc, locType='N').todense() + q = self.current * mkvc(q) + return q + diff --git a/SimPEG/EM/Static/SIP/SurveySIP.py b/SimPEG/EM/Static/SIP/SurveySIP.py index de375794..9362a35c 100644 --- a/SimPEG/EM/Static/SIP/SurveySIP.py +++ b/SimPEG/EM/Static/SIP/SurveySIP.py @@ -1,9 +1,11 @@ import SimPEG from SimPEG.EM.Base import BaseEMSurvey -from SimPEG import np, sp, Survey +from SimPEG import np, sp, Survey, Utils from SimPEG.Utils import Zero, Identity -from SimPEG.EM.Static.DC.SrcDC import BaseSrc +from SimPEG.EM.Static.SIP.SrcSIP import BaseSrc from SimPEG.EM.Static.SIP.RxSIP import BaseRx +import uuid + class Survey(BaseEMSurvey): rxPair = BaseRx @@ -29,4 +31,72 @@ class Survey(BaseEMSurvey): .. math:: d_\\text{pred} = Pf(m) """ - return self.prob.Jvec(m, m, f=f) + return self.prob.forward(m, f=f) + + +class Data(SimPEG.Survey.Data): + """Fancy data storage by Src and Rx""" + + def __init__(self, survey, v=None): + self.uid = str(uuid.uuid4()) + self.survey = survey + self._dataDict = {} + for src in self.survey.srcList: + self._dataDict[src] = {} + for rx in src.rxList: + self._dataDict[src][rx] = {} + + if v is not None: + self.fromvec(v) + + def _ensureCorrectKey(self, key): + if type(key) is tuple: + if len(key) is not 3: + raise KeyError('Key must be [Src, Rx, tInd]') + if key[0] not in self.survey.srcList: + raise KeyError('Src Key must be a source in the survey.') + if key[1] not in key[0].rxList: + raise KeyError('Rx Key must be a receiver for the source.') + return key + elif isinstance(key, self.survey.srcPair): + if key not in self.survey.srcList: + raise KeyError('Key must be a source in the survey.') + return key, None, None + else: + raise KeyError('Key must be [Src] or [Src,Rx] or [Src, Rx, tInd]') + + def __setitem__(self, key, value): + src, rx, t = self._ensureCorrectKey(key) + assert rx is not None, 'set data using [Src, Rx]' + assert isinstance(value, np.ndarray), 'value must by ndarray' + assert value.size == rx.nD, "value must have the same number of data as the source." + self._dataDict[src][rx][t] = Utils.mkvc(value) + + def __getitem__(self, key): + src, rx, t = self._ensureCorrectKey(key) + if rx is not None: + if rx not in self._dataDict[src]: + raise Exception('Data for receiver has not yet been set.') + return self._dataDict[src][rx][t] + + return np.concatenate([self[src,rx, t] for rx in src.rxList]) + + def tovec(self): + val = [] + for src in self.survey.srcList: + for rx in src.rxList: + for t in rx.times: + val.append(self[src, rx, t]) + return np.concatenate(val) + + + def fromvec(self, v): + v = Utils.mkvc(v) + assert v.size == self.survey.nD, 'v must have the correct number of data.' + indBot, indTop = 0, 0 + for src in self.survey.srcList: + for rx in src.rxList: + for t in rx.times: + indTop += rx.nRx + self[src, rx, t] = v[indBot:indTop] + indBot += rx.nRx diff --git a/SimPEG/EM/Static/SIP/__init__.py b/SimPEG/EM/Static/SIP/__init__.py index f81ceb3c..94fdb591 100644 --- a/SimPEG/EM/Static/SIP/__init__.py +++ b/SimPEG/EM/Static/SIP/__init__.py @@ -1,2 +1,4 @@ from ProblemSIP import Problem3D_CC, Problem3D_N -from SurveySIP import Survey +from SurveySIP import Survey, Data +import SrcSIP as Src #Pole +import RxSIP as Rx diff --git a/tests/em/static/test_SIP_jvecjtvecadj.py b/tests/em/static/test_SIP_jvecjtvecadj.py new file mode 100644 index 00000000..6e6f071d --- /dev/null +++ b/tests/em/static/test_SIP_jvecjtvecadj.py @@ -0,0 +1,154 @@ +import unittest +from SimPEG import * +import SimPEG +from SimPEG import Mesh, Utils, EM, Maps, np, Survey +from SimPEG.EM.Static import SIP, DC, IP +from pymatsolver import MumpsSolver + + +class IPProblemTestsCC(unittest.TestCase): + + def setUp(self): + + cs = 25. + hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hz = [(cs,0, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC) + blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC) + sigma = np.ones(mesh.nC)*1e-2 + eta = np.zeros(mesh.nC) + tau = np.ones_like(sigma)*1. + eta[blkind0] = 0.1 + eta[blkind1] = 0.1 + tau[blkind0] = 0.1 + tau[blkind1] = 0.01 + + x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] + Aloc = np.r_[-200., 0., 0.] + Bloc = np.r_[200., 0., 0.] + M = Utils.ndgrid(x-25.,y, np.r_[0.]) + N = Utils.ndgrid(x+25.,y, np.r_[0.]) + + times = np.arange(10)*1e-3 + 1e-3 + rx = SIP.Rx.Dipole(M, N, times) + src = SIP.Src.Dipole([rx], Aloc, Bloc) + survey = SIP.Survey([src]) + colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))] + problem = SIP.Problem3D_CC(mesh, rho=1./sigma, mapping=colemap) + problem.Solver = MumpsSolver + problem.pair(survey) + mSynth = np.r_[eta, 1./tau] + survey.makeSyntheticData(mSynth) + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC*2) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-10 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + +class IPProblemTestsN(unittest.TestCase): + + def setUp(self): + + cs = 25. + hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hz = [(cs,0, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN") + blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC) + blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC) + sigma = np.ones(mesh.nC)*1e-2 + eta = np.zeros(mesh.nC) + tau = np.ones_like(sigma)*1. + eta[blkind0] = 0.1 + eta[blkind1] = 0.1 + tau[blkind0] = 0.1 + tau[blkind1] = 0.01 + + x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] + Aloc = np.r_[-200., 0., 0.] + Bloc = np.r_[200., 0., 0.] + M = Utils.ndgrid(x-25.,y, np.r_[0.]) + N = Utils.ndgrid(x+25.,y, np.r_[0.]) + + times = np.arange(10)*1e-3 + 1e-3 + rx = SIP.Rx.Dipole(M, N, times) + src = SIP.Src.Dipole([rx], Aloc, Bloc) + survey = SIP.Survey([src]) + colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))] + problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap) + problem.Solver = MumpsSolver + problem.pair(survey) + mSynth = np.r_[eta, 1./tau] + survey.makeSyntheticData(mSynth) + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + reg = Regularization.Tikhonov(mesh) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC*2) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-8 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() From c4c97ae054602376cf2d5df596b3d5729a723709 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Fri, 20 May 2016 00:04:20 -0700 Subject: [PATCH 3/6] Ad MultiRegularization for inverting multiple parameters --- SimPEG/EM/Static/SIP/Regularization.py | 105 +++++++++++++++++++++++++ SimPEG/EM/Static/SIP/__init__.py | 1 + 2 files changed, 106 insertions(+) create mode 100644 SimPEG/EM/Static/SIP/Regularization.py diff --git a/SimPEG/EM/Static/SIP/Regularization.py b/SimPEG/EM/Static/SIP/Regularization.py new file mode 100644 index 00000000..4a067929 --- /dev/null +++ b/SimPEG/EM/Static/SIP/Regularization.py @@ -0,0 +1,105 @@ +from SimPEG import Utils, Maps, Mesh, sp, np +from SimPEG.Regularization import BaseRegularization, Simple + +class MultiRegularization(Simple): + """ + **MultiRegularization Class** + + This is used to regularize the model space + having multiple models [m1, m2, m3, ...] :: + + reg = Regularization(mesh) + + """ + nModels = None # Number of models + ratios = None + crossgrad = False + + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): + BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) + if self.nModels == None: + raise Exception("Put nModels as a initial input!") + if self.ratios == None: + self.ratios = [1. for imodel in range(self.nModels)] + + @property + def Wsmall(self): + """Regularization matrix Wsmall""" + if getattr(self,'_Wsmall', None) is None: + vecs = [] + for imodel in range(self.nModels): + vecs.append((self.regmesh.vol*self.alpha_s*self.wght*self.ratios[imodel])**0.5) + self._Wsmall = Utils.sdiag(np.hstack(vecs)) + return self._Wsmall + + @property + def Wx(self): + """Regularization matrix Wx""" + if getattr(self, '_Wx', None) is None: + mats = [] + for imodel in range(self.nModels): + mats.append(Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*self.ratios[imodel]*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil) + self._Wx = sp.block_diag(mats) + return self._Wx + + @property + def Wy(self): + """Regularization matrix Wy""" + if getattr(self, '_Wy', None) is None: + mats = [] + for imodel in range(self.nModels): + mats.append(Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol*self.alpha_y*self.ratios[imodel]*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil) + self._Wy = sp.block_diag(mats) + return self._Wy + + @property + def Wz(self): + """Regularization matrix Wz""" + if getattr(self, '_Wz', None) is None: + mats = [] + for imodel in range(self.nModels): + mats.append(Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*self.ratios[imodel]*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil) + self._Wz = sp.block_diag(mats) + return self._Wz + + @property + def Wsmooth(self): + """Full smoothness regularization matrix W""" + if getattr(self, '_Wsmooth', None) is None: + wlist = (self.Wx,) + if self.regmesh.dim > 1: + wlist += (self.Wy,) + if self.regmesh.dim > 2: + wlist += (self.Wz,) + self._Wsmooth = sp.vstack(wlist) + return self._Wsmooth + + @property + def W(self): + """Full regularization matrix W""" + if getattr(self, '_W', None) is None: + wlist = (self.Wsmall, self.Wsmooth) + self._W = sp.vstack(wlist) + return self._W + + @Utils.timeIt + def _evalSmall(self, m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmooth(self, m): + if self.mrefInSmooth == True: + r = self.Wsmooth * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wsmooth * ( self.mapping * m) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalCross(self, m): + if self.crossgrad == False: + return 0. + elif self.crossgrad == True: + r = self.Wcross * ( self.mapping * m) + return 0.5 * r.dot(r) + diff --git a/SimPEG/EM/Static/SIP/__init__.py b/SimPEG/EM/Static/SIP/__init__.py index 94fdb591..1de46fcf 100644 --- a/SimPEG/EM/Static/SIP/__init__.py +++ b/SimPEG/EM/Static/SIP/__init__.py @@ -2,3 +2,4 @@ from ProblemSIP import Problem3D_CC, Problem3D_N from SurveySIP import Survey, Data import SrcSIP as Src #Pole import RxSIP as Rx +from Regularization import MultiRegularization From 0179631fe31db5bfa2b6800210cdc973f8bcb683 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Fri, 20 May 2016 01:57:40 -0700 Subject: [PATCH 4/6] Starting Cross gradient ... --- SimPEG/EM/Static/SIP/Regularization.py | 47 +++++++++++++++++++------- 1 file changed, 35 insertions(+), 12 deletions(-) diff --git a/SimPEG/EM/Static/SIP/Regularization.py b/SimPEG/EM/Static/SIP/Regularization.py index 4a067929..a327ce4d 100644 --- a/SimPEG/EM/Static/SIP/Regularization.py +++ b/SimPEG/EM/Static/SIP/Regularization.py @@ -12,8 +12,12 @@ class MultiRegularization(Simple): """ nModels = None # Number of models - ratios = None - crossgrad = False + ratios = None # Ratio for different models + crossgrad = False # Use cross gradient or not + betacross = 1. + wx = [] + wy = [] + wz = [] def __init__(self, mesh, mapping=None, indActive=None, **kwargs): BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) @@ -38,7 +42,8 @@ class MultiRegularization(Simple): if getattr(self, '_Wx', None) is None: mats = [] for imodel in range(self.nModels): - mats.append(Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*self.ratios[imodel]*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil) + self.wx.append(Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*self.ratios[imodel]*(self.regmesh.aveCC2Fx*self.wght))**0.5)) + mats.append(self.wx[imodel]*self.regmesh.cellDiffxStencil) self._Wx = sp.block_diag(mats) return self._Wx @@ -48,7 +53,8 @@ class MultiRegularization(Simple): if getattr(self, '_Wy', None) is None: mats = [] for imodel in range(self.nModels): - mats.append(Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol*self.alpha_y*self.ratios[imodel]*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil) + self.wy.append(Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol*self.alpha_y*self.ratios[imodel]*(self.regmesh.aveCC2Fy*self.wght))**0.5)) + mats.append(self.wy[imodel]*self.regmesh.cellDiffyStencil) self._Wy = sp.block_diag(mats) return self._Wy @@ -58,7 +64,8 @@ class MultiRegularization(Simple): if getattr(self, '_Wz', None) is None: mats = [] for imodel in range(self.nModels): - mats.append(Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*self.ratios[imodel]*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil) + self.wz.append(Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*self.ratios[imodel]*(self.regmesh.aveCC2Fz*self.wght))**0.5)) + mats.append(self.wz[imodel]*self.regmesh.cellDiffzStencil) self._Wz = sp.block_diag(mats) return self._Wz @@ -95,11 +102,27 @@ class MultiRegularization(Simple): r = self.Wsmooth * ( self.mapping * m) return 0.5 * r.dot(r) - @Utils.timeIt - def _evalCross(self, m): - if self.crossgrad == False: - return 0. - elif self.crossgrad == True: - r = self.Wcross * ( self.mapping * m) - return 0.5 * r.dot(r) + # TODO: Implement Cross Gradients.. + # @Utils.timeIt + # def _evalCross(self, m): + # if self.crossgrad == False: + # return 0. + # elif self.crossgrad == True: + # M = (self.mapping * m).reshape((self.regmesh.nC, self.nModels), order="F") + + # for imodel in range(self.nModels): + # ux.append(self.regmesh.aveFx2CC*self.regmesh.wx[imodel]*M[:,imodel]) + # uy.append(self.regmesh.aveFy2CC*self.regmesh.wy[imodel]*M[:,imodel]) + # uz.append(self.regmesh.aveFz2CC*self.regmesh.wz[imodel]*M[:,imodel]) + + # ax, ay, az = ux[0], uy[0], uz[0] + # for imodel in range(1,self.nModels): + # bx, by, bz = ux[imodel], uy[imodel], uz[imodel] + # cx = ay*bz - az*by + # cy = az*bx - ax*bz + # cz = ax*by - ay*bx + # ax, ay, az = cx.copy(), cy.copy(), cz.copy() + # r = np.r_[ax, ay, az]*np.sqrt(self.betacross) + + # return 0.5 * r.dot(r) From f20fcb4504ef3f229d951982d6bc39a09c58da6a Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Tue, 24 May 2016 08:45:12 -0700 Subject: [PATCH 5/6] Minor type error based upon numpy version Cross gradient? --- SimPEG/EM/Static/DC/ProblemDC.py | 2 +- SimPEG/EM/Static/IP/ProblemIP.py | 2 +- SimPEG/EM/Static/SIP/Regularization.py | 116 ++++++++++++++++++++----- 3 files changed, 98 insertions(+), 22 deletions(-) diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index 0dd4f259..51b528a5 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -77,7 +77,7 @@ class BaseDCProblem(BaseEMProblem): dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT - Jtv += df_dmT + du_dmT + Jtv += (df_dmT + du_dmT).astype(float) return Utils.mkvc(Jtv) diff --git a/SimPEG/EM/Static/IP/ProblemIP.py b/SimPEG/EM/Static/IP/ProblemIP.py index 4ff81990..a637516a 100644 --- a/SimPEG/EM/Static/IP/ProblemIP.py +++ b/SimPEG/EM/Static/IP/ProblemIP.py @@ -89,7 +89,7 @@ class BaseIPProblem(BaseEMProblem): dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT - Jtv += df_dmT + du_dmT + Jtv += (df_dmT + du_dmT).astype(float) # Conductivity ((d u / d log sigma).T) if self._formulation is 'EB': return -Utils.mkvc(Jtv) diff --git a/SimPEG/EM/Static/SIP/Regularization.py b/SimPEG/EM/Static/SIP/Regularization.py index a327ce4d..3c20e6fe 100644 --- a/SimPEG/EM/Static/SIP/Regularization.py +++ b/SimPEG/EM/Static/SIP/Regularization.py @@ -89,6 +89,11 @@ class MultiRegularization(Simple): self._W = sp.vstack(wlist) return self._W + + @Utils.timeIt + def eval(self, m): + return self._evalSmall(m) + self._evalSmooth(m) + @Utils.timeIt def _evalSmall(self, m): r = self.Wsmall * ( self.mapping * (m - self.mref) ) @@ -102,27 +107,98 @@ class MultiRegularization(Simple): r = self.Wsmooth * ( self.mapping * m) return 0.5 * r.dot(r) - # TODO: Implement Cross Gradients.. - # @Utils.timeIt - # def _evalCross(self, m): - # if self.crossgrad == False: - # return 0. - # elif self.crossgrad == True: - # M = (self.mapping * m).reshape((self.regmesh.nC, self.nModels), order="F") + def cross(a,b): + ax, ay, az = a[0], a[1], a[2] + bx, by, bz = b[0], b[1], b[2] + cx = ay*bz - az*by + cy = az*bx - ax*bz + cz = ax*by - ay*bx + return [cx, cy, cz] - # for imodel in range(self.nModels): - # ux.append(self.regmesh.aveFx2CC*self.regmesh.wx[imodel]*M[:,imodel]) - # uy.append(self.regmesh.aveFy2CC*self.regmesh.wy[imodel]*M[:,imodel]) - # uz.append(self.regmesh.aveFz2CC*self.regmesh.wz[imodel]*M[:,imodel]) + # TODO: Implement Cross Gradients.. + @Utils.timeIt + def _evalCross(self, m): + if self.crossgrad == False: + return 0. + elif self.crossgrad == True: + M = (self.mapping * m).reshape((self.regmesh.nC, self.nModels), order="F") + + ax = self.regmesh.aveFx2CC*self.regmesh.wx[0]*M[:,0] + ay = self.regmesh.aveFy2CC*self.regmesh.wy[0]*M[:,0] + az = self.regmesh.aveFz2CC*self.regmesh.wz[0]*M[:,0] + bx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1] + by = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1] + bz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1] + #ab + out_ab = cross([ax, ay, az], [bx, by, bz]) + r = np.r_[out_ab[0], out_ab[1], out_ab[2]]*np.sqrt(self.betacross) + + if self.nModels == 3: + cx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1] + cy = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1] + cz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1] + #ac + out_ac = cross([ax, ay, az], [cx, cy, cz]) + #bc + out_bc = cross([bx, by, bz], [cx, cy, cz]) + r = np.r_[r, np.hstack(out_ac)*np.sqrt(self.betacross), np.hstack(out_bc)*np.sqrt(self.betacross)] + + return 0.5 * r.dot(r) + + @Utils.timeIt + def evalDeriv(self, m): + """ + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} + + So the derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} + + """ + deriv = self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + if self.crossgrad==True: + deriv += self._evalCrossDeriv(m) + return deriv + + @Utils.timeIt + def _evalCrossDeriv(self,m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) + + @Utils.timeIt + def eval2Deriv(self, m, v=None): + """ + Second derivative + + :param numpy.array m: geophysical model + :param numpy.array v: vector to multiply + :rtype: scipy.sparse.csr_matrix or numpy.ndarray + :return: WtW or WtW*v + + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} + + So the second derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top W} + + """ + mD = self.mapping.deriv(m - self.mref) + if v is None: + return mD.T * self.W.T * self.W * mD + + return mD.T * ( self.W.T * ( self.W * ( mD * v) ) ) - # ax, ay, az = ux[0], uy[0], uz[0] - # for imodel in range(1,self.nModels): - # bx, by, bz = ux[imodel], uy[imodel], uz[imodel] - # cx = ay*bz - az*by - # cy = az*bx - ax*bz - # cz = ax*by - ay*bx - # ax, ay, az = cx.copy(), cy.copy(), cz.copy() - # r = np.r_[ax, ay, az]*np.sqrt(self.betacross) - # return 0.5 * r.dot(r) From 21d817d9a257138d90daa647bc67d9ec0ebbc58e Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Tue, 24 May 2016 21:53:18 -0700 Subject: [PATCH 6/6] fix bug for adjoint problem --- SimPEG/EM/Static/SIP/ProblemSIP.py | 17 ++++-- tests/em/static/test_SIP_jvecjtvecadj.py | 78 ++++++++++++++++++++++++ 2 files changed, 90 insertions(+), 5 deletions(-) diff --git a/SimPEG/EM/Static/SIP/ProblemSIP.py b/SimPEG/EM/Static/SIP/ProblemSIP.py index e4ab33d8..88b85a7b 100644 --- a/SimPEG/EM/Static/SIP/ProblemSIP.py +++ b/SimPEG/EM/Static/SIP/ProblemSIP.py @@ -34,13 +34,20 @@ class BaseSIPProblem(BaseEMProblem): peta = self.curModel.eta*np.exp(-self.curModel.taui*t) return peta - def EtaDeriv(self, t, v): + def EtaDeriv(self, t, v, adjoint=False): v = np.array(v, dtype=float) - return np.exp(-self.curModel.taui*t) * (self.curModel.etaDeriv*v) + if adjoint: + return self.curModel.etaDeriv.T * (np.exp(-self.curModel.taui*t)*v) + else: + return np.exp(-self.curModel.taui*t) * (self.curModel.etaDeriv*v) - def TauiDeriv(self, t, v): + + def TauiDeriv(self, t, v, adjoint=False): v = np.array(v, dtype=float) - return -self.curModel.eta*t*np.exp(-self.curModel.taui*t) * (self.curModel.tauiDeriv*v) + if adjoint: + return -self.curModel.tauiDeriv.T * (self.curModel.eta*t*np.exp(-self.curModel.taui*t)*v) + else: + return -self.curModel.eta*t*np.exp(-self.curModel.taui*t) * (self.curModel.tauiDeriv*v) def fields(self, m): self.curModel = m @@ -152,7 +159,7 @@ class BaseSIPProblem(BaseEMProblem): dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT - Jtv += np.r_[self.EtaDeriv(self.survey.times[tind], du_dmT), self.TauiDeriv(self.survey.times[tind], du_dmT)] + Jtv += np.r_[self.EtaDeriv(self.survey.times[tind], du_dmT, adjoint=True), self.TauiDeriv(self.survey.times[tind], du_dmT, adjoint=True)] # Conductivity ((d u / d log sigma).T) if self._formulation is 'EB': diff --git a/tests/em/static/test_SIP_jvecjtvecadj.py b/tests/em/static/test_SIP_jvecjtvecadj.py index 6e6f071d..7cdb6def 100644 --- a/tests/em/static/test_SIP_jvecjtvecadj.py +++ b/tests/em/static/test_SIP_jvecjtvecadj.py @@ -150,5 +150,83 @@ class IPProblemTestsN(unittest.TestCase): passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) self.assertTrue(passed) +class IPProblemTestsN_air(unittest.TestCase): + + def setUp(self): + + cs = 25. + hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hz = [(cs,0, -1.3),(cs,20),(cs,0, 1.3)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCC") + blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC) + blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC) + sigma = np.ones(mesh.nC)*1e-2 + airind = mesh.gridCC[:,2]>0. + sigma[airind] = 1e-8 + eta = np.zeros(mesh.nC) + tau = np.ones_like(sigma)*1. + eta[blkind0] = 0.1 + eta[blkind1] = 0.1 + tau[blkind0] = 0.1 + tau[blkind1] = 0.01 + + actmapeta = Maps.InjectActiveCells(mesh, ~airind, 0.) + actmaptau = Maps.InjectActiveCells(mesh, ~airind, 1.) + + x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] + Aloc = np.r_[-200., 0., 0.] + Bloc = np.r_[200., 0., 0.] + M = Utils.ndgrid(x-25.,y, np.r_[0.]) + N = Utils.ndgrid(x+25.,y, np.r_[0.]) + + times = np.arange(10)*1e-3 + 1e-3 + rx = SIP.Rx.Dipole(M, N, times) + src = SIP.Src.Dipole([rx], Aloc, Bloc) + survey = SIP.Survey([src]) + colemap = [("eta", Maps.IdentityMap(mesh)*actmapeta), ("taui", Maps.IdentityMap(mesh)*actmaptau)] + problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap) + problem.Solver = MumpsSolver + problem.pair(survey) + mSynth = np.r_[eta[~airind], 1./tau[~airind]] + survey.makeSyntheticData(mSynth) + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + regmap = Maps.IdentityMap(nP=int(mSynth[~airind].size*2)) + reg = SIP.MultiRegularization(mesh, mapping=regmap, nModels=2, indActive=~airind) + opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4) + inv = Inversion.BaseInversion(invProb) + + self.inv = inv + self.reg = reg + self.p = problem + self.mesh = mesh + self.m0 = mSynth + self.survey = survey + self.dmis = dmis + + def test_misfit(self): + derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + + def test_adjoint(self): + # Adjoint Test + u = np.random.rand(self.mesh.nC*self.survey.nSrc) + v = np.random.rand(self.mesh.nC) + w = np.random.rand(self.survey.dobs.shape[0]) + wtJv = w.dot(self.p.Jvec(self.m0, v)) + vtJtw = v.dot(self.p.Jtvec(self.m0, w)) + passed = np.abs(wtJv - vtJtw) < 1e-8 + print 'Adjoint Test', np.abs(wtJv - vtJtw), passed + self.assertTrue(passed) + + def test_dataObj(self): + derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)] + passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) + self.assertTrue(passed) + if __name__ == '__main__': unittest.main()