From c77a0279c8b4f90ed27a3f9e2bc8e020340383a2 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 29 Apr 2014 11:36:35 -0700 Subject: [PATCH] Documentation updates. --- docs/api_FDEM.rst | 1 - docs/api_TDEM.rst | 5 +- simpegEM/TDEM/SurveyTDEM.py | 126 ------------------ simpegEM/TDEM/TDEM_b.py | 110 +++++++++++++-- simpegEM/Tests/test_TDEM_b_DerivAdjoint.py | 14 +- .../Tests/test_TDEM_b_MultiTx_DerivAdjoint.py | 2 +- 6 files changed, 109 insertions(+), 149 deletions(-) diff --git a/docs/api_FDEM.rst b/docs/api_FDEM.rst index faa93cd0..2c063216 100644 --- a/docs/api_FDEM.rst +++ b/docs/api_FDEM.rst @@ -107,4 +107,3 @@ FDEM Survey :show-inheritance: :members: :undoc-members: - diff --git a/docs/api_TDEM.rst b/docs/api_TDEM.rst index 6dbd5446..a77ed227 100644 --- a/docs/api_TDEM.rst +++ b/docs/api_TDEM.rst @@ -55,13 +55,12 @@ TDEM - B formulation :show-inheritance: :members: :undoc-members: - :inherited-members: Field Storage ============= -.. automodule:: simpegEM.TDEM.FieldsTDEM +.. autoclass:: simpegEM.TDEM.SurveyTDEM.FieldsTDEM :show-inheritance: :members: :undoc-members: @@ -71,7 +70,7 @@ Field Storage TDEM Survey Classes =================== -.. automodule:: simpegEM.TDEM.SurveyTDEM +.. autoclass:: simpegEM.TDEM.SurveyTDEM.SurveyTDEM :show-inheritance: :members: :undoc-members: diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 3df191cd..34bbfe42 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -135,129 +135,3 @@ class SurveyTDEM(Survey.BaseSurvey): return f - -# class SurveyTDEM1D(BaseSurvey): -# """ -# docstring for SurveyTDEM1D -# """ - -# txLoc = None #: txLoc -# txType = None #: txType -# rxLoc = None #: rxLoc -# rxType = None #: rxType -# timeCh = None #: timeCh -# nTx = 1 #: Number of transmitters - -# @property -# def nTimeCh(self): -# """Number of time channels""" -# return self.timeCh.size - -# def __init__(self, **kwargs): -# BaseSurvey.__init__(self, **kwargs) -# Utils.setKwargs(self, **kwargs) - -# def projectFields(self, u): -# #TODO: this is hardcoded to 1Tx -# return self.Qrx.dot(u.b[:,:,0].T).T - -# def projectFieldsAdjoint(self, d): -# # TODO: make the following self.nTimeCh -# d = d.reshape((self.prob.nT, self.nTx), order='F') -# #TODO: *Qtime.T need to multiply by a time projection. (outside for loop??) -# ii = 0 -# F = FieldsTDEM(self.prob.mesh, self.nTx, self.prob.nT, 'b') -# for ii in range(self.prob.nT): -# b = self.Qrx.T*d[ii,:] -# F.set_b(b, ii) -# F.set_e(np.zeros((self.prob.mesh.nE,self.nTx)), ii) -# return F - -# #################################################### -# # Interpolation Matrices -# #################################################### - -# @property -# def Qrx(self): -# if self._Qrx is None: -# if self.rxType == 'bz': -# locType = 'Fz' -# self._Qrx = self.prob.mesh.getInterpolationMat(self.rxLoc, locType=locType) -# return self._Qrx -# _Qrx = None - - -# class FieldsTDEM_OLD(object): -# """docstring for FieldsTDEM""" - -# phi0 = None #: Initial electric potential -# A0 = None #: Initial magnetic vector potential -# e0 = None #: Initial electric field -# b0 = None #: Initial magnetic flux density -# j0 = None #: Initial current density -# h0 = None #: Initial magnetic field - -# phi = None #: Electric potential -# A = None #: Magnetic vector potential -# e = None #: Electric field -# b = None #: Magnetic flux density -# j = None #: Current density -# h = None #: Magnetic field - -# def __init__(self, mesh, nTx, nT, store='b'): - -# self.nT = nT #: Number of times -# self.nTx = nTx #: Number of transmitters -# self.mesh = mesh - -# def update(self, newFields, tInd): -# self.set_b(newFields['b'], tInd) -# self.set_e(newFields['e'], tInd) - -# def fieldVec(self): -# u = np.ndarray((0, self.nTx)) -# for i in range(self.nT): -# u = np.r_[u, self.get_b(i), self.get_e(i)] -# if self.nTx == 1: -# u = u.flatten() -# return u - -# #################################################### -# # Get Methods -# #################################################### - -# def get_b(self, ind): -# if ind == -1: -# return self.b0 -# else: -# return self.b[ind,:,:] - -# def get_e(self, ind): -# if ind == -1: -# return self.e0 -# else: -# return self.e[ind,:,:] - -# #################################################### -# # Set Methods -# #################################################### - -# def set_b(self, b, ind): -# if self.b is None: -# self.b = np.zeros((self.nT, np.sum(self.mesh.nF), self.nTx)) -# self.b[:] = np.nan -# if len(b.shape) == 1: -# b = b[:, np.newaxis] -# self.b[ind,:,:] = b - -# def set_e(self, e, ind): -# if self.e is None: -# self.e = np.zeros((self.nT, np.sum(self.mesh.nE), self.nTx)) -# self.e[:] = np.nan -# if len(e.shape) == 1: -# e = e[:, np.newaxis] -# self.e[ind,:,:] = e - - -# def __contains__(self, key): -# return key in self.children diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 207baf8a..65695d00 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -18,7 +18,7 @@ class ProblemTDEM_b(BaseTDEMProblem): def __init__(self, mesh, mapping=None, **kwargs): BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs) - solType = 'b' + solType = 'b' #: Type of the solution, in this case the 'b' field surveyPair = SurveyTDEM @@ -47,16 +47,42 @@ class ProblemTDEM_b(BaseTDEMProblem): #################################################### def Jvec(self, m, v, u=None): - if u is None: - u = self.fields(m) + """ + :param numpy.array m: Conductivity model + :param numpy.ndarray v: vector (model object) + :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m + :rtype: numpy.ndarray + :return: w (data object) + + Multiplying \\\(\\\mathbf{J}\\\) onto a vector can be broken into three steps + + * Compute \\\(\\\\vec{p} = \\\mathbf{G}v\\\) + * Solve \\\(\\\hat{\\\mathbf{A}} \\\\vec{y} = \\\\vec{p}\\\) + * Compute \\\(\\\\vec{w} = -\\\mathbf{Q} \\\\vec{y}\\\) + + """ + u = u or self.fields(m) p = self.Gvec(m, v, u) y = self.solveAh(m, p) Jv = self.survey.projectFieldsDeriv(u, v=y) return - mkvc(Jv) def Jtvec(self, m, v, u=None): - if u is None: - u = self.fields(m) + """ + :param numpy.array m: Conductivity model + :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object) + :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m + :rtype: numpy.ndarray + :return: w (model object) + + Multiplying \\\(\\\mathbf{J}^\\\\top\\\) onto a vector can be broken into three steps + + * Compute \\\(\\\\vec{p} = \\\mathbf{Q}^\\\\top \\\\vec{v}\\\) + * Solve \\\(\\\hat{\\\mathbf{A}}^\\\\top \\\\vec{y} = \\\\vec{p}\\\) + * Compute \\\(\\\\vec{w} = -\\\mathbf{G}^\\\\top y\\\) + + """ + u = u or self.fields(m) if not isinstance(v, self.dataPair): v = self.dataPair(self.survey, v) @@ -76,8 +102,7 @@ class ProblemTDEM_b(BaseTDEMProblem): Multiply G by a vector """ - if u is None: - u = self.fields(m) + u = u or self.fields(m) # Note: Fields has shape (nF/E, nTx, nT+1) # However, p will only really fill (:,:,1:nT+1) @@ -104,8 +129,7 @@ class ProblemTDEM_b(BaseTDEMProblem): Multiply G.T by a vector """ - if u is None: - u = self.fields(m) + u = u or self.fields(m) nTx, nE = self.survey.nTx, self.mesh.nE tmp = np.zeros(nE) # Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go. @@ -120,6 +144,38 @@ class ProblemTDEM_b(BaseTDEMProblem): return p def solveAh(self, m, p): + """ + :param numpy.array m: Conductivity model + :param simpegEM.TDEM.FieldsTDEM p: Fields object + :rtype: simpegEM.TDEM.FieldsTDEM + :return: y + + Solve the block-matrix system \\\(\\\hat{A} \\\hat{y} = \\\hat{p}\\\): + + .. math:: + \mathbf{\hat{A}} = \left[ + \\begin{array}{cccc} + A & 0 & & \\\\ + B & A & & \\\\ + & \ddots & \ddots & \\\\ + & & B & A + \end{array} + \\right] \\\\ + \mathbf{A} = + \left[ + \\begin{array}{cc} + \\frac{1}{\delta t} \MfMui & \MfMui\dcurl \\\\ + \dcurl^\\top \MfMui & -\MeSig + \end{array} + \\right] \\\\ + \mathbf{B} = + \left[ + \\begin{array}{cc} + -\\frac{1}{\delta t} \MfMui & 0 \\\\ + 0 & 0 + \end{array} + \\right] \\\\ + """ def AhRHS(tInd, y): rhs = self.MfMui*self.mesh.edgeCurl*self.MeSigmaI*p[:,'e',tInd+1] + p[:,'b',tInd+1] @@ -139,6 +195,38 @@ class ProblemTDEM_b(BaseTDEMProblem): return self.forward(m, AhRHS, AhCalcFields) def solveAht(self, m, p): + """ + :param numpy.array m: Conductivity model + :param simpegEM.TDEM.FieldsTDEM p: Fields object + :rtype: simpegEM.TDEM.FieldsTDEM + :return: y + + Solve the block-matrix system \\\(\\\hat{A}^\\\\top \\\hat{y} = \\\hat{p}\\\): + + .. math:: + \mathbf{\hat{A}}^\\top = \left[ + \\begin{array}{cccc} + A & B & & \\\\ + & \ddots & \ddots & \\\\ + & & A & B \\\\ + & & 0 & A + \end{array} + \\right] \\\\ + \mathbf{A} = + \left[ + \\begin{array}{cc} + \\frac{1}{\delta t} \MfMui & \MfMui\dcurl \\\\ + \dcurl^\\top \MfMui & -\MeSig + \end{array} + \\right] \\\\ + \mathbf{B} = + \left[ + \\begin{array}{cc} + -\\frac{1}{\delta t} \MfMui & 0 \\\\ + 0 & 0 + \end{array} + \\right] \\\\ + """ # Mini Example: # @@ -184,7 +272,7 @@ class ProblemTDEM_b(BaseTDEMProblem): # Functions for tests #################################################### - def AhVec(self, m, vec): + def _AhVec(self, m, vec): """ :param numpy.array m: Conductivity model :param simpegEM.TDEM.FieldsTDEM vec: Fields object @@ -229,7 +317,7 @@ class ProblemTDEM_b(BaseTDEMProblem): f[:,'e',i] = self.mesh.edgeCurl.T*self.MfMui*vec[:,'b',i] - self.MeSigma*vec[:,'e',i] return f - def AhtVec(self, m, vec): + def _AhtVec(self, m, vec): """ :param numpy.array m: Conductivity model :param simpegEM.TDEM.FieldsTDEM vec: Fields object diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 36176944..db1a3dfb 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -48,7 +48,7 @@ class TDEM_bDerivTests(unittest.TestCase): sigma = self.sigma u = prb.fields(sigma) - Ahu = prb.AhVec(sigma, u) + Ahu = prb._AhVec(sigma, u) V1 = Ahu[:,'b',1] V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0] @@ -87,7 +87,7 @@ class TDEM_bDerivTests(unittest.TestCase): f = prb.fields(sigma) u1 = A*f.tovec() - u2 = prb.AhVec(sigma,f).tovec() + u2 = prb._AhVec(sigma,f).tovec() self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12) @@ -130,7 +130,7 @@ class TDEM_bDerivTests(unittest.TestCase): for i in range(prb.nT): f[:,'e', i] = np.random.rand(mesh.nE, 1) - Ahf = prb.AhVec(sigma, f) + Ahf = prb._AhVec(sigma, f) f_test = prb.solveAh(sigma, Ahf) u1 = f.tovec() @@ -149,7 +149,7 @@ class TDEM_bDerivTests(unittest.TestCase): dm = 1000*np.random.rand(self.prb.mapping.nP) h = 0.01 - derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] + derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] print '\ntest_DerivG' passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) self.assertTrue(passed) @@ -221,8 +221,8 @@ class TDEM_bDerivTests(unittest.TestCase): f2[:,'b',i] = np.random.rand(mesh.nF, 1) f2[:,'e',i] = np.random.rand(mesh.nE, 1) - V1 = f2.tovec().dot(prb.AhVec(sigma, f1).tovec()) - V2 = f1.tovec().dot(prb.AhtVec(sigma, f2).tovec()) + V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec()) + V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec()) self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) # def test_solveAhtVsAhtVec(self): @@ -236,7 +236,7 @@ class TDEM_bDerivTests(unittest.TestCase): # f1[:,'e',i] = np.random.rand(mesh.nE, 1) # f2 = prb.solveAht(sigma, f1) - # f3 = prb.AhtVec(sigma, f2) + # f3 = prb._AhtVec(sigma, f2) # if True: # import matplotlib.pyplot as plt diff --git a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py index 477952e1..1d7ea8f2 100644 --- a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py @@ -53,7 +53,7 @@ class TDEM_bDerivTests(unittest.TestCase): dm = 1000*np.random.rand(self.prb.mapping.nP) h = 0.01 - derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] + derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] print '\ntest_DerivG' passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) self.assertTrue(passed)