From e8bd78f63d8ac84e749e15471e69e40a8174ccd8 Mon Sep 17 00:00:00 2001 From: seogi_macbook Date: Thu, 19 May 2016 02:09:48 +0900 Subject: [PATCH] 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()