try using list for Jv instead of datapair - has been seen to cause memory leaks

- allow vector containing values in the inactive set to be passed (not just nC in length)

SimPEG/Directives.py

@@ -144,7 +144,6 @@ class BetaSchedule(InversionDirective):
             if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
             self.invProb.beta /= self.coolingFactor
 
-
 class TargetMisfit(InversionDirective):
 
     chifact = 1.
@@ -243,6 +242,12 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         # Save the file as a npz
         np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
 
+
+# class UpdateReferenceModel(Parameter):
+
+#     mref0 = None
+
+#     def nextIter(self):
 #         mref = getattr(self, 'm_prev', None)
 #         if mref is None:
 #             if self.debug: print 'UpdateReferenceModel is using mref0'
@@ -253,138 +258,56 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 class Update_IRLS(InversionDirective):
 
     eps_min = None
-    eps_p = None
-    eps_q = None
-    norms = [2.,2.,2.,2.]
     factor = None
     gamma = None
     phi_m_last = None
     phi_d_last = None
-    f_old = None
-    f_min_change = 1e-2
-    beta_tol = 5e-2
 
-    # Solving parameter for IRLS (mode:2)
-    IRLSiter   = 0
-    minGNiter = 5
-    maxIRLSiter = 10
-    iterStart = 0
-
-    # Beta schedule
-    coolingFactor = 2.
-    coolingRate = 1
-
-    mode = 1
-
-    @property
-    def target(self):
-        if getattr(self, '_target', None) is None:
-            self._target = self.survey.nD*0.5
-        return self._target
-    @target.setter
-    def target(self, val):
-        self._target = val
 
     def initialize(self):
 
-        if self.mode == 1:
-            self.reg.norms = [2., 2., 2., 2.]
+        # Scale the regularization for changes in norm
+        if getattr(self, 'phi_m_last', None) is not None:
+
+            self.reg.curModel = self.invProb.curModel
+            self.reg.gamma = 1.
+            phim_new = self.reg.eval(self.invProb.curModel)
+            self.gamma = self.phi_m_last / phim_new
+
+        self.reg.curModel = self.invProb.curModel
+        self.reg.gamma = self.gamma
+
+        if getattr(self, 'phi_d_last', None) is None:
+            self.phi_d_last = self.invProb.phi_d
 
     def endIter(self):
+        # Cool the threshold parameter if required
+        if getattr(self, 'factor', None) is not None:
+            eps = self.reg.eps / self.factor
 
-        # After reaching target misfit with l2-norm, switch to IRLS (mode:2)
-        if self.invProb.phi_d < self.target and self.mode == 1:
-            print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
-
-            self.mode = 2
-            print self.eps_p, self.eps_q, self.norms
-            self.reg.eps_p = self.eps_p
-            self.reg.eps_q = self.eps_q
-            self.reg.norms = self.norms
-            self.coolingFactor = 1.
-            self.coolingRate = 1
-            self.iterStart = self.opt.iter
-            self.phi_d_last = self.invProb.phi_d
-            self.phi_m_last = self.invProb.phi_m_last
-
-            self.reg.l2model = self.invProb.curModel
-            self.reg.curModel = self.invProb.curModel
-
-            if getattr(self, 'f_old', None) is None:
-                self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False)
-
-        # Beta Schedule
-        if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
-            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
-            self.invProb.beta /= self.coolingFactor
-
-
-        # Only update after GN iterations
-        if (self.opt.iter-self.iterStart) % self.minGNiter == 0 and self.mode==2:
-
-            self.IRLSiter += 1
-
-            phim_new = self.reg.eval(self.invProb.curModel)
-            self.f_change = np.abs(self.f_old - phim_new) / self.f_old
-
-            print "Regularization decrease: %6.3e" % (self.f_change)
-
-            # Check for maximum number of IRLS cycles
-            if self.IRLSiter == self.maxIRLSiter:
-                print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
-                self.opt.stopNextIteration = True
-                return
-
-            # Check if the function has changed enough
-            if self.f_change < self.f_min_change and self.IRLSiter > 1:
-                print "Minimum decrease in regularization. End of IRLS"
-                self.opt.stopNextIteration = True
-                return
+            if getattr(self, 'eps_min', None) is not None:
+                self.reg.eps = np.max([self.eps_min,eps])
             else:
-                self.f_old = phim_new
+                self.reg.eps = eps
 
-            # Cool the threshold parameter if required
-            if getattr(self, 'factor', None) is not None:
-                eps = self.reg.eps / self.factor
+        # Get phi_m at the end of current iteration
+        self.phi_m_last = self.invProb.phi_m_last
 
-                if getattr(self, 'eps_min', None) is not None:
-                    self.reg.eps = np.max([self.eps_min,eps])
-                else:
-                    self.reg.eps = eps
+         # Update the model used for the IRLS weights
+        self.reg.curModel = self.invProb.curModel
 
-            # Get phi_m at the end of current iteration
-            self.phi_m_last = self.invProb.phi_m_last
+        # Temporarely set gamma to 1. to get raw phi_m
+        self.reg.gamma = 1.
 
-            # Reset the regularization matrices so that it is
-            # recalculated for current model
-            self.reg._Wsmall = None
-            self.reg._Wx = None
-            self.reg._Wy = None
-            self.reg._Wz = None
+        # Compute new model objective function value
+        phim_new = self.reg.eval(self.invProb.curModel)
 
-             # Update the model used for the IRLS weights
-            self.reg.curModel = self.invProb.curModel
+        # Update gamma to scale the regularization between IRLS iterations
+        self.reg.gamma = self.phi_m_last / phim_new
 
-            # Temporarely set gamma to 1. to get raw phi_m
-            self.reg.gamma = 1.
-
-            # Compute new model objective function value
-            phim_new = self.reg.eval(self.invProb.curModel)
-
-            # Update gamma to scale the regularization between IRLS iterations
-            self.reg.gamma = self.phi_m_last / phim_new
-
-            # Reset the regularization matrices again for new gamma
-            self.reg._Wsmall = None
-            self.reg._Wx = None
-            self.reg._Wy = None
-            self.reg._Wz = None
-
-            # Check if misfit is within the tolerance, otherwise scale beta
-            val = self.invProb.phi_d / (self.survey.nD*0.5)
-
-            if np.abs(1.-val) > self.beta_tol:
-                self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
+        # Set the weighting matrix to None so that it is recomputed next time
+        # it is called in the inversion
+        self.reg._W = None
 
 class Update_lin_PreCond(InversionDirective):
     """
@@ -437,3 +360,19 @@ class Update_Wj(InversionDirective):
             JtJdiag = JtJdiag / max(JtJdiag)
 
             self.reg.wght = JtJdiag
+
+class Scale_Beta(InversionDirective):
+    """
+        Instead of a linear cooling schedule, beta is allowed to change based
+        on the ratio between the target misfit and the current data misfit. The
+        update is done only if the misfit is outside some threshold bounds.
+    """
+    tol = 0.05
+
+    def endIter(self):
+
+        # Check if misfit is within the tolerance, otherwise adjust beta
+        val = self.invProb.phi_d / (self.survey.nD*0.5)
+
+        if np.abs(1.-val) > self.tol:
+            self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -60,20 +60,6 @@ class Fields(SimPEG.Problem.Fields):
 
         return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
 
-    def _bSecondary(self, solution, srcList):
-        """
-        Total magnetic flux density is sum of primary and secondary
-
-        :param numpy.ndarray solution: field we solved for
-        :param list srcList: list of sources
-        :rtype: numpy.ndarray
-        :return: total magnetic flux density
-        """
-        if getattr(self, '_bSecondary', None) is None:
-            raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
-
-        return  self._bSecondary(solution, srcList)
-
     def _h(self, solution, srcList):
         """
         Total magnetic field is sum of primary and secondary
@@ -138,21 +124,6 @@ class Fields(SimPEG.Problem.Fields):
             return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
         return np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex)
 
-    def _bSecondaryDeriv(self, src, du_dm_v, v, adjoint = False):
-        """
-        Total derivative of b with respect to the inversion model. Returns :math:`d\mathbf{b}/d\mathbf{m}` for forward and (:math:`d\mathbf{b}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
-
-        :param Src src: sorce
-        :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
-        :param numpy.ndarray v: vector to take sensitivity product with
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: derivative times a vector (or tuple for adjoint)
-        """
-        # TODO: modify when primary field is dependent on m
-
-        return self._bDeriv(src, du_dm_v, v, adjoint = adjoint)
-
     def _hDeriv(self, src, du_dm_v, v, adjoint = False):
         """
         Total derivative of h with respect to the inversion model. Returns :math:`d\mathbf{h}/d\mathbf{m}` for forward and (:math:`d\mathbf{h}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
@@ -286,7 +257,7 @@ class Fields3D_e(Fields):
         """
 
         # assuming primary does not depend on the model
-        return Zero()
+        return src.ePrimaryDeriv(self.prob, v, adjoint) #Zero()
 
     def _bPrimary(self, eSolution, srcList):
         """
@@ -500,8 +471,6 @@ class Fields3D_b(Fields):
             return 'E'
         elif fieldType == 'b':
             return 'F'
-        elif fieldType == 'bSecondary':
-            return 'F'
         elif (fieldType == 'h') or (fieldType == 'j'):
             return'CCV'
         else:
@@ -631,8 +600,8 @@ class Fields3D_b(Fields):
 
 
         if adjoint:
-            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv
-        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv
+            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv + src.ePrimaryDeriv(self.prob, v, adjoint)
+        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv + src.ePrimaryDeriv(self.prob, v, adjoint)
 
     def _j(self, bSolution, srcList):
         """

SimPEG/EM/FDEM/ProblemFDEM.py

@@ -74,7 +74,8 @@ class BaseFDEMProblem(BaseEMProblem):
 
         self.curModel = m
 
-        Jv = self.dataPair(self.survey)
+        # Jv = self.dataPair(self.survey)
+        Jv = []
 
         for freq in self.survey.freqs:
             A = self.getA(freq)
@@ -89,9 +90,9 @@ class BaseFDEMProblem(BaseEMProblem):
                 for rx in src.rxList:
                     df_dmFun = getattr(f, '_{0}Deriv'.format(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)
+                    Jv.append(rx.evalDeriv(src, self.mesh, f, df_dm_v))
             Ainv.clean()
-        return Utils.mkvc(Jv)
+        return np.hstack(Jv)
 
     def Jtvec(self, m, v, f=None):
         """
@@ -166,7 +167,6 @@ class BaseFDEMProblem(BaseEMProblem):
 
         for i, src in enumerate(Srcs):
             smi, sei = src.eval(self)
-            #Why are you adding?
             s_m[:,i] = s_m[:,i] + smi
             s_e[:,i] = s_e[:,i] + sei
 

SimPEG/EM/FDEM/RxFDEM.py

@@ -97,19 +97,6 @@ class Point_b(BaseRx):
         self.projField = 'b'
         super(Point_b, self).__init__(locs, orientation, component)
 
-class Point_bSecondary(BaseRx):
-    """
-    Magnetic flux FDEM receiver
-
-    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
-    :param string orientation: receiver orientation 'x', 'y' or 'z'
-    :param string component: real or imaginary component 'real' or 'imag'
-    """
-
-    def __init__(self, locs, orientation=None, component=None):
-        self.projField = 'bSecondary'
-        super(Point_bSecondary, self).__init__(locs, orientation, component)
-
 
 class Point_h(BaseRx):
     """

SimPEG/EM/FDEM/SrcFDEM.py

@@ -60,6 +60,18 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._bPrimary
 
+    def bPrimaryDeriv(self, prob, v, adjoint=False):
+        """
+        Derivative of the primary magnetic flux density
+
+        :param Problem prob: FDEM Problem
+        :param numpy.ndarray v: vector
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: primary magnetic flux density
+        """
+        return Zero()
+
     def hPrimary(self, prob):
         """
         Primary magnetic field
@@ -72,6 +84,18 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._hPrimary
 
+    def hPrimaryDeriv(self, prob, v, adjoint=False):
+        """
+        Derivative of the primary magnetic field
+
+        :param Problem prob: FDEM Problem
+        :param numpy.ndarray v: vector
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: primary magnetic flux density
+        """
+        return Zero()
+
     def ePrimary(self, prob):
         """
         Primary electric field
@@ -84,6 +108,18 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._ePrimary
 
+    def ePrimaryDeriv(self, prob, v, adjoint=False):
+        """
+        Derivative of the primary electric field
+
+        :param Problem prob: FDEM Problem
+        :param numpy.ndarray v: vector
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: primary magnetic flux density
+        """
+        return Zero()
+
     def jPrimary(self, prob):
         """
         Primary current density
@@ -96,6 +132,18 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._jPrimary
 
+    def jPrimaryDeriv(self, prob, v, adjoint=False):
+        """
+        Derivative of the primary current density
+
+        :param Problem prob: FDEM Problem
+        :param numpy.ndarray v: vector
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: primary magnetic flux density
+        """
+        return Zero()
+
     def s_m(self, prob):
         """
         Magnetic source term
@@ -614,5 +662,155 @@ class CircularLoop(BaseSrc):
             return -C.T * (MMui_s * self.bPrimary(prob))
 
 
+class PrimSecSigma(BaseSrc):
+
+    def __init__(self, rxList, freq, sigBack, ePrimary, **kwargs):
+        self.sigBack = sigBack
+
+        BaseSrc.__init__(self, rxList, freq=freq, _ePrimary=ePrimary, **kwargs)
+
+    def s_e(self, prob):
+        return (prob.MeSigma -  prob.mesh.getEdgeInnerProduct(self.sigBack)) * self.ePrimary(prob)
+
+    def s_eDeriv(self, prob, v, adjoint=False):
+        if adjoint:
+            return prob.MeSigmaDeriv(self.ePrimary(prob)).T * v
+        return prob.MeSigmaDeriv(self.ePrimary(prob)) * v
+
+
+class PrimSecMappedSigma(BaseSrc):
+
+    """
+    Primary-Secondary Source in which a mapping is provided to put the current model
+    onto the primary mesh. This is solved on every model update.
+
+    There are a lot of layers to the derivatives here!
+
+    **Required**
+    :param list rxList: Receiver List
+    :param float freq: frequency
+    :param ProblemFDEM primaryProblem: FDEM primary problem
+    :param SurveyFDEM primarySurvey: FDEM primary survey
+
+    **Optional**
+    :param Mapping map2meshSecondary: mapping current model to act as primary model on the secondary mesh
+
+    """
+
+    def __init__(self, rxList, freq, primaryProblem, primarySurvey, map2meshSecondary = None ,**kwargs):
+
+        self.primaryProblem = primaryProblem
+        self.primarySurvey = primarySurvey
+
+        if self.primaryProblem.ispaired is False:
+            self.primaryProblem.pair(self.primarySurvey)
+
+        self.map2meshSecondary = map2meshSecondary
+
+        BaseSrc.__init__(self, rxList, freq=freq, **kwargs)
+
+    def _ProjPrimary(self, prob):
+        # if getattr(self, '__ProjPrimary', None) is None:
+        return self.primaryProblem.mesh.getInterpolationMatCartMesh(prob.mesh, locType='F', locTypeTo='E')
+        # return self.__ProjPrimary
+
+
+    def _primaryFields(self, prob, fieldType=None):
+
+        # TODO: cache and check if prob.curModel has changed
+        fields = self.primaryProblem.fields(prob.curModel.sigmaModel)
+
+        if fieldType is not None:
+            return fields[:,fieldType]
+        return fields
+
+    def _primaryFieldsDeriv(self, prob, v, adjoint=False, f=None):
+        if adjoint:
+            raise NotImplementedError
+
+        # TODO: this should not be hard-coded for j
+        # jp = self._primaryFields(prob)[:,'j']
+
+        # TODO: pull apart Jvec so that don't have to copy paste this code in
+        # A = self.primaryProblem.getA(self.freq)
+        # Ainv = self.primaryProblem.Solver(A, **self.primaryProblem.solverOpts) # create the concept of Ainv (actually a solve)
+
+        if f is None:
+           f = self._primaryFields(prob.curModel.sigmaModel)
+
+        freq = self.freq
+
+        A = self.primaryProblem.getA(freq)
+        Ainv = self.primaryProblem.Solver(A, **self.primaryProblem.solverOpts) # create the concept of Ainv (actually a solve)
+
+        src = self.primarySurvey.srcList[0]
+        # for src in self.survey.getSrcByFreq(freq):
+        u_src = Utils.mkvc(f[src, self.primaryProblem._solutionType])
+        dA_dm_v = self.primaryProblem.getADeriv(freq, u_src, v)
+        dRHS_dm_v = self.primaryProblem.getRHSDeriv(freq, src, v)
+        du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
+
+        df_dmFun = getattr(f, '_{0}Deriv'.format('j'), 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)
+        Ainv.clean()
+
+        return df_dm_v
+
+        # return self.primaryProblem.Jvec(prob.curModel, v, f=f)
+
+    def ePrimary(self, prob, f=None):
+        if f is None:
+            f = self._primaryFields(prob)
+
+        ep = self._ProjPrimary(prob) * (
+                self.primaryProblem.MfI * (
+                self.primaryProblem.MfRho * f[:,'j'])
+                )
+
+        return Utils.mkvc(ep)
+
+    def ePrimaryDeriv(self, prob, v, adjoint=False, f=None):
+
+        if adjoint is True:
+            raise NotImplementedError
+
+        if f is None:
+            f = self._primaryFields(prob)
+
+        epDeriv = self._ProjPrimary(prob) * (
+                    self.primaryProblem.MfI * (
+                    (self.primaryProblem.MfRhoDeriv(f[:,'j']) * v)
+                    +
+                    (self.primaryProblem.MfRho * self._primaryFieldsDeriv(prob, v, f=f))
+                    )
+                    )
+
+        return Utils.mkvc(epDeriv)
+
+
+    def s_e(self, prob):
+        sigmaPrimary = self.map2meshSecondary * prob.curModel.sigmaModel
+
+        return Utils.mkvc((prob.MeSigma -  prob.mesh.getEdgeInnerProduct(sigmaPrimary)) * self.ePrimary(prob))
+
+
+    def s_eDeriv(self, prob, v, adjoint=False):
+        if adjoint:
+            raise NotImplementedError
+            return prob.MeSigmaDeriv(self.ePrimary(prob)).T * v
+
+        sigmaPrimary = self.map2meshSecondary * prob.curModel.sigmaModel
+        sigmaPrimaryDeriv = self.map2meshSecondary.deriv(prob.curModel.sigmaModel)
+
+        f = self._primaryFields(prob)
+        ePrimary = self.ePrimary(prob,f=f)
+
+        return (prob.MeSigmaDeriv(ePrimary) * v
+                - prob.mesh.getEdgeInnerProductDeriv(sigmaPrimary)(ePrimary) * sigmaPrimaryDeriv * v
+                + (prob.MeSigma -  prob.mesh.getEdgeInnerProduct(sigmaPrimary)) * self.ePrimaryDeriv(prob, v, None, f=f)
+                )
+
+
 
 

SimPEG/EM/Static/DC/RxDC.py

@@ -43,7 +43,14 @@ class BaseRx(SimPEG.Survey.BaseRx):
         elif adjoint:
             return P.T*v
 
-# DC.Rx.Dipole(locs)
+# DC.Rx.Pole(locs)
+class Pole(BaseRx):
+
+    def __init__(self, locs, rxType = 'phi', **kwargs):
+        BaseRx.__init__(self, locs, rxType)
+
+
+# DC.Rx.Dipole(locsM, locsN)
 class Dipole(BaseRx):
 
     def __init__(self, locsM, locsN, rxType = 'phi', **kwargs):

SimPEG/EM/TDEM/FieldsTDEM.py

@@ -1,142 +0,0 @@
-import numpy as np
-import scipy.sparse as sp
-import SimPEG
-from SimPEG import Utils
-from SimPEG.EM.Utils import omega
-from SimPEG.Utils import Zero, Identity
-
-class Fields(SimPEG.Problem.TimeFields):
-    """
-
-    Fancy Field Storage for a TDEM survey. Only one field type is stored for
-    each problem, the rest are computed. The fields obejct acts like an array and is indexed by
-
-    .. code-block:: python
-
-        f = problem.fields(m)
-        e = f[srcList,'e']
-        b = f[srcList,'b']
-
-    If accessing all sources for a given field, use the :code:`:`
-
-    .. code-block:: python
-
-        f = problem.fields(m)
-        e = f[:,'e']
-        b = f[:,'b']
-
-    The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies)
-    """
-
-    knownFields = {}
-    dtype = float
-
-    def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
-        if adjoint is True:
-            return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint)
-        return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v)
-
-    def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
-        if adjoint is True:
-            return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
-        return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
-
-class Fields_Derivs(Fields):
-    knownFields = {
-                    'bDeriv': 'F',
-                    'eDeriv': 'E',
-                    'hDeriv': 'E',
-                    'jDeriv': 'F'
-                  }
-
-
-class Fields_b(Fields):
-    """Fancy Field Storage for a TDEM survey."""
-    knownFields = {'bSolution': 'F'}
-    aliasFields = {
-                    'b': ['bSolution', 'F', '_b'],
-                    'e': ['bSolution', 'E', '_e'],
-                  }
-
-    def startup(self):
-        self.MeSigmaI      = self.survey.prob.MeSigmaI
-        self.MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
-        self.edgeCurl      = self.survey.prob.mesh.edgeCurl
-        self.MfMui         = self.survey.prob.MfMui
-
-    def _b(self, bSolution, srcList, tInd):
-        return bSolution
-
-    def _bDeriv_u(self, tInd, src, dun_dm_v, adjoint=False):
-        return Identity()*dun_dm_v
-
-    def _bDeriv_m(self, tInd, src, v, adjoint=False):
-        return Zero()
-
-    # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
-    #     if adjoint is True:
-    #         return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
-    #     return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
-
-    def _e(self, bSolution, srcList, tInd):
-        e = self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) )
-        for i, src in enumerate(srcList):
-            _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
-            e[:,i] = e[:,i] - self.MeSigmaI * S_e
-        return e
-
-    def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False):
-        if adjoint is True:
-            return self.MfMui.T * ( self.edgeCurl * ( self.MeSigmaI.T * dun_dm_v ) )
-        return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * dun_dm_v ) )
-
-    def _eDeriv_m(self, tInd, src, v, adjoint = False):
-        _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
-        bSolution = self[[src],'bSolution',tInd]
-
-        _, S_eDeriv = src.evalDeriv(self.survey.prob.times[tInd], self, adjoint=adjoint)
-
-        if adjoint is True:
-            return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution ) ).T * v - S_eDeriv(self.MeSigmaI.T * v)
-
-        return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution)) * v - self.MeSigmaI * S_eDeriv(v)
-
-
-
-class Fields_e(Fields):
-    """Fancy Field Storage for a TDEM survey."""
-    knownFields = {'eSolution': 'E'}
-    aliasFields = {
-                    'e': ['eSolution', 'E', '_e'],
-                    'b': ['eSolution', 'F', '_b'],
-                  }
-
-    def startup(self):
-        self.MeSigmaI      = self.survey.prob.MeSigmaI
-        self.MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
-        self.edgeCurl      = self.survey.prob.mesh.edgeCurl
-        self.MfMui         = self.survey.prob.MfMui
-
-
-    def _e(self, eSolution, srcList, tInd):
-        return eSolution
-
-    def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False):
-        return dun_dm_v
-
-    def _eDeriv_m(self, tInd, src, v, adjoint = False):
-        return Zero()
-
-    def _b(self, eSolution, srcList, tInd):
-        raise NotImplementedError
-
-    def _bDeriv_u(self, tInd, src, dun_dm_v, adjoint=False):
-        raise NotImplementedError
-
-    def _bDeriv_m(self, tInd, src, v, adjoint=False):
-        raise NotImplementedError
-
-    # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
-    #     if adjoint is True:
-    #         return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
-    #     return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)

SimPEG/EM/TDEM/SrcTDEM.py

@@ -1,246 +0,0 @@
-import SimPEG
-from SimPEG import np, Utils
-from SimPEG.Utils import Zero, Identity
-from scipy.constants import mu_0
-from SimPEG.EM.Utils import *
-
-####################################################
-# Sources
-####################################################
-
-class BaseWaveform(object):
-
-    def __init__(self, offTime=0., hasInitialFields=False):
-        self.offTime = offTime
-        self.hasInitialFields = hasInitialFields
-
-    def _assertMatchesPair(self, pair):
-        assert (isinstance(self, pair)
-            ), "Waveform object must be an instance of a %s BaseWaveform class."%(pair.__name__)
-
-    def eval(self, time):
-        raise NotImplementedError
-
-    def evalDeriv(self, time):
-        raise NotImplementedError # needed for E-formulation
-
-
-class StepOffWaveform(BaseWaveform):
-
-    def __init__(self, offTime=0.):
-        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
-
-    def eval(self, time):
-        return 0.
-
-
-class RawWaveform(BaseWaveform):
-
-    def __init__(self, offTime=0.):
-        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
-
-    def eval(self, time):
-        raise NotImplementedError('RawWaveform has not been implemented, you should write it!')
-
-
-class TriangularWaveform(BaseWaveform):
-
-    def __init__(self, offTime=0.):
-        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
-
-    def eval(self, time):
-        raise NotImplementedError('TriangularWaveform has not been implemented, you should write it!')
-
-
-
-class BaseSrc(SimPEG.Survey.BaseSrc):
-
-    # rxPair = Rx
-    integrate = True
-    waveformPair = BaseWaveform
-
-    @property
-    def waveform(self):
-        "A waveform instance is not None"
-        return getattr(self, '_waveform', None)
-    @waveform.setter
-    def waveform(self, val):
-        if self.waveform is None:
-            val._assertMatchesPair(self.waveformPair)
-            self._mapping = val
-        else:
-            self._mapping = self.PropMap(val)
-
-
-    def __init__(self, rxList, waveform = StepOffWaveform(), **kwargs):
-        self.waveform = waveform
-        SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs)
-
-
-    def bInitial(self, prob):
-        return Zero()
-
-    def bInitialDeriv(self, prob, v=None, adjoint=False):
-        return Zero()
-
-    def eInitial(self, prob):
-        return Zero()
-
-    def eInitialDeriv(self, prob, v=None, adjoint=False):
-        return Zero()
-
-    def eval(self, prob, time):
-        S_m = self.S_m(prob, time)
-        S_e = self.S_e(prob, time)
-        return S_m, S_e
-
-    def evalDeriv(self, prob, time, v=None, adjoint=False):
-        if v is not None:
-            return self.S_mDeriv(prob, time, v, adjoint), self.S_eDeriv(prob, time, v, adjoint)
-        else:
-            return lambda v: self.S_mDeriv(prob, time, v, adjoint), lambda v: self.S_eDeriv(prob, time, v, adjoint)
-
-    def S_m(self, prob, time):
-        return Zero()
-
-    def S_e(self, prob, time):
-        return Zero()
-
-    def S_mDeriv(self, prob, time, v=None, adjoint=False):
-        return Zero()
-
-    def S_eDeriv(self, prob, time, v=None, adjoint=False):
-        return Zero()
-
-
-class MagDipole(BaseSrc):
-
-    waveform = None
-    loc = None
-    orientation = 'Z'
-    moment = 1.
-    mu = mu_0
-
-    def __init__(self, rxList, **kwargs):
-        assert self.orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
-        self.integrate = False
-        BaseSrc.__init__(self, rxList, **kwargs)
-
-    def _bfromVectorPotential(self, prob):
-        if prob._eqLocs is 'FE':
-            gridX = prob.mesh.gridEx
-            gridY = prob.mesh.gridEy
-            gridZ = prob.mesh.gridEz
-            C = prob.mesh.edgeCurl
-
-        elif prob._eqLocs is 'EF':
-            gridX = prob.mesh.gridFx
-            gridY = prob.mesh.gridFy
-            gridZ = prob.mesh.gridFz
-            C = prob.mesh.edgeCurl.T
-
-
-        if prob.mesh._meshType is 'CYL':
-            if not prob.mesh.isSymmetric:
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
-
-        else:
-            srcfct = MagneticDipoleVectorPotential
-            ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
-            ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
-            az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
-            a = np.concatenate((ax, ay, az))
-
-        return C*a
-
-
-    def bInitial(self, prob):
-
-        if self.waveform.hasInitialFields is False:
-            return Zero()
-
-        return self._bfromVectorPotential(prob)
-
-    def eInitial(self, prob):
-
-        if self.waveform.hasInitialFields is False:
-            return Zero()
-
-        b = self.bInitial(prob)
-        MeSigmaI = prob.MeSigmaI
-        MfMui = prob.MfMui
-        C = prob.mesh.edgeCurl
-
-        return MeSigmaI * (C.T * (MfMui * b))
-
-    def eInitialDeriv(self, prob, v=None, adjoint=False):
-
-        if self.waveform.hasInitialFields is False:
-            return Zero()
-
-        b = self.bInitial(prob)
-        MeSigmaIDeriv = prob.MeSigmaIDeriv
-        MfMui = prob.MfMui
-        C = prob.mesh.edgeCurl
-        S_e = self.S_e(prob, prob.t0)
-
-        # S_e doesn't depend on the model
-
-        if adjoint:
-            return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ).T * v
-
-        return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ) * v
-
-
-    def S_m(self, prob, time):
-        if self.waveform.hasInitialFields is False:
-            raise NotImplementedError
-        return Zero()
-
-    def S_e(self, prob, time):
-        if self.waveform.hasInitialFields is False:
-            raise NotImplementedError
-        return Zero()
-
-class CircularLoop(MagDipole):
-
-    waveform = None
-    loc = None
-    orientation = 'Z'
-    radius = None
-    mu = mu_0
-
-    def __init__(self, rxList, **kwargs):
-        assert self.orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
-        self.integrate = False
-        BaseSrc.__init__(self, rxList, **kwargs)
-
-    def _bfromVectorPotential(self, prob):
-        if prob._eqLocs is 'FE':
-            gridX = prob.mesh.gridEx
-            gridY = prob.mesh.gridEy
-            gridZ = prob.mesh.gridEz
-            C = prob.mesh.edgeCurl
-
-        elif prob._eqLocs is 'EF':
-            gridX = prob.mesh.gridFx
-            gridY = prob.mesh.gridFy
-            gridZ = prob.mesh.gridFz
-            C = prob.mesh.edgeCurl.T
-
-
-        if prob.mesh._meshType is 'CYL':
-            if not prob.mesh.isSymmetric:
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', radius=self.radius, mu=self.mu)
-
-        else:
-            srcfct = MagneticLoopVectorPotential
-            ax = srcfct(self.loc, gridX, 'x', mu=self.mu, radius=self.radius)
-            ay = srcfct(self.loc, gridY, 'y', mu=self.mu, radius=self.radius)
-            az = srcfct(self.loc, gridZ, 'z', mu=self.mu, radius=self.radius)
-            a = np.concatenate((ax, ay, az))
-
-        return C*a
-

SimPEG/EM/TDEM/SurveyTDEM.py

@@ -1,16 +1,10 @@
-import SimPEG
-from SimPEG import np, Utils
-from SimPEG.Utils import Zero, Identity
-from scipy.constants import mu_0
+from SimPEG import Utils, Survey, np
+from SimPEG.Survey import BaseSurvey
 from SimPEG.EM.Utils import *
-import SrcTDEM as Src
+from BaseTDEM import FieldsTDEM
 
 
-####################################################
-# Receivers
-####################################################
-
-class Rx(SimPEG.Survey.BaseTimeRx):
+class RxTDEM(Survey.BaseTimeRx):
 
     knownRxTypes = {
                     'ex':['e', 'Ex', 'N'],
@@ -27,7 +21,7 @@ class Rx(SimPEG.Survey.BaseTimeRx):
                    }
 
     def __init__(self, locs, times, rxType):
-        SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType)
+        Survey.BaseTimeRx.__init__(self, locs, times, rxType)
 
     @property
     def projField(self):
@@ -62,60 +56,144 @@ class Rx(SimPEG.Survey.BaseTimeRx):
         u_part = Utils.mkvc(u[src, self.projField, :])
         return P*u_part
 
-    def evalDeriv(self, src, mesh, timeMesh, v, adjoint=False):
+    def evalDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
         P = self.getP(mesh, timeMesh)
+
         if not adjoint:
-            return P * v #Utils.mkvc(v[src, self.projField+'Deriv', :])
+            return P * Utils.mkvc(v[src, self.projField, :])
         elif adjoint:
-            # dP_dF_T = P.T * v #[src, self]
-            # newshape = (len(dP_dF_T)/timeMesh.nN, timeMesh.nN )
-            return P.T * v #np.reshape(dP_dF_T, newshape, order='F')
+            return P.T * v[src, self]
 
 
-####################################################
-# Survey
-####################################################
+class SrcTDEM(Survey.BaseSrc):
+    rxPair = RxTDEM
+    radius = None
 
-class Survey(SimPEG.Survey.BaseSurvey):
+    def getInitialFields(self, mesh):
+        F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
+        return F0
+
+    def getJs(self, mesh, time):
+        return None
+
+
+class SrcTDEM_VMD_MVP(SrcTDEM):
+
+    def __init__(self,rxList,loc,waveformType="STEPOFF"):
+        self.loc = loc
+        self.waveformType = waveformType
+        SrcTDEM.__init__(self,rxList)
+
+    def getInitialFields(self, mesh):
+        """Vertical magnetic dipole, magnetic vector potential"""
+        if self.waveformType == "STEPOFF":
+            print ">> Step waveform: Non-zero initial condition"
+            if mesh._meshType is 'CYL':
+                if mesh.isSymmetric:
+                    MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
+                else:
+                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            elif mesh._meshType is 'TENSOR':
+                MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
+            else:
+                raise Exception('Unknown mesh for VMD')
+            return {"b": mesh.edgeCurl*MVP}
+        elif self.waveformType == "GENERAL":
+            print ">> General waveform: Zero initial condition"
+            return {"b": np.zeros(mesh.nF)}
+        else:
+            raise NotImplementedError("Only use STEPOFF or GENERAL")
+
+    def getMeS(self, mesh, MfMui):
+        if mesh._meshType is 'CYL':
+            if mesh.isSymmetric:
+                MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
+            else:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+        elif mesh._meshType is 'TENSOR':
+            MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
+        else:
+            raise Exception('Unknown mesh for VMD')
+        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
+
+
+class SrcTDEM_CircularLoop_MVP(SrcTDEM):
+    def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
+        self.loc = loc
+        self.radius = radius
+        self.waveformType = waveformType
+        SrcTDEM.__init__(self,rxList)
+
+    def getInitialFields(self, mesh):
+        """Circular Loop, magnetic vector potential"""
+        if self.waveformType == "STEPOFF":
+            print ">> Step waveform: Non-zero initial condition"
+            if mesh._meshType is 'CYL':
+                if mesh.isSymmetric:
+                    MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
+                else:
+                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            elif mesh._meshType is 'TENSOR':
+                MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
+            else:
+                raise Exception('Unknown mesh for CircularLoop')
+            return {"b": mesh.edgeCurl*MVP}
+        elif self.waveformType == "GENERAL":
+            print ">> General waveform: Zero initial condition"
+            return {"b": np.zeros(mesh.nF)}
+        else:
+            raise NotImplementedError("Only use STEPOFF or GENERAL")
+
+    def getMeS(self, mesh, MfMui):
+        if mesh._meshType is 'CYL':
+            if mesh.isSymmetric:
+                MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
+            else:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+        elif mesh._meshType is 'TENSOR':
+            MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
+        else:
+            raise Exception('Unknown mesh for CircularLoop')
+        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
+
+
+class SurveyTDEM(Survey.BaseSurvey):
     """
-    Time domain electromagnetic survey
+        docstring for SurveyTDEM
     """
-
-    srcPair = Src.BaseSrc
-    rxPair = Rx
+    srcPair = SrcTDEM
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
         self.srcList = srcList
-        SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
+        Survey.BaseSurvey.__init__(self, **kwargs)
 
     def eval(self, u):
-        data = SimPEG.Survey.Data(self)
+        data = Survey.Data(self)
         for src in self.srcList:
             for rx in src.rxList:
                 data[src, rx] = rx.eval(src, self.mesh, self.prob.timeMesh, u)
         return data
 
     def evalDeriv(self, u, v=None, adjoint=False):
-        raise Exception('Use Receivers to project fields deriv.')
-        # assert v is not None, 'v to multiply must be provided.'
+        assert v is not None, 'v to multiply must be provided.'
 
-        # if not adjoint:
-        #     data = SimPEG.Survey.Data(self)
-        #     for src in self.srcList:
-        #         for rx in src.rxList:
-        #             data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
-        #     return data
-        # else:
-        #     f = FieldsTDEM(self.mesh, self)
-        #     for src in self.srcList:
-        #         for rx in src.rxList:
-        #             Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
-        #             Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
-        #             if rx.projField not in f: # first time we are projecting
-        #                 f[src, rx.projField, :] = Ptv
-        #             else: # there are already fields, so let's add to them!
-        #                 f[src, rx.projField, :] += Ptv
-        #     return f
+        if not adjoint:
+            data = Survey.Data(self)
+            for src in self.srcList:
+                for rx in src.rxList:
+                    data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
+            return data
+        else:
+            f = FieldsTDEM(self.mesh, self)
+            for src in self.srcList:
+                for rx in src.rxList:
+                    Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
+                    Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
+                    if rx.projField not in f: # first time we are projecting
+                        f[src, rx.projField, :] = Ptv
+                    else: # there are already fields, so let's add to them!
+                        f[src, rx.projField, :] += Ptv
+            return f
 
 

SimPEG/EM/TDEM/TDEM.py

@@ -1,553 +0,0 @@
-from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
-from SimPEG.EM.Base import BaseEMProblem
-from SimPEG.EM.TDEM.SurveyTDEM import Survey as SurveyTDEM
-from SimPEG.EM.TDEM.FieldsTDEM import *
-from scipy.constants import mu_0
-import time
-
-class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
-    """
-    We start with the first order form of Maxwell's equations
-    """
-    surveyPair = SurveyTDEM
-    fieldsPair = Fields
-
-    def __init__(self, mesh, mapping=None, **kwargs):
-        Problem.BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
-
-    def fields(self, m):
-        """
-        Solve the forward problem for the fields.
-
-        :param numpy.array m: inversion model (nP,)
-        :rtype numpy.array:
-        :return F: fields
-        """
-
-        tic = time.time()
-        self.curModel = m
-
-        F = self.fieldsPair(self.mesh, self.survey)
-
-        # set initial fields
-        F[:,self._fieldType+'Solution',0] = self.getInitialFields()
-
-        # timestep to solve forward
-        if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
-        Ainv = None
-        for tInd, dt in enumerate(self.timeSteps):
-            if Ainv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same
-                Ainv.clean()
-                Ainv = None
-
-            if Ainv is None:
-                A = self.getAdiag(tInd)
-                if self.verbose: print 'Factoring...   (dt = %e)'%dt
-                Ainv = self.Solver(A, **self.solverOpts)
-                if self.verbose: print 'Done'
-
-            rhs = self.getRHS(tInd+1) # this is on the nodes of the time mesh
-            Asubdiag = self.getAsubdiag(tInd)
-
-            if self.verbose: print ('    Solving...   (tInd = %i)')% (tInd+1)
-            sol = Ainv * (rhs - Asubdiag * F[:,self._fieldType+'Solution',tInd]) # taking a step
-
-            if self.verbose: print '    Done...'
-
-            if sol.ndim == 1:
-                sol.shape = (sol.size,1)
-            F[:,self._fieldType+'Solution',tInd+1] = sol
-        if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
-        Ainv.clean()
-        return F
-
-
-    def Jvec(self, m, v, f=None):
-        """
-        Jvec computes the sensitivity times a vector
-
-        .. math::
-            \mathbf{J} \mathbf{v} = \\frac{d\mathbf{P}}{d\mathbf{F}} \left( \\frac{d\mathbf{F}}{d\mathbf{u}} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{\partial\mathbf{F}}{\partial\mathbf{m}} \\right) \mathbf{v}
-
-        where
-
-        .. math::
-            \mathbf{A} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} = \\frac{d \mathbf{RHS}}{d \mathbf{m}}
-        """
-
-        if f is None:
-           f = self.fields(m)
-
-        ftype = self._fieldType + 'Solution' # the thing we solved for
-        self.curModel = m
-
-        # mat to store previous time-step's solution deriv times a vector for each source
-        # size: nu x nSrc
-
-        # this is a bit silly
-
-        # if self._fieldType is 'b' or self._fieldType is 'j':
-        #     ifields = np.zeros((self.mesh.nF, len(Srcs)))
-        # elif self._fieldType is 'e' or self._fieldType is 'h':
-        #     ifields = np.zeros((self.mesh.nE, len(Srcs)))
-
-        # for i, src in enumerate(self.survey.srcList):
-        dun_dm_v = np.hstack([Utils.mkvc(self.getInitialFieldsDeriv(src,v),2) for src in self.survey.srcList]) # can over-write this at each timestep
-        #
-        df_dm_v = Fields_Derivs(self.mesh, self.survey) # store the field derivs we need to project to calc full deriv
-
-        Adiaginv = None
-
-        for tInd, dt in zip(range(self.nT), self.timeSteps):
-            if Adiaginv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same
-                Adiaginv.clean()
-                Adiaginv = None
-
-            if Adiaginv is None:
-                A = self.getAdiag(tInd)
-                Adiaginv = self.Solver(A, **self.solverOpts)
-
-            Asubdiag = self.getAsubdiag(tInd)
-
-            for i, src in enumerate(self.survey.srcList):
-
-                # here, we are lagging by a timestep, so filling in as we go
-                for projField in set([rx.projField for rx in src.rxList]):
-                    # Seogi: df_duFun?
-                    df_dmFun = getattr(f, '_%sDeriv'%projField, None)
-                    # df_dm_v is dense, but we only need the times at (rx.P.T * ones > 0)
-                    # This should be called rx.footprint
-                    df_dm_v[src, '%sDeriv'%projField , tInd] = df_dmFun(tInd, src, dun_dm_v[:,i], v)
-
-                un_src = f[src,ftype,tInd+1]
-
-                dA_dm_v   = self.getAdiagDeriv(tInd, un_src, v) # cell centered on time mesh
-                dRHS_dm_v = self.getRHSDeriv(tInd+1, src, v) # on nodes of time mesh
-
-                dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, f[src,ftype,tInd], v)
-
-                JRHS = dRHS_dm_v - dAsubdiag_dm_v - dA_dm_v
-
-                # step in time and overwrite
-                if tInd != len(self.timeSteps+1):
-                    dun_dm_v[:,i] = Adiaginv * (JRHS - Asubdiag * dun_dm_v[:,i])
-
-        # Seogi: suspcious spot
-        # Jv = self.dataPair(self.survey)
-        Jv = []
-        for src in self.survey.srcList:
-            for rx in src.rxList:
-                # Looping over data class append memory as well!!
-                # Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(df_dm_v[src,'%sDeriv'%rx.projField,:]))
-                Jv.append(rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(df_dm_v[src,'%sDeriv'%rx.projField,:])))
-        Adiaginv.clean()
-        # del df_dm_v, dun_dm_v, Asubdiag
-        # return Utils.mkvc(Jv)
-        return np.hstack(Jv)
-
-    def Jtvec(self, m, v, f=None):
-
-        """
-        Jvec computes the adjoint of the sensitivity times a vector
-
-        .. math::
-            \mathbf{J}^\\top \mathbf{v} =  \left( \\frac{d\mathbf{u}}{d\mathbf{m}} ^ \\top \\frac{d\mathbf{F}}{d\mathbf{u}} ^ \\top  + \\frac{\partial\mathbf{F}}{\partial\mathbf{m}} ^ \\top \\right) \\frac{d\mathbf{P}}{d\mathbf{F}} ^ \\top \mathbf{v}
-
-        where
-
-        .. math::
-            \\frac{d\mathbf{u}}{d\mathbf{m}} ^\\top \mathbf{A}^\\top  + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} ^ \\top = \\frac{d \mathbf{RHS}}{d \mathbf{m}} ^ \\top
-        """
-
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-        ftype = self._fieldType + 'Solution' # the thing we solved for
-
-        # Ensure v is a data object.
-        if not isinstance(v, self.dataPair):
-            v = self.dataPair(self.survey, v)
-
-        df_duT_v = Fields_Derivs(self.mesh, self.survey)
-        ATinv_df_duT_v = np.zeros((len(self.survey.srcList), len(f[self.survey.srcList[0],ftype,0])), dtype=float) # same size as fields at a single timestep
-
-        JTv = np.zeros(m.shape, dtype=float)
-
-        # Loop over sources and receivers to create a fields object: PT_v, df_duT_v, df_dmT_v
-        PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources)
-        for src in self.survey.srcList:
-            # Looping over initializing field class is appending memory!
-            # PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources)
-            # initialize size
-            df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(f[src, self._fieldType, :])
-
-            for rx in src.rxList:
-                print ('_%sDeriv')%(rx.projField)
-                PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) # this is +=
-
-                # PT_v = np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F')
-                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
-
-                for tInd in range(self.nT+1):
-                    cur = df_duTFun(tInd, src, None, Utils.mkvc(PT_v[src,'%sDeriv'%rx.projField,tInd]), adjoint=True)
-                    df_duT_v[src, '%sDeriv'%self._fieldType, tInd] = df_duT_v[src, '%sDeriv'%self._fieldType, tInd] + Utils.mkvc(cur[0],2)
-                    JTv = cur[1] + JTv
-
-        del PT_v # no longer need this
-
-        AdiagTinv = None
-
-        # Do the back-solve through time
-        for tIndP in reversed(range(self.nT + 1)):
-            tInd = tIndP - 1
-            if AdiagTinv is not None and (tInd <= self.nT and self.timeSteps[tInd] != self.timeSteps[tInd+1]): # if the previous timestep is the same --> no need to refactor the matrix
-                AdiagTinv.clean()
-                AdiagTinv = None
-
-            # refactor if we need to
-            if AdiagTinv is None and tInd > -1:
-                Adiag = self.getAdiag(tInd)
-                AdiagTinv = self.Solver(Adiag.T, **self.solverOpts)
-
-            dAsubdiag_dm_v = Zero()
-
-            if tInd < self.nT - 1:
-                Asubdiag = self.getAsubdiag(tInd+1)
-
-
-            for isrc, src in enumerate(self.survey.srcList):
-                # solve against df_duT_v
-                if tInd >= self.nT-1:
-                    # last timestep (first to be solved)
-                    ATinv_df_duT_v[isrc,:] = AdiagTinv * df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]
-                elif tInd > -1:
-                # else:
-                    ATinv_df_duT_v[isrc,:] = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
-                else:
-                    # AdiagTinv = I
-                    ATinv_df_duT_v[isrc,:] = Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])
-                    # - Utils.mkvc(Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
-                    # (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
-
-                if tInd < self.nT - 1:
-                    dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], ATinv_df_duT_v[isrc,:], adjoint = True)
-
-                if tInd > -1:
-                    un_src = f[src,ftype,tInd+1]
-                    dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
-                    dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
-
-                    JTv = JTv + Utils.mkvc(- dAT_dm_v - dAsubdiag_dm_v + dRHST_dm_v)
-                else:
-                    # dA_dm_v = self.getInitialFieldsDeriv(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True)
-                    # print np.linalg.norm(self.getInitialFieldsDeriv(src, df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True))
-                    # print np.linalg.norm(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
-                    # vec = - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]) + Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
-                    # dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint = True)
-                    dRHST_dm_v = Utils.mkvc(self.getInitialFieldsDeriv(src, Utils.mkvc(ATinv_df_duT_v[isrc,:]) , adjoint=True))
-
-                    JTv = JTv + Utils.mkvc( -dAsubdiagT_dm_v + dRHST_dm_v) #
-
-
-
-                    # # dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
-                    # dRHST_dm_v0 = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
-                    # dRHST_dm_v1 = self.getInitialFieldsDeriv( Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
-                    # JTv = JTv + Utils.mkvc(dRHST_dm_v0 + dRHST_dm_v1)
-
-                    # print 'here'
-                    # inFields = self.getInitialFieldsDeriv(f[src,ftype,tInd+1], Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
-                    #  # - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint=True)
-                    # print inFields.shape
-                    # JTv = JTv + inFields
-                # dAsubdiag_dm_v = 0
-
-
-
-        # Missing the 0 step
-
-        # adding du_dm^T * dF_du^T * P^T vfor time 0 (no dRHS_dm_v at time 0)
-        # Asubdiag = self.getAsubdiag(0)
-        # for src in self.survey.srcList:
-        #     for projField in set(rx.projField):
-        #         v = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
-        #         JTv = JTv - Utils.mkvc(self.getAdiagDeriv(0, f[src, ftype, tInd], v, adjoint = True))
-        #         # JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0] - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])), adjoint=True)
-
-        # del df_duT_v, ATinv_df_duT_v, A, Asubdiag
-        if AdiagTinv is not None:
-            AdiagTinv.clean()
-
-        return Utils.mkvc(JTv).astype(float)
-
-
-
-    def getSourceTerm(self, tInd):
-
-        Srcs = self.survey.srcList
-
-        if self._eqLocs is 'FE':
-            S_m = np.zeros((self.mesh.nF,len(Srcs)))
-            S_e = np.zeros((self.mesh.nE,len(Srcs)))
-        elif self._eqLocs is 'EF':
-            S_m = np.zeros((self.mesh.nE,len(Srcs)))
-            S_e = np.zeros((self.mesh.nF,len(Srcs)))
-
-        for i, src in enumerate(Srcs):
-            smi, sei = src.eval(self, self.times[tInd])
-            S_m[:,i] = S_m[:,i] + smi
-            S_e[:,i] = S_e[:,i] + sei
-
-        return S_m, S_e
-
-    def getInitialFields(self):
-
-        Srcs = self.survey.srcList
-
-        if self._fieldType is 'b' or self._fieldType is 'j':
-            ifields = np.zeros((self.mesh.nF, len(Srcs)))
-        elif self._fieldType is 'e' or self._fieldType is 'h':
-            ifields = np.zeros((self.mesh.nE, len(Srcs)))
-
-        for i,src in enumerate(Srcs):
-            ifields[:,i] = ifields[:,i] + getattr(src, '%sInitial'%self._fieldType, None)(self)
-
-        return ifields
-
-    def getInitialFieldsDeriv(self, src, v, adjoint=False):
-
-        if adjoint is False:
-            if self._fieldType is 'b' or self._fieldType is 'j':
-                ifieldsDeriv = np.zeros(self.mesh.nF)
-            elif self._fieldType is 'e' or self._fieldType is 'h':
-                ifieldsDeriv = np.zeros(self.mesh.nE)
-
-        elif adjoint is True:
-            ifieldsDeriv = np.zeros(self.mapping.nP)
-
-        ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv
-
-            # ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv
-            # ifieldsDeriv = self.getAdiagDeriv(None, u, v, adjoint)
-            # ifieldsDeriv = ifieldsDeriv.sum()
-
-        return ifieldsDeriv
-
-
-##########################################################################################
-################################ E-B Formulation #########################################
-##########################################################################################
-
-# ------------------------------- Problem_b -------------------------------------------- #
-
-class Problem_b(BaseTDEMProblem):
-    """
-    Starting from the quasi-static E-B formulation of Maxwell's equations (semi-discretized)
-
-    .. math::
-
-        \mathbf{C} \mathbf{e} + \\frac{\partial \mathbf{b}}{\partial t} = \mathbf{s_m} \\\\
-        \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}
-
-    where :math:`\mathbf{s_e}` is an integrated quantity, we eliminate :math:`\mathbf{e}` using
-
-    .. math::
-        \mathbf{e} = \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}
-
-    to obtain a second order semi-discretized system in :math:`\mathbf{b}`
-
-    .. math::
-        \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b}  + \\frac{\partial \mathbf{b}}{\partial t} = \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e} + \mathbf{s_m}
-
-    and moving everything except the time derivative to the rhs gives
-
-    .. math::
-        \\frac{\partial \mathbf{b}}{\partial t} = -\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} + \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e} + \mathbf{s_m}
-
-    For the time discretization, we use backward euler. To solve for the :math:`n+1`th time step, we have
-
-    .. math::
-        \\frac{\mathbf{b}^{n+1} - \mathbf{b}^{n}}{\mathbf{dt}} = -\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b}^{n+1} + \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}^{n+1} + \mathbf{s_m}^{n+1}
-
-    re-arranging to put :math:`\mathbf{b}^{n+1}` on the left hand side gives
-
-    .. math::
-        (\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}) \mathbf{b}^{n+1} = \mathbf{b}^{n} + \mathbf{dt}(\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}^{n+1} + \mathbf{s_m}^{n+1})
-
-    :param Mesh mesh: mesh
-    :param Mapping mapping: mapping
-    """
-
-    _fieldType = 'b'
-    _eqLocs    = 'FE'
-    fieldsPair = Fields_b
-    surveyPair = SurveyTDEM
-
-    def __init__(self, mesh, mapping=None, **kwargs):
-        BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs)
-
-    def getAdiag(self, tInd):
-        """
-        System matrix at a given time index
-
-        .. math::
-            (\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f})
-
-        """
-        assert tInd >= 0 and tInd < self.nT
-
-        dt = self.timeSteps[tInd]
-        C = self.mesh.edgeCurl
-        MeSigmaI = self.MeSigmaI
-        MfMui = self.MfMui
-        I = Utils.speye(self.mesh.nF)
-
-        A = 1./dt * I + ( C * ( MeSigmaI * (C.T * MfMui ) ) )
-
-        if self._makeASymmetric is True:
-            return MfMui.T * A
-        return A
-
-    def getAdiagDeriv(self, tInd, u, v, adjoint=False):
-        C = self.mesh.edgeCurl
-        MeSigmaIDeriv = lambda x: self.MeSigmaIDeriv(x)
-        MfMui = self.MfMui
-
-        if adjoint:
-            if self._makeASymmetric is True:
-                v = MfMui * v
-            return  MeSigmaIDeriv(C.T * ( MfMui * u )).T * ( C.T * v )
-
-        ADeriv = ( C * ( MeSigmaIDeriv(C.T * ( MfMui * u )) * v ) )
-        if self._makeASymmetric is True:
-            return MfMui.T * ADeriv
-        return ADeriv
-
-
-    def getAsubdiag(self, tInd):
-
-        dt = self.timeSteps[tInd]
-        MfMui = self.MfMui
-        Asubdiag = - 1./dt * sp.eye(self.mesh.nF)
-
-        if self._makeASymmetric is True:
-            return MfMui.T * Asubdiag
-
-        return Asubdiag
-
-    def getAsubdiagDeriv(self, tInd, u, v, adjoint=False):
-        return Zero() * v
-
-
-
-    def getRHS(self, tInd):
-        C = self.mesh.edgeCurl
-        MeSigmaI = self.MeSigmaI
-        MfMui = self.MfMui
-
-        S_m, S_e = self.getSourceTerm(tInd)
-
-        rhs =  (C * (MeSigmaI * S_e) + S_m)
-        if self._makeASymmetric is True:
-            return MfMui.T * rhs
-        return rhs
-
-    def getRHSDeriv(self, tInd, src, v, adjoint=False):
-
-        C = self.mesh.edgeCurl
-        MeSigmaI = self.MeSigmaI
-        MeSigmaIDeriv = lambda u: self.MeSigmaIDeriv(u)
-        MfMui = self.MfMui
-
-        _, S_e = src.eval(tInd, self)
-        S_mDeriv, S_eDeriv = src.evalDeriv(self.times[tInd], self, adjoint=adjoint)
-
-        if adjoint:
-            if self._makeASymmetric is True:
-                v = self.MfMui * v
-            if isinstance(S_e, Utils.Zero):
-                MeSigmaIDerivT_v = Utils.Zero()
-            else:
-                MeSigmaIDerivT_v = MeSigmaIDeriv(S_e).T * v
-            RHSDeriv = MeSigmaIDerivT_v + S_eDeriv( MeSigmaI.T *  ( C.T * v ) ) + S_mDeriv(v)
-            return RHSDeriv
-
-        if isinstance(S_e, Utils.Zero):
-            MeSigmaIDeriv_v = Utils.Zero()
-        else:
-            MeSigmaIDeriv_v = MeSigmaIDeriv(S_e) * v
-
-        RHSDeriv = (C * (MeSigmaIDeriv_v + MeSigmaI * S_eDeriv(v) + S_mDeriv(v)))
-
-        if self._makeASymmetric is True:
-            return self.MfMui.T * RHSDeriv
-        return RHSDeriv
-
-
-# ------------------------------- Problem_e -------------------------------------------- #
-
-class Problem_e(BaseTDEMProblem):
-
-    _fieldType = 'e'
-    _eqLocs    = 'FE'
-    fieldsPair = Fields_e
-    surveyPair = SurveyTDEM
-
-    def __init__(self, mesh, mapping=None, **kwargs):
-        BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs)
-
-    def getAdiag(self, tInd):
-        """
-        System matrix at a given time index
-
-        """
-        assert tInd >= 0 and tInd < self.nT
-
-        dt = self.timeSteps[tInd]
-        C = self.mesh.edgeCurl
-        MfMui = self.MfMui
-        MeSigma = self.MeSigma
-
-        return C.T * ( MfMui * C ) + 1./dt * MeSigma
-
-
-    def getAdiagDeriv(self, tInd, u, v, adjoint=False):
-        assert tInd >= 0 and tInd < self.nT
-
-        dt = self.timeSteps[tInd]
-        C = self.mesh.edgeCurl
-        MfMui = self.MfMui
-        MeSigmaDeriv = self.MeSigmaDeriv(u)
-
-        if adjoint:
-            return 1./dt * MeSigmaDeriv.T * v
-
-        return 1./dt * MeSigmaDeriv * v
-
-
-    def getAsubdiag(self, tInd):
-        assert tInd >= 0 and tInd < self.nT
-
-        dt = self.timeSteps[tInd]
-
-        return - 1./dt * self.MeSigma
-
-    def getAsubdiagDeriv(self, tInd, u, v, adjoint=False):
-        dt = self.timeSteps[tInd]
-
-        if adjoint:
-            return - 1./dt * self.MeSigmaDeriv(u).T * v
-
-        return - 1./dt * self.MeSigmaDeriv(u) * v
-
-    def getRHS(self, tInd):
-        return Zero()
-
-    def getRHSDeriv(self, tInd, src, v, adjoint=False):
-        return Zero()
-
-
-
-
-

SimPEG/EM/TDEM/__init__.py

@@ -1,3 +1,3 @@
-from TDEM import BaseTDEMProblem, Problem_b, Problem_e
-from FieldsTDEM import Fields, Fields_b
-from SurveyTDEM import Survey, Src, Rx
+from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
+from BaseTDEM import BaseTDEMProblem, FieldsTDEM
+from TDEM_b import ProblemTDEM_b

SimPEG/EM/TDEM_old/SurveyTDEM.py

@@ -1,199 +0,0 @@
-from SimPEG import Utils, Survey, np
-from SimPEG.Survey import BaseSurvey
-from SimPEG.EM.Utils import *
-from BaseTDEM import FieldsTDEM
-import SrcTDEM as Src
-
-class RxTDEM(Survey.BaseTimeRx):
-
-    knownRxTypes = {
-                    'ex':['e', 'Ex', 'N'],
-                    'ey':['e', 'Ey', 'N'],
-                    'ez':['e', 'Ez', 'N'],
-
-                    'bx':['b', 'Fx', 'N'],
-                    'by':['b', 'Fy', 'N'],
-                    'bz':['b', 'Fz', 'N'],
-
-                    'dbxdt':['b', 'Fx', 'CC'],
-                    'dbydt':['b', 'Fy', 'CC'],
-                    'dbzdt':['b', 'Fz', 'CC'],
-                   }
-
-    def __init__(self, locs, times, rxType):
-        Survey.BaseTimeRx.__init__(self, locs, times, rxType)
-
-    @property
-    def projField(self):
-        """Field Type projection (e.g. e b ...)"""
-        return self.knownRxTypes[self.rxType][0]
-
-    @property
-    def projGLoc(self):
-        """Grid Location projection (e.g. Ex Fy ...)"""
-        return self.knownRxTypes[self.rxType][1]
-
-    @property
-    def projTLoc(self):
-        """Time Location projection (e.g. CC N)"""
-        return self.knownRxTypes[self.rxType][2]
-
-    def getTimeP(self, timeMesh):
-        """
-            Returns the time projection matrix.
-
-            .. note::
-
-                This is not stored in memory, but is created on demand.
-        """
-        if self.rxType in ['dbxdt','dbydt','dbzdt']:
-            return timeMesh.getInterpolationMat(self.times, self.projTLoc)*timeMesh.faceDiv
-        else:
-            return timeMesh.getInterpolationMat(self.times, self.projTLoc)
-
-    def eval(self, src, mesh, timeMesh, u):
-        P = self.getP(mesh, timeMesh)
-        u_part = Utils.mkvc(u[src, self.projField, :])
-        return P*u_part
-
-    def evalDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
-        P = self.getP(mesh, timeMesh)
-
-        if not adjoint:
-            return P * Utils.mkvc(v[src, self.projField, :])
-        elif adjoint:
-            return P.T * v[src, self]
-
-
-class SrcTDEM(Survey.BaseSrc):
-    rxPair = RxTDEM
-    radius = None
-
-    def getInitialFields(self, mesh):
-        F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
-        return F0
-
-    def getJs(self, mesh, time):
-        return None
-
-
-class SrcTDEM_VMD_MVP(SrcTDEM):
-
-    def __init__(self,rxList,loc,waveformType="STEPOFF"):
-        self.loc = loc
-        self.waveformType = waveformType
-        SrcTDEM.__init__(self,rxList)
-
-    def getInitialFields(self, mesh):
-        """Vertical magnetic dipole, magnetic vector potential"""
-        if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"
-            if mesh._meshType is 'CYL':
-                if mesh.isSymmetric:
-                    MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
-                else:
-                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            elif mesh._meshType is 'TENSOR':
-                MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
-            else:
-                raise Exception('Unknown mesh for VMD')
-            return {"b": mesh.edgeCurl*MVP}
-        elif self.waveformType == "GENERAL":
-            print ">> General waveform: Zero initial condition"
-            return {"b": np.zeros(mesh.nF)}
-        else:
-            raise NotImplementedError("Only use STEPOFF or GENERAL")
-
-    def getMeS(self, mesh, MfMui):
-        if mesh._meshType is 'CYL':
-            if mesh.isSymmetric:
-                MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
-            else:
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-        elif mesh._meshType is 'TENSOR':
-            MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
-        else:
-            raise Exception('Unknown mesh for VMD')
-        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
-
-
-class SrcTDEM_CircularLoop_MVP(SrcTDEM):
-    def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
-        self.loc = loc
-        self.radius = radius
-        self.waveformType = waveformType
-        SrcTDEM.__init__(self,rxList)
-
-    def getInitialFields(self, mesh):
-        """Circular Loop, magnetic vector potential"""
-        if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"
-            if mesh._meshType is 'CYL':
-                if mesh.isSymmetric:
-                    MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
-                else:
-                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            elif mesh._meshType is 'TENSOR':
-                MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
-            else:
-                raise Exception('Unknown mesh for CircularLoop')
-            return {"b": mesh.edgeCurl*MVP}
-        elif self.waveformType == "GENERAL":
-            print ">> General waveform: Zero initial condition"
-            return {"b": np.zeros(mesh.nF)}
-        else:
-            raise NotImplementedError("Only use STEPOFF or GENERAL")
-
-    def getMeS(self, mesh, MfMui):
-        if mesh._meshType is 'CYL':
-            if mesh.isSymmetric:
-                MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
-            else:
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-        elif mesh._meshType is 'TENSOR':
-            MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
-        else:
-            raise Exception('Unknown mesh for CircularLoop')
-        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
-
-
-class SurveyTDEM(Survey.BaseSurvey):
-    """
-        docstring for SurveyTDEM
-    """
-    srcPair = SrcTDEM
-
-    def __init__(self, srcList, **kwargs):
-        # Sort these by frequency
-        self.srcList = srcList
-        Survey.BaseSurvey.__init__(self, **kwargs)
-
-    def projectFields(self, u):
-        data = Survey.Data(self)
-        for src in self.srcList:
-            for rx in src.rxList:
-                data[src, rx] = rx.projectFields(src, self.mesh, self.prob.timeMesh, u)
-        return data
-
-    def projectFieldsDeriv(self, u, v=None, adjoint=False):
-        assert v is not None, 'v to multiply must be provided.'
-
-        if not adjoint:
-            data = Survey.Data(self)
-            for src in self.srcList:
-                for rx in src.rxList:
-                    data[src, rx] = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v)
-            return data
-        else:
-            f = FieldsTDEM(self.mesh, self)
-            for src in self.srcList:
-                for rx in src.rxList:
-                    Ptv = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
-                    Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
-                    if rx.projField not in f: # first time we are projecting
-                        f[src, rx.projField, :] = Ptv
-                    else: # there are already fields, so let's add to them!
-                        f[src, rx.projField, :] += Ptv
-            return f
-
-

SimPEG/EM/TDEM_old/__init__.py

@@ -1,3 +0,0 @@
-from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
-from BaseTDEM import BaseTDEMProblem, FieldsTDEM
-from TDEM_b import ProblemTDEM_b

SimPEG/Examples/EM_TDEM_1D_Inversion.py

@@ -42,10 +42,10 @@ def run(plotIt=True):
 
 
     rxOffset=1e-3
-    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
-    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 80]))
-    survey = EM.TDEM.Survey([src])
-    prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+    rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
+    src = EM.TDEM.SrcTDEM_VMD_MVP([rx], np.array([0., 0., 80]))
+    survey = EM.TDEM.SurveyTDEM([src])
+    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
 
     prb.Solver = SolverLU
     prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)]
@@ -53,9 +53,9 @@ def run(plotIt=True):
 
     # create observed data
     std = 0.05
-
+    
     survey.dobs = survey.makeSyntheticData(mtrue,std)
-    survey.std = std
+    survey.std = std 
     survey.eps = 1e-5*np.linalg.norm(survey.dobs)
 
     if plotIt:

SimPEG/Examples/Forward_BasicDirectCurrent.py

@@ -1,25 +1,22 @@
 from SimPEG import Mesh, Utils, np, SolverLU
 
+## 2D DC forward modeling example with Tensor and Curvilinear Meshes
+
 def run(plotIt=True):
-
-    """
-        Mesh: Basic Forward 2D DC Resistivity
-        =====================================
-
-        2D DC forward modeling example with Tensor and Curvilinear Meshes
-    """
-
     # Step1: Generate Tensor and Curvilinear Mesh
     sz = [40,40]
+    # Tensor Mesh
     tM = Mesh.TensorMesh(sz)
+    # Curvilinear Mesh
     rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
-
     # Step2: Direct Current (DC) operator
     def DCfun(mesh, pts):
         D = mesh.faceDiv
+        G = D.T
         sigma = 1e-2*np.ones(mesh.nC)
-        MsigI = mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True)
-        A = -D*MsigI*D.T
+        Msigi = mesh.getFaceInnerProduct(1./sigma)
+        MsigI = Utils.sdInv(Msigi)
+        A = D*MsigI*G
         A[-1,-1] /= mesh.vol[-1] # Remove null space
         rhs = np.zeros(mesh.nC)
         txind = Utils.meshutils.closestPoints(mesh, pts)
@@ -40,17 +37,39 @@ def run(plotIt=True):
     if not plotIt: return
 
     import matplotlib.pyplot as plt
+    import matplotlib
+    from matplotlib.mlab import griddata
 
     #Step4: Making Figure
     fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
+    label = ["(a)", "(b)"]
+    opts = {}
     vmin, vmax = phitM.min(), phitM.max()
     dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True)
-    dat = rM.plotImage(phirM, ax=axes[1], clim=(vmin, vmax), grid=True)
+
+    #TODO: At the moment Curvilinear Mesh do not have plotimage
+
+    Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F')
+    Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F')
+    PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear')
+    axes[1].contourf(Xi, Yi, PHIrM, 100, vmin=vmin, vmax=vmax)
+
     cb = plt.colorbar(dat[0], ax=axes[0]); cb.set_label("Voltage (V)")
     cb = plt.colorbar(dat[0], ax=axes[1]); cb.set_label("Voltage (V)")
 
+    tM.plotGrid(ax=axes[0], **opts)
     axes[0].set_title('TensorMesh')
+    rM.plotGrid(ax=axes[1], **opts)
     axes[1].set_title('CurvilinearMesh')
+    for i in range(2):
+        axes[i].set_xlim(0.025, 0.975)
+        axes[i].set_ylim(0.025, 0.975)
+        axes[i].text(0., 1.0, label[i], fontsize=20)
+        if i==0:
+            axes[i].set_ylabel("y")
+        else:
+            axes[i].set_ylabel(" ")
+        axes[i].set_xlabel("x")
     plt.show()
 
 

SimPEG/Examples/Inversion_IRLS.py

@@ -1,7 +1,7 @@
 from SimPEG import *
 
 
-def run(N=100, plotIt=True):
+def run(N=200, plotIt=True):
     """
         Inversion: Linear Problem
         =========================
@@ -18,8 +18,6 @@ def run(N=100, plotIt=True):
     mesh = Mesh.TensorMesh([N])
 
     m0 = np.ones(mesh.nC) * 1e-4
-    mref = np.zeros(mesh.nC)
-
     nk = 10
     jk = np.linspace(1.,nk,nk)
     p = -2.
@@ -52,47 +50,57 @@ def run(N=100, plotIt=True):
     wr = np.sum(prob.G**2.,axis=0)**0.5
     wr = ( wr/np.max(wr) )
 
-#    reg = Regularization.Simple(mesh)
-#    reg.mref = mref
-#    reg.cell_weights = wr
-#
+    reg = Regularization.Simple(mesh)
+    reg.wght = wr
+
     dmis = DataMisfit.l2_DataMisfit(survey)
     dmis.Wd = 1./wd
-#
-#    opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4)
-#    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
-#    invProb.curModel = m0
-#
-#    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
-#    target = Directives.TargetMisfit()
-#
+
+    opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4)
+    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
+    invProb.curModel = m0
+
+    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
+    target = Directives.TargetMisfit()
+
     betaest = Directives.BetaEstimate_ByEig()
-#    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
-#
-#
-#    mrec = inv.run(m0)
-#    ml2 = mrec
-#    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
-#
-#    # Switch regularization to sparse
-#    phim = invProb.phi_m_last
-#    phid =  invProb.phi_d
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
+
+
+    mrec = inv.run(m0)
+    ml2 = mrec
+    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
+
+    # Switch regularization to sparse
+    phim = invProb.phi_m_last
+    phid =  invProb.phi_d
 
     reg = Regularization.Sparse(mesh)
-    reg.mref = mref
-    reg.cell_weights = wr
 
+#==============================================================================
+#     fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
+#     dmdx = reg.mesh.cellDiffxStencil * mrec
+#     plt.plot(np.sort(dmdx))
+#==============================================================================
+
+    #reg.recModel = mrec
+    reg.wght = np.ones(mesh.nC)
     reg.mref = np.zeros(mesh.nC)
-    eps_p = 5e-2
-    eps_q = 5e-2
-    norms   = [0., 0., 2., 2.]
+    reg.eps_p = 5e-2
+    reg.eps_q = 1e-2
+    reg.norms   = [0., 0., 2., 2.]
+    reg.wght = wr
 
-    opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
-    update_Jacobi = Directives.Update_lin_PreCond()
-    IRLS = Directives.Update_IRLS( norms=norms,  eps_p=eps_p, eps_q=eps_q)
+    opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3)
+    invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.)
+    beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
+    #betaest = Directives.BetaEstimate_ByEig()
+    target = Directives.TargetMisfit()
+    IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
 
-    inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
+
+    m0 = mrec
 
     # Run inversion
     mrec = inv.run(m0)
@@ -109,7 +117,7 @@ def run(N=100, plotIt=True):
         axes[0].set_title('Columns of matrix G')
 
         axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
-        axes[1].plot(mesh.vectorCCx, reg.l2model, 'r-')
+        axes[1].plot(mesh.vectorCCx, ml2, 'r-')
         #axes[1].legend(('True Model', 'Recovered Model'))
         axes[1].set_ylim(-1.0,1.25)
 

SimPEG/Examples/__init__.py

@@ -8,9 +8,9 @@ import EM_FDEM_Analytic_MagDipoleWholespace
 import EM_Schenkel_Morrison_Casing
 import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
+import Forward_BasicDirectCurrent
 import Inversion_IRLS
 import Inversion_Linear
-import Mesh_Basic_ForwardDC
 import Mesh_Basic_PlotImage
 import Mesh_Basic_Types
 import Mesh_Operators_CahnHilliard
@@ -22,7 +22,7 @@ import MT_1D_ForwardAndInversion
 import MT_3D_Foward
 import Utils_surface2ind_topo
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Examples/simpegEMpaper_example_1.py

@@ -1,157 +0,0 @@
-from SimPEG import np, Mesh, Maps, Utils, DataMisfit, Regularization, Optimization, Inversion, InvProblem, Directives
-from SimPEG import SolverLU
-from SimPEG.EM import FDEM, TDEM, mu_0
-import matplotlib.pyplot as plt
-import matplotlib
-matplotlib.rcParams['font.size'] = 14
-
-def run(plotIt=True):
-    # Set up cylindrically symmeric mesh
-    cs, ncx, ncz, npad = 10., 15, 25, 13  # padded cyl mesh
-    hx = [(cs,ncx), (cs,npad,1.3)]
-    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
-    mesh = Mesh.CylMesh([hx,1,hz], '00C')
-
-    # Conductivity model
-    layerz = np.r_[-200., -100.]
-    layer = (mesh.vectorCCz>=layerz[0]) & (mesh.vectorCCz<=layerz[1])
-    active = mesh.vectorCCz<0.
-    sig_half  = 1e-2  # Half-space conductivity
-    sig_air   = 1e-8  # Air conductivity
-    sig_layer = 5e-2  # Layer conductivity
-    sigma = np.ones(mesh.nCz)*sig_air
-    sigma[active] = sig_half
-    sigma[layer] = sig_layer
-
-    # Mapping
-    actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
-    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
-    mtrue = np.log(sigma[active])
-
-    # FDEM problem & survey
-    rxlocs = Utils.ndgrid([np.r_[50.], np.r_[0], np.r_[0.]])
-    bzi = FDEM.Rx.Point_bSecondary(rxlocs, 'z', 'real')
-    bzr = FDEM.Rx.Point_bSecondary(rxlocs, 'z', 'imag')
-
-    freqs = np.logspace(2, 3, 5)
-    srcLoc = np.array([0., 0., 0.])
-
-    print 'min skin depth = ', 500./np.sqrt(freqs.max() * sig_half), 'max skin depth = ', 500./np.sqrt(freqs.min() * sig_half)
-    print 'max x ', mesh.vectorCCx.max(), 'min z ', mesh.vectorCCz.min(), 'max z ', mesh.vectorCCz.max()
-
-    srcList = []
-    [srcList.append(FDEM.Src.MagDipole([bzr, bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
-
-    surveyFD = FDEM.Survey(srcList)
-    prbFD = FDEM.Problem3D_b(mesh, mapping=mapping)
-    prbFD.pair(surveyFD)
-    std = 0.03
-    surveyFD.makeSyntheticData(mtrue, std)
-    surveyFD.eps = np.linalg.norm(surveyFD.dtrue)*1e-5
-
-    # FDEM inversion
-    np.random.seed(1)
-    dmisfit = DataMisfit.l2_DataMisfit(surveyFD)
-    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
-    reg = Regularization.Simple(regMesh)
-    opt = Optimization.InexactGaussNewton(maxIterCG=10, maxIter=4)
-    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
-    # Inversion Directives
-    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=3)
-    # betaest = Directives.BetaEstimate_ByEig(beta0_ratio=10.)
-    invProb.beta = 1.
-    target = Directives.TargetMisfit()
-
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta,target])
-    m0 = np.log(np.ones(mtrue.size)*sig_half)
-    reg.alpha_s = 5e-1
-    reg.alpha_x = 1.
-    prbFD.counter = opt.counter = Utils.Counter()
-    opt.LSshorten = 0.5
-    opt.tolG = 1e-10
-    opt.eps = 1e-10
-    opt.remember('xc')
-    moptFD = inv.run(m0)
-
-    # TDEM problem
-    times = np.logspace(-4, np.log10(2e-3), 10)
-    print 'min diffusion distance ', 1.28*np.sqrt(times.min()/(sig_half*mu_0)), 'max diffusion distance ', 1.28*np.sqrt(times.max()/(sig_half*mu_0))
-    rx = TDEM.Rx(rxlocs, times, 'bz')
-    src = TDEM.Src.MagDipole([rx], waveform=TDEM.Src.StepOffWaveform(), loc=srcLoc) # same src location as FDEM problem
-
-    surveyTD = TDEM.Survey([src])
-    prbTD = TDEM.Problem_b(mesh, mapping=mapping)
-    prbTD.timeSteps = [(5e-5, 10),(1e-4, 10),(5e-4, 10)]
-    prbTD.pair(surveyTD)
-    prbTD.Solver = SolverLU
-
-    std = 0.03
-    surveyTD.makeSyntheticData(mtrue, std)
-    surveyTD.std = std
-    surveyTD.eps = np.linalg.norm(surveyTD.dtrue)*1e-5
-
-    # TDEM inversion
-    dmisfit = DataMisfit.l2_DataMisfit(surveyTD)
-    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
-    reg = Regularization.Simple(regMesh)
-    opt = Optimization.InexactGaussNewton(maxIterCG=10, maxIter=4)
-    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
-
-    # Inversion Directives
-    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=3)
-    invProb.beta = 1.
-    # betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1.)
-    target = Directives.TargetMisfit()
-
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta, target])
-    m0 = np.log(np.ones(mtrue.size)*sig_half)
-    reg.alpha_s = 5e-1
-    reg.alpha_x = 1.
-    prbTD.counter = opt.counter = Utils.Counter()
-    opt.LSshorten = 0.5
-    opt.remember('xc')
-    moptTD = inv.run(m0)
-
-    if plotIt:
-        fig, ax = plt.subplots(1,1, figsize = (4, 6))
-        plt.semilogx(sigma[active], mesh.vectorCCz[active], 'k-', lw=2)
-        plt.semilogx(np.exp(moptFD), mesh.vectorCCz[active], 'ko', ms=3)
-        plt.semilogx(np.exp(moptTD), mesh.vectorCCz[active], 'k*')
-        ax.set_ylim(-1000, 0)
-        ax.set_xlim(5e-3, 1e-1)
-
-        ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
-        ax.set_ylabel('Depth (m)', fontsize = 14)
-        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
-        plt.legend(['True', 'Pred (FD)', 'Pred (TD)'], fontsize=13, loc=4)
-        plt.show()
-
-        fig = plt.figure(figsize = (10*1.3, 5*1.3))
-        ax2 = plt.subplot(122)
-        ax2.plot(times, surveyTD.dobs, 'k-', lw=2)
-        ax2.plot(times, surveyTD.dpred(moptTD), 'ko', ms=4)
-        ax2.set_xscale('log')
-        ax2.set_yscale('log')
-        ax2.set_xlim(times.min(), times.max())
-        ax1 = plt.subplot(121)
-        ax1.plot(freqs, -surveyFD.dobs[::2], 'k-', lw=2)
-        ax1.plot(freqs, -surveyFD.dobs[1::2], 'k--', lw=2)
-        dpredFD = surveyFD.dpred(moptTD)
-        ax1.plot(freqs, -dpredFD[::2], 'ko', ms=4)
-        ax1.plot(freqs, -dpredFD[1::2], 'k+', markeredgewidth=2., ms=10)
-        ax1.set_xscale('log')
-        ax1.set_yscale('log')
-        ax2.set_xlabel('Time (s)', fontsize = 14)
-        ax1.set_xlabel('Frequency (Hz)', fontsize = 14)
-        ax1.set_ylabel('Vertical magnetic field (T)', fontsize = 14)
-        ax2.grid(True,which='minor')
-        ax1.grid(True,which='minor')
-        ax2.set_title("(b) TD observed vs. predicted", fontsize = 14)
-        ax1.set_title("(a) FD observed vs. predicted", fontsize = 14)
-        ax2.legend(("Obs", "Pred"), fontsize = 12)
-        ax1.legend(("Obs", "Pred (real)", "Pred (imag)"), fontsize = 12, loc=3)
-        ax1.set_xlim(freqs.max(), freqs.min())
-        plt.show()
-
-if __name__ == '__main__':
-    run()

SimPEG/Maps.py

@@ -1,3 +1,4 @@
+from __future__ import division
 import Utils, numpy as np, scipy.sparse as sp
 from scipy.sparse.linalg import LinearOperator
 from Tests import checkDerivative
@@ -5,6 +6,7 @@ from PropMaps import PropMap, Property
 from numpy.polynomial import polynomial
 from scipy.interpolate import UnivariateSpline
 import warnings
+from SimPEG.Utils import Zero
 
 class IdentityMap(object):
     """
@@ -17,7 +19,7 @@ class IdentityMap(object):
         Utils.setKwargs(self, **kwargs)
 
         if nP is not None:
-            assert type(nP) in [int, long], ' Number of parameters must be an integer.'
+            assert type(nP) in [int, long, np.int64], ' Number of parameters must be an integer.'
 
         self.mesh = mesh
         self._nP  = nP
@@ -129,7 +131,15 @@ class IdentityMap(object):
 
 
 class ComboMap(IdentityMap):
-    """Combination of various maps."""
+    """
+        Combination of various maps.
+
+        The ComboMap holds the information for multiplying and combining
+        maps. It also uses the chain rule to create the derivative.
+        Remember, any time that you make your own combination of mappings
+        be sure to test that the derivative is correct.
+
+    """
 
     def __init__(self, maps, **kwargs):
         IdentityMap.__init__(self, None, **kwargs)
@@ -178,6 +188,12 @@ class ComboMap(IdentityMap):
 
 class ExpMap(IdentityMap):
     """
+        Electrical conductivity varies over many orders of magnitude, so it is a common
+        technique when solving the inverse problem to parameterize and optimize in terms
+        of log conductivity. This makes sense not only because it ensures all conductivities
+        will be positive, but because this is fundamentally the space where conductivity
+        lives (i.e. it varies logarithmically).
+
         Changes the model into the physical property.
 
         A common example of this is to invert for electrical conductivity
@@ -449,6 +465,32 @@ class Mesh2Mesh(IdentityMap):
     """
         Takes a model on one mesh are translates it to another mesh.
 
+        .. plot::
+
+            from SimPEG import *
+            import matplotlib.pyplot as plt
+            M = Mesh.TensorMesh([100,100])
+            h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
+            h1 = h1/h1.sum()
+            M2 = Mesh.TensorMesh([h1,h1])
+            V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
+            v = Utils.mkvc(V)
+            modh = Maps.Mesh2Mesh([M,M2])
+            modH = Maps.Mesh2Mesh([M2,M])
+            H = modH * v
+            h = modh * H
+            ax = plt.subplot(131)
+            M.plotImage(v, ax=ax)
+            ax.set_title('Fine Mesh (Original)')
+            ax = plt.subplot(132)
+            M2.plotImage(H,clim=[0,1],ax=ax)
+            ax.set_title('Course Mesh')
+            ax = plt.subplot(133)
+            M.plotImage(h,clim=[0,1],ax=ax)
+            ax.set_title('Fine Mesh (Interpolated)')
+            plt.show()
+
+
     """
 
     def __init__(self, meshes, **kwargs):
@@ -501,11 +543,19 @@ class InjectActiveCells(IdentityMap):
         self.indInactive = np.logical_not(indActive)
         if Utils.isScalar(valInactive):
             self.valInactive = np.ones(self.nC)*float(valInactive)
+            self.valInactive[self.indActive] = 0.
         else:
-            self.valInactive = valInactive.copy()
-        self.valInactive[self.indActive] = 0
+            if len(valInactive) == sum(self.indInactive):
+                self.valInactive = np.zeros(nC)
+                self.valInactive[self.indInactive] = valInactive.copy()
+            else:
+                assert len(self.valInactive) == self.nC, 'valInactive must be the size of nC or nInactive'
+                self.valInactive = valInactive.copy()
+                if any(self.valInactive[self.indActive] != 0.):
+                    warnings.warn('the inactive has non-zero values in the active set.')
 
         inds = np.nonzero(self.indActive)[0]
+        # inds[self.indActive]
         self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
 
     @property
@@ -574,6 +624,37 @@ class Weighting(IdentityMap):
     def deriv(self, m):
         return self.P
 
+class Projection(IdentityMap):
+    """
+        A map to rearrange parameters
+
+    """
+
+
+    def __init__(self, indTo, indFrom, shape, mesh=None, **kwargs):
+
+        assert len(indTo) == len(indFrom)
+
+        self.P = sp.csr_matrix((np.ones(len(indTo)), (indTo, indFrom)), shape=shape)
+        self._shape = shape
+
+        super(Projection, self).__init__(mesh, **kwargs)
+
+    @property
+    def shape(self):
+        return self._shape
+
+    @property
+    def nP(self):
+        """Number of parameters in the model."""
+        return self.shape[1]
+
+    def _transform(self, m):
+        return self.P*m
+
+    def deriv(self, m):
+        return self.P
+
 
 class ComplexMap(IdentityMap):
     """ComplexMap
@@ -616,13 +697,13 @@ class CircleMap(IdentityMap):
 
         Parameterize the model space using a circle in a wholespace.
 
-        ..math::
+        .. math::
 
             \sigma(m) = \sigma_1 + (\sigma_2 - \sigma_1)\left(\\arctan\left(100*\sqrt{(\\vec{x}-x_0)^2 + (\\vec{y}-y_0)}-r\\right) \pi^{-1} + 0.5\\right)
 
         Define the model as:
 
-        ..math::
+        .. math::
 
             m = [\sigma_1, \sigma_2, x_0, y_0, r]
 
@@ -975,7 +1056,697 @@ class SplineMap(IdentityMap):
         return sp.csr_matrix(np.c_[g1,g2,g3])
 
 
+class ParametrizedLayer(IdentityMap):
+    """
+        Parametrized Layer Space
+
+        m = [val_background, val_layer, layer_center, layer_thickness]
+
+
+        .. plot::
+            :include-source:
+
+            from SimPEG import Mesh, Maps, np
+            import matplotlib.pyplot as plt
+
+            fig, ax = plt.subplots(1,1,figsize=(2,3))
+
+            mesh = Mesh.TensorMesh([50,50],x0='CC')
+            mapping = Maps.ParametrizedLayer(mesh)
+            m = np.hstack(np.r_[1., 2., -0.1, 0.2])
+            rho = mapping._transform(m)
+            mesh.plotImage(rho, ax=ax)
+
+        **Required**
+
+        :param Mesh mesh: SimPEG Mesh, 2D or 3D
+
+        **Optional**
+
+        :param float slopeFact: arctan slope factor - divided by the minimum h spacing to give the slope of the arctan functions
+        :param float slope: slope of the arctan function
+        :param numpy.ndarray indActive: bool vector with
+
+    """
+
+    slopeFact = 1e2 # will be scaled by the mesh.
+    slope = None
+    indActive = None
+
+    def __init__(self, mesh, **kwargs):
+
+        super(ParametrizedLayer, self).__init__(mesh, **kwargs)
+
+
+        if self.slope is None:
+            self.slope = self.slopeFact / np.hstack(self.mesh.h).min()
+
+        self.x = [self.mesh.gridCC[:,0] if self.indActive is None else self.mesh.gridCC[self.indActive,0]][0]
+
+        if self.mesh.dim > 1:
+            self.y = [self.mesh.gridCC[:,1] if self.indActive is None else self.mesh.gridCC[self.indActive,1]][0]
+
+        if self.mesh.dim > 2:
+            self.z = [self.mesh.gridCC[:,2] if self.indActive is None else self.mesh.gridCC[self.indActive,2]][0]
+
+    @property
+    def nP(self):
+        return 4
+
+    @property
+    def shape(self):
+        if self.indActive is not None:
+            return (sum(self.indActive), self.nP)
+        return (self.mesh.nC, self.nP)
+
+    def mDict(self, m):
+        return {
+                'val_background': m[0],
+                'val_layer': m[1],
+                'layer_center': m[2],
+                'layer_thickness': m[3],
+                }
+
+    def _atanfct(self, xyz, xyzi, slope):
+        return np.arctan(slope * (xyz - xyzi))/np.pi + 0.5
+
+    def _atanfctDeriv(self, xyz, xyzi, slope):
+        # d/dx(atan(x)) = 1/(1+x**2)
+        x = slope * (xyz - xyzi)
+        dx = - slope
+        return (1./(1 + x**2))/np.pi * dx
+
+    def _atanLayer(self, mDict):
+        if self.mesh.dim == 2:
+            z = self.y
+        elif self.mesh.dim == 3:
+            z = self.z
+
+        layer_bottom = mDict['layer_center'] - mDict['layer_thickness'] / 2.
+        layer_top = mDict['layer_center'] + mDict['layer_thickness'] / 2.
+        return self._atanfct(z, layer_bottom, self.slope)*self._atanfct(z, layer_top, -self.slope)
+
+    def _atanLayerDeriv_layer_center(self, mDict):
+        if self.mesh.dim == 2:
+            z = self.y
+        elif self.mesh.dim == 3:
+            z = self.z
+
+        layer_bottom = mDict['layer_center'] - mDict['layer_thickness'] / 2.
+        layer_top = mDict['layer_center'] + mDict['layer_thickness'] / 2.
+
+        return (self._atanfctDeriv(z, layer_bottom, self.slope)*self._atanfct(z, layer_top, -self.slope)
+                + self._atanfct(z, layer_bottom, self.slope)*self._atanfctDeriv(z, layer_top, -self.slope))
+
+    def _atanLayerDeriv_layer_thickness(self, mDict):
+        if self.mesh.dim == 2:
+            z = self.y
+        elif self.mesh.dim == 3:
+            z = self.z
+
+        layer_bottom = mDict['layer_center'] - mDict['layer_thickness'] / 2.
+        layer_top = mDict['layer_center'] + mDict['layer_thickness'] / 2.
+
+        return (-0.5*self._atanfctDeriv(z, layer_bottom, self.slope)*self._atanfct(z, layer_top, -self.slope)
+                + 0.5*self._atanfct(z, layer_bottom, self.slope)*self._atanfctDeriv(z, layer_top, -self.slope))
+
+    def layer_cont(self, mDict):
+        return mDict['val_background'] + (mDict['val_layer'] - mDict['val_background'])*self._atanLayer(mDict)
+
+    def _transform(self, m):
+        mDict = self.mDict(m)
+        return self.layer_cont(mDict)
+
+    def _deriv_val_background(self, mDict):
+        return np.ones_like(self.x) - self._atanLayer(mDict)
+
+    def _deriv_val_layer(self, mDict):
+        return self._atanLayer(mDict)
+
+    def _deriv_layer_center(self, mDict):
+        return (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_center(mDict)
+
+    def _deriv_layer_thickness(self, mDict):
+        return (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_thickness(mDict)
+
+    def deriv(self, m):
+
+        mDict = self.mDict(m)
+
+        return sp.csr_matrix(np.vstack([
+            self._deriv_val_background(mDict),
+            self._deriv_val_layer(mDict),
+            self._deriv_layer_center(mDict),
+            self._deriv_layer_thickness(mDict),
+            ]).T)
+
+
+class ParametrizedCasingAndLayer(ParametrizedLayer):
+    """
+        Parametrized layered space with casing.
+
+        m = [val_background, val_layer, val_casing, val_insideCasing, layer_center, layer_thickness, casing_radius, casing_thickness, casing_bottom, casing_top]
+    """
+
+    def __init__(self, mesh, **kwargs):
+
+        assert mesh._meshType == 'CYL', 'Parametrized Casing in a layer map only works for a cyl mesh.'
+
+        super(ParametrizedCasingAndLayer, self).__init__(mesh, **kwargs)
+
+
+    @property
+    def nP(self):
+        return 10
+
+    @property
+    def shape(self):
+        if self.indActive is not None:
+            return (sum(self.indActive), self.nP)
+        return (self.mesh.nC, self.nP)
+
+    def mDict(self, m):
+        #m = [val_background, val_layer, val_casing, val_insideCasing, layer_center, layer_thickness, casing_radius, casing_thickness, casing_bottom, casing_top]
+        return {
+                    'val_background': m[0],
+                    'val_layer': m[1],
+                    'val_casing': m[2],
+                    'val_insideCasing': m[3],
+                    'layer_center': m[4],
+                    'layer_thickness': m[5],
+                    'casing_radius': m[6],
+                    'casing_thickness': m[7],
+                    'casing_bottom': m[8],
+                    'casing_top': m[9]
+               }
+
+    def _atanCasingLength(self, mDict):
+        return (self._atanfct(self.z, mDict['casing_top'], -self.slope)
+                * self._atanfct(self.z, mDict['casing_bottom'], self.slope))
+
+    def _atanCasingLengthDeriv_casing_top(self, mDict):
+        return (self._atanfctDeriv(self.z, mDict['casing_top'], -self.slope)
+                * self._atanfct(self.z, mDict['casing_bottom'], self.slope))
+
+    def _atanCasingLengthDeriv_casing_bottom(self, mDict):
+        return (self._atanfct(self.z, mDict['casing_top'], -self.slope)
+                * self._atanfctDeriv(self.z, mDict['casing_bottom'], self.slope))
+
+    def _atanInsideCasing(self, mDict):
+        casing_a = mDict['casing_radius'] - 0.5*mDict['casing_thickness']
+        return (self._atanCasingLength(mDict)
+                * self._atanfct(self.x, casing_a, -self.slope))
+
+    def _atanInsideCasingDeriv_casing_radius(self, mDict):
+        casing_a = mDict['casing_radius'] - 0.5*mDict['casing_thickness']
+        return (self._atanCasingLength(mDict)
+                * self._atanfctDeriv(self.x, casing_a, -self.slope))
+
+    def _atanInsideCasingDeriv_casing_thickness(self, mDict):
+        casing_a = mDict['casing_radius'] - 0.5*mDict['casing_thickness']
+        return (self._atanCasingLength(mDict)
+                * - 0.5*self._atanfctDeriv(self.x, casing_a, -self.slope))
+
+    def _atanInsideCasingDeriv_casing_top(self, mDict):
+        casing_a = mDict['casing_radius'] - 0.5*mDict['casing_thickness']
+        return (self._atanCasingLengthDeriv_casing_top(mDict)
+                * self._atanfct(self.x, casing_a, -self.slope))
+
+    def _atanInsideCasingDeriv_casing_bottom(self, mDict):
+        casing_a = mDict['casing_radius'] - 0.5*mDict['casing_thickness']
+        return (self._atanCasingLengthDeriv_casing_bottom(mDict)
+                * self._atanfct(self.x, casing_a, -self.slope))
+
+    def _atanCasing(self, mDict):
+        casing_a, casing_b = mDict['casing_radius'] - 0.5*mDict['casing_thickness'], mDict['casing_radius'] + 0.5*mDict['casing_thickness']
+        return (self._atanCasingLength(mDict)
+                * self._atanfct(self.x, casing_a, self.slope)
+                * self._atanfct(self.x, casing_b, -self.slope))
+
+    def _atanCasingDeriv_casing_radius(self, mDict):
+        casing_a, casing_b = mDict['casing_radius'] - 0.5*mDict['casing_thickness'], mDict['casing_radius'] + 0.5*mDict['casing_thickness']
+        return (self._atanCasingLength(mDict) * (
+                self._atanfctDeriv(self.x, casing_a, self.slope)
+                * self._atanfct(self.x, casing_b, -self.slope)
+                +
+                self._atanfct(self.x, casing_a, self.slope)
+                * self._atanfctDeriv(self.x, casing_b, -self.slope)
+                ))
+
+    def _atanCasingDeriv_casing_thickness(self, mDict):
+        casing_a, casing_b = mDict['casing_radius'] - 0.5*mDict['casing_thickness'], mDict['casing_radius'] + 0.5*mDict['casing_thickness']
+        return (self._atanCasingLength(mDict) * (
+                - 0.5*self._atanfctDeriv(self.x, casing_a, self.slope)
+                * 0.5*self._atanfct(self.x, casing_b, -self.slope)
+                +
+                - 0.5*self._atanfct(self.x, casing_a, self.slope)
+                * 0.5*self._atanfctDeriv(self.x, casing_b, -self.slope)
+                ))
+
+    def _atanCasingDeriv_casing_bottom(self, mDict):
+        casing_a, casing_b = mDict['casing_radius'] - 0.5*mDict['casing_thickness'], mDict['casing_radius'] + 0.5*mDict['casing_thickness']
+        return (self._atanCasingLengthDeriv_casing_bottom(mDict)
+                * self._atanfct(self.x, casing_a, self.slope)
+                * self._atanfct(self.x, casing_b, -self.slope))
+
+    def _atanCasingDeriv_casing_top(self, mDict):
+        casing_a, casing_b = mDict['casing_radius'] - 0.5*mDict['casing_thickness'], mDict['casing_radius'] + 0.5*mDict['casing_thickness']
+        return (self._atanCasingLengthDeriv_casing_top(mDict)
+                * self._atanfct(self.x, casing_a, self.slope)
+                * self._atanfct(self.x, casing_b, -self.slope))
+
+    def layer_cont(self, mDict):
+        return mDict['val_background'] + (mDict['val_layer']-mDict['val_background']) * self._atanLayer(mDict) # contribution from the layered background
+
+
+    def _transform(self, m):
+
+        mDict = self.mDict(m)
+
+        # assemble the model
+        layer = self.layer_cont(mDict)
+        casing = (mDict['val_casing'] - layer) * self._atanCasing(mDict)
+        insideCasing = (mDict['val_insideCasing'] - layer) * self._atanInsideCasing(mDict)
+
+        return layer + casing + insideCasing
+
+
+    def _deriv_val_background(self, mDict):
+        d_layer_cont_dval_background = 1. - self._atanLayer(mDict) # contribution from the layered background
+        d_casing_cont_dval_background = -1. * d_layer_cont_dval_background * self._atanCasing(mDict)
+        d_insideCasing_cont_dval_background = -1. * d_layer_cont_dval_background * self._atanInsideCasing(mDict)
+        return d_layer_cont_dval_background + d_casing_cont_dval_background + d_insideCasing_cont_dval_background
+
+    def _deriv_val_layer(self, mDict):
+        d_layer_cont_dval_layer = self._atanLayer(mDict)
+        d_casing_cont_dval_layer = -1. * d_layer_cont_dval_layer * self._atanCasing(mDict)
+        d_insideCasing_cont_dval_layer = -1. * d_layer_cont_dval_layer * self._atanInsideCasing(mDict)
+        return d_layer_cont_dval_layer + d_casing_cont_dval_layer + d_insideCasing_cont_dval_layer
+
+    def _deriv_val_casing(self, mDict):
+        d_layer_cont_dval_casing = 0.
+        d_casing_cont_dval_casing = self._atanCasing(mDict)
+        d_insideCasing_cont_dval_casing = 0.
+        return d_layer_cont_dval_casing + d_casing_cont_dval_casing + d_insideCasing_cont_dval_casing
+
+    def _deriv_val_insideCasing(self, mDict):
+        d_layer_cont_dval_insideCasing = 0.
+        d_casing_cont_dval_insideCasing = 0.
+        d_insideCasing_cont_dval_insideCasing = self._atanInsideCasing(mDict)
+        return d_layer_cont_dval_insideCasing + d_casing_cont_dval_insideCasing + d_insideCasing_cont_dval_insideCasing
+
+    def _deriv_layer_center(self, mDict):
+        d_layer_cont_dlayer_center = (mDict['val_layer'] - mDict['val_background']) * self._atanLayerDeriv_layer_center(mDict)
+        d_casing_cont_dlayer_center =  - d_layer_cont_dlayer_center * self._atanCasing(mDict)
+        d_insideCasing_cont_dlayer_center = - d_layer_cont_dlayer_center * self._atanInsideCasing(mDict)
+        return d_layer_cont_dlayer_center + d_casing_cont_dlayer_center + d_insideCasing_cont_dlayer_center
+
+    def _deriv_layer_thickness(self, mDict):
+        d_layer_cont_dlayer_thickness = (mDict['val_layer']-mDict['val_background']) * self._atanLayerDeriv_layer_thickness(mDict)
+        d_casing_cont_dlayer_thickness =  - d_layer_cont_dlayer_thickness * self._atanCasing(mDict)
+        d_insideCasing_cont_dlayer_thickness = - d_layer_cont_dlayer_thickness * self._atanInsideCasing(mDict)
+        return d_layer_cont_dlayer_thickness + d_casing_cont_dlayer_thickness + d_insideCasing_cont_dlayer_thickness
+
+    def _deriv_casing_radius(self, mDict):
+        layer = self.layer_cont(mDict)
+        d_layer_cont_dcasing_radius = 0.
+        d_casing_cont_dcasing_radius = (mDict['val_casing'] - layer) * self._atanCasingDeriv_casing_radius(mDict)
+        d_insideCasing_cont_dcasing_radius = (mDict['val_insideCasing'] - layer) * self._atanInsideCasingDeriv_casing_radius(mDict)
+        return d_layer_cont_dcasing_radius + d_casing_cont_dcasing_radius + d_insideCasing_cont_dcasing_radius
+
+    def _deriv_casing_thickness(self, mDict):
+        d_layer_cont_dcasing_thickness = 0.
+        d_casing_cont_dcasing_thickness = (mDict['val_casing'] - self.layer_cont(mDict)) * self._atanCasingDeriv_casing_thickness(mDict)
+        d_insideCasing_cont_dcasing_thickness = (mDict['val_insideCasing'] - self.layer_cont(mDict)) * self._atanInsideCasingDeriv_casing_thickness(mDict)
+        return d_layer_cont_dcasing_thickness + d_casing_cont_dcasing_thickness + d_insideCasing_cont_dcasing_thickness
+
+    def _deriv_casing_bottom(self, mDict):
+        d_layer_cont_dcasing_bottom = 0.
+        d_casing_cont_dcasing_bottom = (mDict['val_casing'] - self.layer_cont(mDict)) * self._atanCasingDeriv_casing_bottom(mDict)
+        d_insideCasing_cont_dcasing_bottom = (mDict['val_insideCasing'] - self.layer_cont(mDict)) * self._atanInsideCasingDeriv_casing_bottom(mDict)
+        return d_layer_cont_dcasing_bottom + d_casing_cont_dcasing_bottom + d_insideCasing_cont_dcasing_bottom
+
+    def _deriv_casing_top(self, mDict):
+        d_layer_cont_dcasing_top = 0.
+        d_casing_cont_dcasing_top = (mDict['val_casing'] - self.layer_cont(mDict)) * self._atanCasingDeriv_casing_top(mDict)
+        d_insideCasing_cont_dcasing_top = (mDict['val_insideCasing'] - self.layer_cont(mDict)) * self._atanInsideCasingDeriv_casing_top(mDict)
+        return d_layer_cont_dcasing_top + d_casing_cont_dcasing_top + d_insideCasing_cont_dcasing_top
+
+
+    def deriv(self, m):
+
+        mDict = self.mDict(m)
+
+        return sp.csr_matrix(np.vstack([
+                self._deriv_val_background(mDict),
+                self._deriv_val_layer(mDict),
+                self._deriv_val_casing(mDict),
+                self._deriv_val_insideCasing(mDict),
+                self._deriv_layer_center(mDict),
+                self._deriv_layer_thickness(mDict),
+                self._deriv_casing_radius(mDict),
+                self._deriv_casing_thickness(mDict),
+                self._deriv_casing_bottom(mDict),
+                self._deriv_casing_top(mDict),
+                ]).T)
 
 
 
+class ParametrizedBlockInLayer(ParametrizedLayer):
+    """
+        Parametrized Block in a Layered Space
+
+        For 2D:
+        m = [val_background, val_layer, val_block, layer_center, layer_thickness, block_x0, block_dx]
+
+        For 3D:
+        m = [val_background, val_layer, val_block, layer_center, layer_thickness, block_x0, block_y0, block_dx, block_dy]
+
+        .. plot::
+            :include-source:
+
+            from SimPEG import Mesh, Maps, np
+            import matplotlib.pyplot as plt
+
+            fig, ax = plt.subplots(1,1,figsize=(2,3))
+
+            mesh = Mesh.TensorMesh([50,50],x0='CC')
+            mapping = Maps.ParametrizedBlockInLayer(mesh)
+            m = np.hstack(np.r_[1., 2., 3., -0.1, 0.2, 0.3, 0.2])
+            rho = mapping._transform(m)
+            mesh.plotImage(rho, ax=ax)
+
+        **Required**
+
+        :param Mesh mesh: SimPEG Mesh, 2D or 3D
+
+        **Optional**
+
+        :param float slopeFact: arctan slope factor - divided by the minimum h spacing to give the slope of the arctan functions
+        :param float slope: slope of the arctan function
+        :param numpy.ndarray indActive: bool vector with
+
+    """
+
+    def __init__(self, mesh, **kwargs):
+
+        super(ParametrizedBlockInLayer, self).__init__(mesh, **kwargs)
+
+    @property
+    def nP(self):
+        if self.mesh.dim == 2:
+            return 7
+        elif self.mesh.dim == 3:
+            return 9
+
+    @property
+    def shape(self):
+        if self.indActive is not None:
+            return (sum(self.indActive), self.nP)
+        return (self.mesh.nC, self.nP)
+
+    def _mDict2d(self, m):
+        return{
+                'val_background': m[0],
+                'val_layer': m[1],
+                'val_block': m[2],
+                'layer_center': m[3],
+                'layer_thickness': m[4],
+                'x0_block': m[5],
+                'dx_block': m[6]
+              }
+
+    def _mDict3d(self, m):
+        return{
+                'val_background': m[0],
+                'val_layer': m[1],
+                'val_block': m[2],
+                'layer_center': m[3],
+                'layer_thickness': m[4],
+                'x0_block': m[5],
+                'y0_block': m[6],
+                'dx_block': m[7],
+                'dy_block': m[8]
+              }
+
+    def mDict(self, m):
+        if self.mesh.dim == 2:
+            return self._mDict2d(m)
+        elif self.mesh.dim == 3:
+            return self._mDict3d(m)
+
+    def _atanBlock2d(self, mDict):
+        return (self._atanLayer(mDict)
+                          * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                          * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
+
+    def _atanBlock2dDeriv_layer_center(self, mDict):
+        return (self._atanLayerDeriv_layer_center(mDict)
+                          * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                          * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
+
+    def _atanBlock2dDeriv_layer_thickness(self, mDict):
+        return (self._atanLayerDeriv_layer_thickness(mDict)
+                          * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                          * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
+
+
+    def _atanBlock2dDeriv_x0(self, mDict):
+        return self._atanLayer(mDict) * (
+                    (self._atanfctDeriv(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
+                    +
+                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                    * self._atanfctDeriv(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
+                    )
+
+    def _atanBlock2dDeriv_dx(self, mDict):
+        return self._atanLayer(mDict) * (
+                    (self._atanfctDeriv(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope) * -0.5
+                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
+                    +
+                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                    * self._atanfctDeriv(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope) * 0.5)
+                    )
+
+    def _atanBlock3d(self, mDict):
+        return (self._atanLayer(mDict)
+                      * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                      * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
+                      * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
+                      * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
+
+
+    def _atanBlock3dDeriv_layer_center(self, mDict):
+        return (self._atanLayerDeriv_layer_center(mDict)
+                      * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                      * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
+                      * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
+                      * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
+
+    def _atanBlock3dDeriv_layer_thickness(self, mDict):
+        return (self._atanLayerDeriv_layer_thickness(mDict)
+                      * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                      * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
+                      * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
+                      * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
+
+
+    def _atanBlock3dDeriv_x0(self, mDict):
+        return self._atanLayer(mDict) * (
+                    (self._atanfctDeriv(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
+                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
+                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
+                    +
+                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                    * self._atanfctDeriv(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
+                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
+                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
+                    )
+
+    def _atanBlock3dDeriv_y0(self, mDict):
+        return self._atanLayer(mDict) * (
+                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
+                    * self._atanfctDeriv(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
+                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
+                    +
+                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
+                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
+                    * self._atanfctDeriv(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
+                    )
+
+    def _atanBlock3dDeriv_dx(self, mDict):
+        return self._atanLayer(mDict) * (
+                    (self._atanfctDeriv(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope) * -0.5
+                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
+                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
+                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
+                    +
+                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                    * self._atanfctDeriv(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope) * 0.5
+                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
+                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
+                    )
+
+    def _atanBlock3dDeriv_dy(self, mDict):
+        return self._atanLayer(mDict) * (
+                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
+                    * self._atanfctDeriv(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope) * -0.5
+                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
+                    +
+                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
+                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
+                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
+                    * self._atanfctDeriv(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope) * 0.5)
+                    )
+
+
+    def _transform2d(self, m):
+        mDict = self.mDict(m)
+        # assemble the model
+        layer_cont = mDict['val_background'] + (mDict['val_layer']-mDict['val_background'])*self._atanLayer(mDict) # contribution from the layered background
+        block_cont = (mDict['val_block']-layer_cont)*self._atanBlock2d(mDict) # perturbation due to the block
+
+        return layer_cont + block_cont
+
+    def _deriv2d_val_background(self, mDict):
+        d_layer_dval_background = np.ones_like(self.x) - self._atanLayer(mDict)
+        d_block_dval_background = (-d_layer_dval_background)*self._atanBlock2d(mDict)
+        return d_layer_dval_background + d_block_dval_background
+
+    def _deriv2d_val_layer(self, mDict):
+        d_layer_dval_layer = self._atanLayer(mDict)
+        d_block_dval_layer = (-d_layer_dval_layer)*self._atanBlock2d(mDict)
+        return d_layer_dval_layer + d_block_dval_layer
+
+    def _deriv2d_val_block(self, mDict):
+        d_layer_dval_block = 0.
+        d_block_dval_block = (1.-d_layer_dval_block)*self._atanBlock2d(mDict)
+        return d_layer_dval_block + d_block_dval_block
+
+    def _deriv2d_layer_center(self, mDict):
+        d_layer_dlayer_center = (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_center(mDict)
+        d_block_dlayer_center = ((mDict['val_block']-self.layer_cont(mDict))*self._atanBlock2dDeriv_layer_center(mDict)
+                                    - d_layer_dlayer_center*self._atanBlock2d(mDict))
+        return d_layer_dlayer_center + d_block_dlayer_center
+
+    def _deriv2d_layer_thickness(self, mDict):
+        d_layer_dlayer_thickness = (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_thickness(mDict)
+        d_block_dlayer_thickness = ((mDict['val_block']-self.layer_cont(mDict))*self._atanBlock2dDeriv_layer_thickness(mDict)
+                                    - d_layer_dlayer_thickness*self._atanBlock2d(mDict))
+        return d_layer_dlayer_thickness + d_block_dlayer_thickness
+
+    def _deriv2d_x0_block(self, mDict):
+        d_layer_dx0 = 0.
+        d_block_dx0 = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock2dDeriv_x0(mDict)
+        return d_layer_dx0 + d_block_dx0
+
+    def _deriv2d_dx_block(self, mDict):
+        d_layer_ddx = 0.
+        d_block_ddx = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock2dDeriv_dx(mDict)
+        return d_layer_ddx + d_block_ddx
+
+    def _deriv2d(self, m):
+        mDict = self.mDict(m)
+
+        return np.vstack([
+                            self._deriv2d_val_background(mDict),
+                            self._deriv2d_val_layer(mDict),
+                            self._deriv2d_val_block(mDict),
+                            self._deriv2d_layer_center(mDict),
+                            self._deriv2d_layer_thickness(mDict),
+                            self._deriv2d_x0_block(mDict),
+                            self._deriv2d_dx_block(mDict)
+                        ]).T
+
+    def _transform3d(self, m):
+        # parse model
+        mDict = self.mDict(m)
+
+        # assemble the model
+        layer_cont = mDict['val_background'] + (mDict['val_layer']-mDict['val_background'])*self._atanLayer(mDict) # contribution from the layered background
+        block_cont = (mDict['val_block']-layer_cont)*self._atanBlock3d(mDict) # perturbation due to the block
+
+        return layer_cont + block_cont
+
+    def _deriv3d_val_background(self, mDict):
+        d_layer_dval_background = np.ones_like(self.x) - self._atanLayer(mDict)
+        d_block_dval_background = (-d_layer_dval_background)*self._atanBlock3d(mDict)
+        return d_layer_dval_background + d_block_dval_background
+
+    def _deriv3d_val_layer(self, mDict):
+        d_layer_dval_layer = self._atanLayer(mDict)
+        d_block_dval_layer = (-d_layer_dval_layer)*self._atanBlock3d(mDict)
+        return d_layer_dval_layer + d_block_dval_layer
+
+    def _deriv3d_val_block(self, mDict):
+        d_layer_dval_block = 0.
+        d_block_dval_block = (1.-d_layer_dval_block)*self._atanBlock3d(mDict)
+        return d_layer_dval_block + d_block_dval_block
+
+    def _deriv3d_layer_center(self, mDict):
+        d_layer_dlayer_center = (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_center(mDict)
+        d_block_dlayer_center = ((mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_layer_center(mDict)
+                                    - d_layer_dlayer_center*self._atanBlock3d(mDict))
+        return d_layer_dlayer_center + d_block_dlayer_center
+
+    def _deriv3d_layer_thickness(self, mDict):
+        d_layer_dlayer_thickness = (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_thickness(mDict)
+        d_block_dlayer_thickness = ((mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_layer_thickness(mDict)
+                                - d_layer_dlayer_thickness*self._atanBlock3d(mDict))
+        return d_layer_dlayer_thickness + d_block_dlayer_thickness
+
+    def _deriv3d_x0_block(self, mDict):
+        d_layer_dx0 = 0.
+        d_block_dx0 = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_x0(mDict)
+        return d_layer_dx0 + d_block_dx0
+
+    def _deriv3d_y0_block(self, mDict):
+        d_layer_dy0 = 0.
+        d_block_dy0 = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_y0(mDict)
+        return d_layer_dy0 + d_block_dy0
+
+    def _deriv3d_dx_block(self, mDict):
+        d_layer_ddx = 0.
+        d_block_ddx = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_dx(mDict)
+        return d_layer_ddx + d_block_ddx
+
+    def _deriv3d_dy_block(self, mDict):
+        d_layer_ddy = 0.
+        d_block_ddy = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_dy(mDict)
+        return d_layer_ddy + d_block_ddy
+
+    def _deriv3d(self, m):
+
+        mDict = self.mDict(m)
+
+        return np.vstack([
+                            self._deriv3d_val_background(mDict),
+                            self._deriv3d_val_layer(mDict),
+                            self._deriv3d_val_block(mDict),
+                            self._deriv3d_layer_center(mDict),
+                            self._deriv3d_layer_thickness(mDict),
+                            self._deriv3d_x0_block(mDict),
+                            self._deriv3d_y0_block(mDict),
+                            self._deriv3d_dx_block(mDict),
+                            self._deriv3d_dy_block(mDict),
+                        ]).T
+
+    def _transform(self, m):
+
+        if self.mesh.dim == 2:
+            return self._transform2d(m)
+        elif self.mesh.dim == 3:
+            return self._transform3d(m)
+
+    def deriv(self, m):
+
+        if self.mesh.dim == 2:
+            return sp.csr_matrix(self._deriv2d(m))
+        elif self.mesh.dim == 3:
+            return sp.csr_matrix(self._deriv3d(m))
 

SimPEG/Mesh/CurvilinearMesh.py

@@ -2,7 +2,6 @@ from SimPEG import Utils, np
 from BaseMesh import BaseRectangularMesh
 from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
-from View import CurvView
 
 # Some helper functions.
 length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
@@ -11,7 +10,7 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
 normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
 
 
-class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView):
+class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
     """
     CurvilinearMesh is a mesh class that deals with curvilinear meshes.
 
@@ -331,6 +330,102 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvVie
 
 
 
+    #############################################
+    #            Plotting Functions             #
+    #############################################
+
+    def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True,  showIt=False):
+        """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
+
+
+        .. plot::
+            :include-source:
+
+            from SimPEG import Mesh, Utils
+            X, Y = Utils.exampleLrmGrid([3,3],'rotate')
+            M = Mesh.CurvilinearMesh([X, Y])
+            M.plotGrid(showIt=True)
+
+        """
+        import matplotlib.pyplot as plt
+        import matplotlib
+        from mpl_toolkits.mplot3d import Axes3D
+        mkvc = Utils.mkvc
+
+        axOpts = {'projection':'3d'} if self.dim == 3 else {}
+        if ax is None: ax = plt.subplot(111, **axOpts)
+
+        NN = self.r(self.gridN, 'N', 'N', 'M')
+        if self.dim == 2:
+
+            if lines:
+                X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
+                Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
+
+                X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
+                Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
+
+                X = np.r_[X1, X2]
+                Y = np.r_[Y1, Y2]
+
+                ax.plot(X, Y, 'b-')
+            if centers:
+                ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
+
+            # Nx = self.r(self.normals, 'F', 'Fx', 'V')
+            # Ny = self.r(self.normals, 'F', 'Fy', 'V')
+            # Tx = self.r(self.tangents, 'E', 'Ex', 'V')
+            # Ty = self.r(self.tangents, 'E', 'Ey', 'V')
+
+            # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
+
+            # nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
+            # ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
+            # ax.plot(nX, nY, 'r-')
+
+            # nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
+            # #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
+            # ax.plot(nX, nY, 'g-')
+
+            # tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
+            # tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
+            # ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
+            # ax.plot(tX, tY, 'r-')
+
+            # nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
+            # #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
+            # ax.plot(nX, nY, 'g-')
+
+        elif self.dim == 3:
+            X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
+            Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
+            Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
+
+            X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten()
+            Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten()
+            Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten()
+
+            X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten()
+            Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten()
+            Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten()
+
+            X = np.r_[X1, X2, X3]
+            Y = np.r_[Y1, Y2, Y3]
+            Z = np.r_[Z1, Z2, Z3]
+
+            ax.plot(X, Y, 'b', zs=Z)
+            ax.set_zlabel('x3')
+
+        ax.grid(True)
+        ax.set_xlabel('x1')
+        ax.set_ylabel('x2')
+
+        if showIt: plt.show()
+
+
 if __name__ == '__main__':
     nc = 5
     h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])

SimPEG/Mesh/View.py

@@ -552,8 +552,7 @@ class CurvView(object):
     def __init__(self):
         pass
 
-
-    def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True,  showIt=False):
+    def plotGrid(self, length=0.05, showIt=False):
         """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
 
 
@@ -561,63 +560,60 @@ class CurvView(object):
             :include-source:
 
             from SimPEG import Mesh, Utils
-            X, Y = Utils.exampleLrmGrid([3,3],'rotate')
+            X, Y = Utils.exampleCurvGird([3,3],'rotate')
             M = Mesh.CurvilinearMesh([X, Y])
             M.plotGrid(showIt=True)
 
         """
-        import matplotlib.pyplot as plt
-        import matplotlib
-        from mpl_toolkits.mplot3d import Axes3D
-
-        axOpts = {'projection':'3d'} if self.dim == 3 else {}
-        if ax is None: ax = plt.subplot(111, **axOpts)
-
         NN = self.r(self.gridN, 'N', 'N', 'M')
         if self.dim == 2:
+            fig = plt.figure(2)
+            fig.clf()
+            ax = plt.subplot(111)
+            X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
+            Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
 
-            if lines:
-                X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
-                Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
+            X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
+            Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
 
-                X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
-                Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
+            X = np.r_[X1, X2]
+            Y = np.r_[Y1, Y2]
 
-                X = np.r_[X1, X2]
-                Y = np.r_[Y1, Y2]
+            plt.plot(X, Y)
 
-                ax.plot(X, Y, 'b-')
-            if centers:
-                ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
+            plt.hold(True)
+            Nx = self.r(self.normals, 'F', 'Fx', 'V')
+            Ny = self.r(self.normals, 'F', 'Fy', 'V')
+            Tx = self.r(self.tangents, 'E', 'Ex', 'V')
+            Ty = self.r(self.tangents, 'E', 'Ey', 'V')
 
-            # Nx = self.r(self.normals, 'F', 'Fx', 'V')
-            # Ny = self.r(self.normals, 'F', 'Fy', 'V')
-            # Tx = self.r(self.tangents, 'E', 'Ex', 'V')
-            # Ty = self.r(self.tangents, 'E', 'Ey', 'V')
+            plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
 
-            # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
+            nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
+            nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
+            plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
+            plt.plot(nX, nY, 'r-')
 
-            # nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
-            # ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
-            # ax.plot(nX, nY, 'r-')
+            nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
+            nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
+            #plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
+            plt.plot(nX, nY, 'g-')
 
-            # nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
-            # #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
-            # ax.plot(nX, nY, 'g-')
+            tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
+            tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
+            plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
+            plt.plot(tX, tY, 'r-')
 
-            # tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
-            # tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
-            # ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
-            # ax.plot(tX, tY, 'r-')
-
-            # nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
-            # #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
-            # ax.plot(nX, nY, 'g-')
+            nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
+            nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
+            #plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
+            plt.plot(nX, nY, 'g-')
+            plt.axis('equal')
 
         elif self.dim == 3:
+            fig = plt.figure(3)
+            fig.clf()
+            ax = fig.add_subplot(111, projection='3d')
             X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
             Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
             Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
@@ -634,50 +630,16 @@ class CurvView(object):
             Y = np.r_[Y1, Y2, Y3]
             Z = np.r_[Z1, Z2, Z3]
 
-            ax.plot(X, Y, 'b', zs=Z)
+            plt.plot(X, Y, 'b', zs=Z)
             ax.set_zlabel('x3')
 
         ax.grid(True)
+        ax.hold(False)
         ax.set_xlabel('x1')
         ax.set_ylabel('x2')
 
         if showIt: plt.show()
 
-    def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
-        if self.dim == 3: raise NotImplementedError('This is not yet done!')
-
-        import matplotlib.pyplot as plt
-        import matplotlib
-        from mpl_toolkits.mplot3d import Axes3D
-        import matplotlib.colors as colors
-        import matplotlib.cm as cmx
-
-        if ax is None: ax = plt.subplot(111)
-        jet = cm = plt.get_cmap('jet')
-        cNorm  = colors.Normalize(
-            vmin=I.min() if clim is None else clim[0],
-            vmax=I.max() if clim is None else clim[1])
-
-        scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
-        # ax.set_xlim((self.x0[0], self.h[0].sum()))
-        # ax.set_ylim((self.x0[1], self.h[1].sum()))
-
-        Nx = self.r(self.gridN[:,0],'N','N','M')
-        Ny = self.r(self.gridN[:,1],'N','N','M')
-        cell = self.r(I,'CC','CC','M')
-
-        for ii in range(self.nCx):
-            for jj in range(self.nCy):
-                I = [ii,ii+1,ii+1,ii]
-                J = [jj,jj,jj+1,jj+1]
-                ax.add_patch(plt.Polygon(np.c_[Nx[I,J],Ny[I,J]], facecolor=scalarMap.to_rgba(cell[ii,jj]), edgecolor='k' if grid else 'none'))
-
-        scalarMap._A = []  # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
-        ax.set_xlabel('x')
-        ax.set_ylabel('y')
-        if showIt: plt.show()
-        return [scalarMap]
-
 
 if __name__ == '__main__':
     from SimPEG import *

SimPEG/Optimization.py

@@ -1008,4 +1008,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
             indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
             delx[indx] = 0.
 
-            return delx
+            return delx

SimPEG/Regularization.py

@@ -1,6 +1,4 @@
-import Utils, Maps, Mesh
-import numpy as np
-import scipy.sparse as sp
+import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
 
 class RegularizationMesh(object):
     """
@@ -41,7 +39,7 @@ class RegularizationMesh(object):
             if self.indActive is None:
                 self._nC = self.mesh.nC
             else:
-                self._nC = sum(self.indActive)
+                self._nC = int(sum(self.indActive))
         return self._nC
 
     @property
@@ -306,7 +304,7 @@ class BaseRegularization(object):
     mesh    = None    #: A SimPEG.Mesh instance.
     mref    = None    #: Reference model.
 
-    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+    def __init__(self, mesh=None, nP=None, mapping=None, indActive=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
         assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
         if indActive is not None and indActive.dtype != 'bool':
@@ -316,11 +314,19 @@ class BaseRegularization(object):
         if indActive is not None and mapping is None:
             mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
 
+        if mesh is None and nP is None:
+            raise Exception, 'either Mesh or number of parameters must be provided to the BaseRegularization'
+
         self.regmesh = RegularizationMesh(mesh,indActive)
-        self.mapping = mapping or self.mapPair(mesh)
-        self.mapping._assertMatchesPair(self.mapPair)
         self.indActive = indActive
 
+        if mesh is not None and nP is None:
+            nP = self.regmesh.nC
+        self.nP = nP
+
+        self.mapping = mapping or self.mapPair(nP=self.nP)
+        self.mapping._assertMatchesPair(self.mapPair)
+
     @property
     def parent(self):
         """This is the parent of the regularization."""
@@ -348,7 +354,7 @@ class BaseRegularization(object):
     @property
     def W(self):
         """Full regularization weighting matrix W."""
-        return sp.identity(self.regmesh.nC)
+        return sp.identity(self.nP)
 
     @Utils.timeIt
     def eval(self, m):
@@ -405,238 +411,7 @@ class BaseRegularization(object):
 
         return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
 
-class Simple(BaseRegularization):
-    """
-    Simple regularization that does not include length scales in the derivatives.
-    """
-
-    mrefInSmooth = False  #: include mref in the smoothness?
-    alpha_s      = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
-    alpha_x      = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'],     "Weight for the first derivative in the x direction")
-    alpha_y      = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'],     "Weight for the first derivative in the y direction")
-    alpha_z      = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'],     "Weight for the first derivative in the z direction")
-    cell_weights = 1.
-
-    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
-
-        if isinstance(self.cell_weights,float):
-            self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
-
-    @property
-    def Wsmall(self):
-        """Regularization matrix Wsmall"""
-        if getattr(self,'_Wsmall', None) is None:
-            self._Wsmall = Utils.sdiag((self.alpha_s*self.cell_weights)**0.5)
-        return self._Wsmall
-
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
-        if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.alpha_x * (self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.regmesh.cellDiffxStencil
-        return self._Wx
-
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
-        if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.alpha_y * (self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.regmesh.cellDiffyStencil
-        return self._Wy
-
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
-        if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.alpha_z * (self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.regmesh.cellDiffzStencil
-        return self._Wz
-
-#    @property
-#    def Wsmooth(self):
-#        """Full smoothness regularization matrix W"""
-#        print 'wtf why are we using Wsmooth'
-#        raise NotImplementedError
-#        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"""
-#        print 'wtf why are we using W'
-#        if getattr(self, '_W', None) is None:
-#            wlist = (self.Wsmall, self.Wx)
-#            if self.regmesh.dim > 1:
-#                wlist += (self.Wy,)
-#            if self.regmesh.dim > 2:
-#                wlist += (self.Wz,)
-#            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 _evalSmallDeriv(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
-
-    @Utils.timeIt
-    def _evalSmall2Deriv(self, m, v = None):
-        rDeriv = self.Wsmall * ( self.mapping.deriv(m - self.mref) )
-        if v is not None:
-            return rDeriv.T * (rDeriv * v)
-        return rDeriv.T * rDeriv
-
-    @Utils.timeIt
-    def _evalSmoothx(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wx * ( self.mapping * (m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wx * ( self.mapping * (m) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmoothy(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wy * ( self.mapping * (m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wy * ( self.mapping * (m) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmoothz(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wz * ( self.mapping * (m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wz * ( self.mapping * (m) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
-        phiSmooth = self._evalSmoothx(m)
-        if self.regmesh.dim > 1:
-            phiSmooth += self._evalSmoothy(m)
-        if self.regmesh.dim > 2:
-            phiSmooth += self._evalSmoothz(m)
-        return phiSmooth
-
-    @Utils.timeIt
-    def _evalSmoothxDeriv(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wx * ( self.mapping * ( m - self.mref ) )
-            return r.T * ( self.Wx * self.mapping.deriv(m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wx * ( self.mapping * m )
-            return r.T * ( self.Wx * self.mapping.deriv(m) )
-
-    @Utils.timeIt
-    def _evalSmoothx2Deriv(self, m, v=None):
-        if self.mrefInSmooth == True:
-            rDeriv = self.Wx * ( self.mapping.deriv( m - self.mref ) )
-        elif self.mrefInSmooth == False:
-            rDeriv = self.Wx * ( self.mapping.deriv(m) )
-
-        if v is not None:
-            return rDeriv.T * ( rDeriv * v )
-        return rDeriv.T * rDeriv
-
-    @Utils.timeIt
-    def _evalSmoothyDeriv(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wy * ( self.mapping * ( m - self.mref ) )
-            return r.T * ( self.Wy * self.mapping.deriv(m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wy * ( self.mapping * m )
-            return r.T * ( self.Wy * self.mapping.deriv(m) )
-
-    @Utils.timeIt
-    def _evalSmoothy2Deriv(self, m, v=None):
-        if self.mrefInSmooth == True:
-            rDeriv = self.Wy * ( self.mapping.deriv( m - self.mref ) )
-        elif self.mrefInSmooth == False:
-            rDeriv = self.Wy * ( self.mapping.deriv(m) )
-
-        if v is not None:
-            return rDeriv.T * ( rDeriv * v )
-        return rDeriv.T * rDeriv
-
-    @Utils.timeIt
-    def _evalSmoothzDeriv(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wz * ( self.mapping * ( m - self.mref ) )
-            return r.T * ( self.Wz * self.mapping.deriv(m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wz * ( self.mapping * m )
-            return r.T * ( self.Wz * self.mapping.deriv(m) )
-
-    @Utils.timeIt
-    def _evalSmoothz2Deriv(self, m, v=None):
-        if self.mrefInSmooth == True:
-            rDeriv = self.Wz * ( self.mapping.deriv( m - self.mref ) )
-        elif self.mrefInSmooth == False:
-            rDeriv = self.Wz * ( self.mapping.deriv(m) )
-
-        if v is not None:
-            return rDeriv.T * ( rDeriv * v )
-        return rDeriv.T * rDeriv
-
-    @Utils.timeIt
-    def _evalSmoothDeriv(self, m):
-        deriv = self._evalSmoothxDeriv(m)
-        if self.regmesh.dim > 1:
-            deriv += self._evalSmoothyDeriv(m)
-        if self.regmesh.dim > 2:
-            deriv += self._evalSmoothzDeriv(m)
-        return deriv
-
-    @Utils.timeIt
-    def _evalSmooth2Deriv(self, m, v=None):
-        deriv = self._evalSmoothx2Deriv(m, v)
-        if self.regmesh.dim > 1:
-            deriv += self._evalSmoothy2Deriv(m, v)
-        if self.regmesh.dim > 2:
-            deriv += self._evalSmoothz2Deriv(m, v)
-        return deriv
-
-
-    @Utils.timeIt
-    def eval(self, m):
-        return self._evalSmall(m) + self._evalSmooth(m)
-
-    @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})}
-
-        """
-        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
-
-    @Utils.timeIt
-    def eval2Deriv(self, m, v=None):
-        return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v)
-
-
-
-class Tikhonov(Simple):
+class Tikhonov(BaseRegularization):
     """
     L2 Tikhonov regularization with both smallness and smoothness (first order
     derivative) contributions.
@@ -726,131 +501,56 @@ class Tikhonov(Simple):
             self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
         return self._Wzz
 
-
     @property
-    def Wsmooth2(self):
+    def Wsmooth(self):
         """Full smoothness regularization matrix W"""
         if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wxx)
+            wlist = (self.Wx, self.Wxx)
             if self.regmesh.dim > 1:
-                wlist += (self.Wyy)
+                wlist += (self.Wy, self.Wyy)
             if self.regmesh.dim > 2:
-                wlist += (self.Wzz)
+                wlist += (self.Wz, self.Wzz)
             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 _evalSmoothxx(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wxx * ( self.mapping * (m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wxx * ( self.mapping * (m) )
+    def _evalSmall(self, m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
         return 0.5 * r.dot(r)
 
     @Utils.timeIt
-    def _evalSmoothyy(self, m):
+    def _evalSmooth(self, m):
         if self.mrefInSmooth == True:
-            r = self.Wyy * ( self.mapping * (m - self.mref) )
+            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
         elif self.mrefInSmooth == False:
-            r = self.Wyy * ( self.mapping * (m) )
+            r = self.Wsmooth * ( self.mapping * (m) )
         return 0.5 * r.dot(r)
 
-    @Utils.timeIt
-    def _evalSmoothzz(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wzz * ( self.mapping * (m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wzz * ( self.mapping * (m) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth2(self, m):
-        phiSmooth2 = self._evalSmoothxx(m)
-        if self.regmesh.dim > 1:
-            phiSmooth2 += self._evalSmoothyy(m)
-        if self.regmesh.dim > 2:
-            phiSmooth2 += self._evalSmoothzz(m)
-        return phiSmooth2
-
-    @Utils.timeIt
-    def _evalSmoothxxDeriv(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wxx * ( self.mapping * ( m - self.mref ) )
-            return r.T * ( self.Wxx * self.mapping.deriv(m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wxx * ( self.mapping * m )
-            return r.T * ( self.Wxx * self.mapping.deriv(m) )
-
-    @Utils.timeIt
-    def _evalSmoothyyDeriv(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wyy * ( self.mapping * ( m - self.mref ) )
-            return r.T * ( self.Wyy * self.mapping.deriv(m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wyy * ( self.mapping * m )
-            return r.T * ( self.Wyy * self.mapping.deriv(m) )
-
-    @Utils.timeIt
-    def _evalSmoothzzDeriv(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wzz * ( self.mapping * ( m - self.mref ) )
-            return r.T * ( self.Wzz * self.mapping.deriv(m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wzz * ( self.mapping * m )
-            return r.T * ( self.Wzz * self.mapping.deriv(m) )
-
-    @Utils.timeIt
-    def _evalSmoothxx2Deriv(self, m, v=None):
-        if self.mrefInSmooth == True:
-            rDeriv = self.Wxx * ( self.mapping.deriv( m - self.mref ) )
-        elif self.mrefInSmooth == False:
-            rDeriv = self.Wxx * self.mapping.deriv(m)
-        if v is not None:
-            return rDeriv.T * (rDeriv * v)
-        return rDeriv.T * rDeriv
-
-    @Utils.timeIt
-    def _evalSmoothyy2Deriv(self, m, v=None):
-        if self.mrefInSmooth == True:
-            rDeriv = self.Wyy * ( self.mapping.deriv( m - self.mref ) )
-        elif self.mrefInSmooth == False:
-            rDeriv = self.Wyy * self.mapping.deriv(m)
-        if v is not None:
-            return rDeriv.T * (rDeriv * v)
-        return rDeriv.T * rDeriv
-
-    @Utils.timeIt
-    def _evalSmoothzz2Deriv(self, m, v=None):
-        if self.mrefInSmooth == True:
-            rDeriv = self.Wzz * ( self.mapping.deriv( m - self.mref ) )
-        elif self.mrefInSmooth == False:
-            rDeriv = self.Wzz * self.mapping.deriv(m)
-        if v is not None:
-            return rDeriv.T * (rDeriv * v)
-        return rDeriv.T * rDeriv
-
-    @Utils.timeIt
-    def _evalSmoothDeriv2(self, m):
-        deriv = self._evalSmoothxxDeriv(m)
-        if self.regmesh.dim > 1:
-            deriv += self._evalSmoothyyDeriv(m)
-        if self.regmesh.dim > 2:
-            deriv += self._evalSmoothzzDeriv(m)
-        return deriv
-
-    @Utils.timeIt
-    def _evalSmooth2Deriv2(self, m, v=None):
-        deriv = self._evalSmoothxx2Deriv(m, v)
-        if self.regmesh.dim > 1:
-            deriv += self._evalSmoothyy2Deriv(m, v)
-        if self.regmesh.dim > 2:
-            deriv += self._evalSmoothzz2Deriv(m, v)
-        return deriv
-
-
     @Utils.timeIt
     def eval(self, m):
-        return self._evalSmall(m) + self._evalSmooth(m) + self._evalSmooth2(m)
+        return self._evalSmall(m) + self._evalSmooth(m)
+
+    @Utils.timeIt
+    def _evalSmallDeriv(self,m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
+        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
+
+    @Utils.timeIt
+    def _evalSmoothDeriv(self,m):
+        if self.mrefInSmooth == True:
+            r = self.Wsmooth * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wsmooth * ( self.mapping * m )
+            return r.T * ( self.Wsmooth * self.mapping.deriv(m) )
 
     @Utils.timeIt
     def evalDeriv(self, m):
@@ -868,134 +568,184 @@ class Tikhonov(Simple):
             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
-        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m)
-
-    def eval2Deriv(self, m, v=None):
-        """
-        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})}
-
-        """
-        return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + self._evalSmooth2Deriv2(m, v)
+        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
 
 
-
-class Sparse(Simple):
+class Simple(Tikhonov):
     """
-        The regularization is:
-    
-        .. math::
-    
-            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R W(m-m_\\text{ref})}
-    
-        where the IRLS weight
-    
-        .. math::
-    
-            R = \eta TO FINISH LATER!!!
-    
-        So the derivative is straight forward:
-    
-        .. math::
-    
-            R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})}
-    
-        The IRLS weights are recomputed after each beta solves.
-        It is strongly recommended to do a few Gauss-Newton iterations
-        before updating.
+    Simple regularization that does not include length scales in the derivatives.
     """
-        
-    # set default values
-    eps_p = 1e-1        # Threshold value for the model norm
-    eps_q = 1e-1        # Threshold value for the model gradient norm
-    curModel = None     # Requires model to compute the weights
-    l2model = None
-    gamma = 1.          # Model norm scaling to smooth out convergence
-    norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz)
-    cell_weights = 1.        # Consider overwriting with sensitivity weights
+
+    mrefInSmooth = False  #: SMOOTH and SMOOTH_MOD_DIF options
+    alpha_s     = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
+    alpha_x     = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
+    alpha_y     = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
+    alpha_z     = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
+    wght        = 1.
 
     def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
 
-        if isinstance(self.cell_weights,float):
-            self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
+        if isinstance(self.wght,float):
+            self.wght = np.ones(self.regmesh.nC) * self.wght
 
     @property
     def Wsmall(self):
         """Regularization matrix Wsmall"""
         if getattr(self,'_Wsmall', None) is None:
-            if getattr(self, 'curModel', None) is None:
-                self.Rs = Utils.speye(self.regmesh.nC)
-
-            else:
-                f_m = self.mapping * (self.curModel - self.reg.mref)
-                self.rs = self.R(f_m , self.eps_p, self.norms[0])
-                self.Rs = Utils.sdiag( self.rs )
-
-            self._Wsmall = Utils.sdiag((self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs
-
+            self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
         return self._Wsmall
 
     @property
     def Wx(self):
         """Regularization matrix Wx"""
-        if getattr(self,'_Wx', None) is None:
-            if getattr(self, 'curModel', None) is None:
-                self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
-
-            else:
-                f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel)
-                self.rx = self.R( f_m , self.eps_q, self.norms[1])
-                self.Rx = Utils.sdiag( self.rx )
-
-            self._Wx = Utils.sdiag(( self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
-
+        if getattr(self, '_Wx', None) is None:
+            self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil
         return self._Wx
 
     @property
     def Wy(self):
         """Regularization matrix Wy"""
-        if getattr(self,'_Wy', None) is None:
-            if getattr(self, 'curModel', None) is None:
-                self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
-
-            else:
-                f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel)
-                self.ry = self.R( f_m , self.eps_q, self.norms[2])
-                self.Ry = Utils.sdiag( self.ry )
-
-            self._Wy = Utils.sdiag((self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
-
+        if getattr(self, '_Wy', None) is None:
+            self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil
         return self._Wy
 
     @property
     def Wz(self):
         """Regularization matrix Wz"""
-        if getattr(self,'_Wz', None) is None:
-            if getattr(self, 'curModel', None) is None:
-                self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
-
-            else:
-                f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel)
-                self.rz = self.R( f_m , self.eps_q, self.norms[3])
-                self.Rz = Utils.sdiag( self.rz )
-
-            self._Wz = Utils.sdiag((self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
-
+        if getattr(self, '_Wz', None) is None:
+            self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil
         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)
+
+
+class Sparse(Simple):
+
+    # set default values
+    eps_p      = 1e-1
+    eps_q      = 1e-1
+    curModel = None # use a model to compute the weights
+    gamma    = 1.
+    norms    = [0., 2., 2., 2.]
+    wght     = 1.
+
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+
+        if isinstance(self.wght,float):
+            self.wght = np.ones(self.regmesh.nC) * self.wght
+
+    @property
+    def Wsmall(self):
+        """Regularization matrix Wsmall"""
+        if getattr(self, 'curModel', None) is None:
+            self.Rs = Utils.speye(self.regmesh.nC)
+
+        else:
+            f_m = self.curModel - self.reg.mref
+            self.rs = self.R(f_m , self.eps_p, self.norms[0])
+            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
+            self.Rs = Utils.sdiag( self.rs )
+
+        return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
+
+
+    @property
+    def Wx(self):
+        """Regularization matrix Wx"""
+
+        if getattr(self, 'curModel', None) is None:
+            self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
+
+        else:
+            f_m = self.regmesh.cellDiffxStencil * self.curModel
+            self.rx = self.R( f_m , self.eps_q, self.norms[1])
+            self.Rx = Utils.sdiag( self.rx )
+
+        return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
+
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
+
+        if getattr(self, 'curModel', None) is None:
+            self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
+
+        else:
+            f_m = self.regmesh.cellDiffyStencil * self.curModel
+            self.ry = self.R( f_m , self.eps_q, self.norms[2])
+            self.Ry = Utils.sdiag( self.ry )
+
+        return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
+
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
+
+        if getattr(self, 'curModel', None) is None:
+            self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
+
+        else:
+            f_m = self.regmesh.cellDiffzStencil * self.curModel
+            self.rz = self.R( f_m , self.eps_q, self.norms[3])
+            self.Rz = Utils.sdiag( self.rz )
+
+        return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
+
+    @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 sp.vstack(wlist)
+
+    @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
+
     def R(self, f_m , eps, exponent):
 
-        # Eta scaling is important for mix-norms...do not mess with it
         eta = (eps**(1.-exponent/2.))**0.5
         r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)
 

SimPEG/Survey.py

@@ -1,5 +1,4 @@
 import Utils, numpy as np, scipy.sparse as sp, uuid
-import gc
 
 class BaseRx(object):
     """SimPEG Receiver Object"""

docs/api_Maps.rst

@@ -122,92 +122,10 @@ When these are used in the inverse problem, this is extremely important!!
 The API
 =======
 
-.. autoclass:: SimPEG.Maps.IdentityMap
+.. automodule:: SimPEG.Maps
     :members:
     :undoc-members:
 
 
-Common Maps
-===========
 
 
-Exponential Map
----------------
-
-Electrical conductivity varies over many orders of magnitude, so it is a common
-technique when solving the inverse problem to parameterize and optimize in terms
-of log conductivity. This makes sense not only because it ensures all conductivities
-will be positive, but because this is fundamentally the space where conductivity
-lives (i.e. it varies logarithmically).
-
-.. autoclass:: SimPEG.Maps.ExpMap
-    :members:
-    :undoc-members:
-
-
-Vertical 1D Map
----------------
-
-.. autoclass:: SimPEG.Maps.Vertical1DMap
-    :members:
-    :undoc-members:
-
-
-Map 2D Cross-Section to 3D Model
---------------------------------
-
-.. autoclass:: SimPEG.Maps.Map2Dto3D
-    :members:
-    :undoc-members:
-
-
-Mesh to Mesh Map
-----------------
-
-
-.. plot::
-
-    from SimPEG import *
-    import matplotlib.pyplot as plt
-    M = Mesh.TensorMesh([100,100])
-    h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
-    h1 = h1/h1.sum()
-    M2 = Mesh.TensorMesh([h1,h1])
-    V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
-    v = Utils.mkvc(V)
-    modh = Maps.Mesh2Mesh([M,M2])
-    modH = Maps.Mesh2Mesh([M2,M])
-    H = modH * v
-    h = modh * H
-    ax = plt.subplot(131)
-    M.plotImage(v, ax=ax)
-    ax.set_title('Fine Mesh (Original)')
-    ax = plt.subplot(132)
-    M2.plotImage(H,clim=[0,1],ax=ax)
-    ax.set_title('Course Mesh')
-    ax = plt.subplot(133)
-    M.plotImage(h,clim=[0,1],ax=ax)
-    ax.set_title('Fine Mesh (Interpolated)')
-    plt.show()
-
-
-.. autoclass:: SimPEG.Maps.Mesh2Mesh
-    :members:
-    :undoc-members:
-
-
-Some Extras
-===========
-
-Combo Map
----------
-
-The ComboMap holds the information for multiplying and combining
-maps. It also uses the chain rule to create the derivative.
-Remember, any time that you make your own combination of mappings
-be sure to test that the derivative is correct.
-
-.. autoclass:: SimPEG.Maps.ComboMap
-    :members:
-    :undoc-members:
-

docs/em/api_TDEM.rst

@@ -347,10 +347,10 @@ and
 
 
 
-TDEM Problem
-============
+TDEM - B formulation
+====================
 
-.. automodule:: SimPEG.EM.TDEM.TDEM
+.. automodule:: SimPEG.EM.TDEM.TDEM_b
     :show-inheritance:
     :members:
     :undoc-members:
@@ -359,7 +359,7 @@ TDEM Problem
 Field Storage
 =============
 
-.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.Fields
+.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.FieldsTDEM
     :show-inheritance:
     :members:
     :undoc-members:
@@ -369,19 +369,19 @@ Field Storage
 TDEM Survey Classes
 ===================
 
-.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.Survey
+.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.SurveyTDEM
     :show-inheritance:
     :members:
     :undoc-members:
     :inherited-members:
 
 
-.. Base Classes
-.. ============
+Base Classes
+============
 
-.. .. automodule:: SimPEG.EM.TDEM.BaseTDEM
-..     :show-inheritance:
-..     :members:
-..     :undoc-members:
-..     :inherited-members:
+.. automodule:: SimPEG.EM.TDEM.BaseTDEM
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:
 

docs/examples/DC_Forward_PseudoSection.rst

@@ -22,6 +22,7 @@ radi    = Radius of spheres [r1,r2]
 param   = Conductivity of background and two spheres [m0,m1,m2]
 surveyType = survey type 'pole-dipole' or 'dipole-dipole'
 unitType = Data type "appResistivity" | "appConductivity"  | "volt"
+
 Created by @fourndo
 
 

docs/examples/Forward_BasicDirectCurrent.rst

@@ -1,4 +1,4 @@
-.. _examples_Mesh_Basic_ForwardDC:
+.. _examples_Forward_BasicDirectCurrent:
 
 .. --------------------------------- ..
 ..                                   ..
@@ -8,18 +8,14 @@
 ..                                   ..
 .. --------------------------------- ..
 
-
-Mesh: Basic Forward 2D DC Resistivity
-=====================================
-
-2D DC forward modeling example with Tensor and Curvilinear Meshes
-
+Forward BasicDirectCurrent
+==========================
 
 .. plot::
 
     from SimPEG import Examples
-    Examples.Mesh_Basic_ForwardDC.run()
+    Examples.Forward_BasicDirectCurrent.run()
 
-.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_ForwardDC.py
+.. literalinclude:: ../../SimPEG/Examples/Forward_BasicDirectCurrent.py
     :language: python
     :linenos:

tests/base/test_maps.py

@@ -5,8 +5,10 @@ from scipy.sparse.linalg import dsolve
 
 TOL = 1e-14
 
-MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull","FullMap","Vertical1DMap"]
-MAPS_TO_TEST_3D = [             "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull","FullMap","Vertical1DMap"]
+MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull", "FullMap", "Vertical1DMap", "ParametrizedLayer", "ParametrizedBlockInLayer"]
+MAPS_TO_TEST_3D = [             "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull", "FullMap", "Vertical1DMap", "ParametrizedLayer", "ParametrizedBlockInLayer"]
+MAPS_TO_TEST_CYL = [             "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull", "FullMap", "Vertical1DMap", "ParametrizedLayer"]
+
 
 class MapTests(unittest.TestCase):
 
@@ -17,6 +19,8 @@ class MapTests(unittest.TestCase):
         self.mesh2 = Mesh.TensorMesh([a, b], x0=np.array([3, 5]))
         self.mesh3 = Mesh.TensorMesh([a, b, [3,4]], x0=np.array([3, 5, 2]))
         self.mesh22 = Mesh.TensorMesh([b, a], x0=np.array([3, 5]))
+        self.meshCyl = Mesh.CylMesh([10.,1.,10.], x0='00C')
+        print self.meshCyl._meshType
 
     def test_transforms2D(self):
         for M in MAPS_TO_TEST_2D:
@@ -28,6 +32,15 @@ class MapTests(unittest.TestCase):
             maps = getattr(Maps, M)(self.mesh3)
             self.assertTrue(maps.test())
 
+    def test_transformsCyl(self):
+        for M in MAPS_TO_TEST_CYL:
+            maps = getattr(Maps, M)(self.meshCyl)
+            self.assertTrue(maps.test())
+
+    def test_ParametricCasingAndLayer(self):
+        mapping = Maps.ParametrizedCasingAndLayer(self.meshCyl)
+        m = np.r_[-2., 1., 6., 2., -0.1, 0.2, 0.5, 0.2, -0.2, 0.2]
+        self.assertTrue(mapping.test(m))
 
     def test_transforms_logMap_reciprocalMap(self):
         # Note that log/reciprocal maps can be kinda finicky, so we are being explicit about the random seed.

tests/em/tdem/test_TDEM_b_DerivAdjoint.py

@@ -3,218 +3,310 @@ from SimPEG import *
 from SimPEG import EM
 
 plotIt = False
+tol = 1e-6
 
-testDeriv   = True
-testAdjoint = True
+class TDEM_bDerivTests(unittest.TestCase):
 
-TOL = 1e-5
+    def setUp(self):
 
-def setUp(prbtype='b', rxcomp='bz'):
-    cs = 5.
-    ncx = 20
-    ncy = 15
-    npad = 20
-    hx = [(cs,ncx), (cs,npad,1.3)]
-    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
-    mesh = Mesh.CylMesh([hx,1,hy], '00C')
-#
-    active = mesh.vectorCCz<0.
-    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
-    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
+        cs = 5.
+        ncx = 20
+        ncy = 6
+        npad = 20
+        hx = [(cs,ncx), (cs,npad,1.3)]
+        hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
+        mesh = Mesh.CylMesh([hx,1,hy], '00C')
 
-    rxOffset = 10.
-    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,]
-    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 0.]))
+        active = mesh.vectorCCz<0.
+        activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+        mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
 
-    survey = EM.TDEM.Survey([src])
+        rxOffset = 40.
+        rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
+        src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
 
-    if prbtype == 'b':
-        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
-    elif prbtype == 'e':
-        prb = EM.TDEM.Problem_e(mesh, mapping=mapping)
+        survey = EM.TDEM.SurveyTDEM([src])
 
-    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-    # prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)]
+        self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
+        # self.prb.timeSteps = [1e-5]
+        self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+        # self.prb.timeSteps = [(1e-05, 100)]
 
-    try:
-        from pymatsolver import MumpsSolver
-        prb.Solver = MumpsSolver
-    except ImportError, e:
-        prb.Solver = SolverLU
+        try:
+            from pymatsolver import MumpsSolver
+            self.prb.Solver = MumpsSolver
+        except ImportError, e:
+            self.prb.Solver = SolverLU
 
-    m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP)
+        self.sigma = np.ones(mesh.nCz)*1e-8
+        self.sigma[mesh.vectorCCz<0] = 1e-1
+        self.sigma = np.log(self.sigma[active])
 
-    prb.pair(survey)
-    mesh = mesh
+        self.prb.pair(survey)
+        self.mesh = mesh
 
-    return prb, m, mesh
+    def test_AhVec(self):
+        """
+            Test that fields and AhVec produce consistent results
+        """
+
+        prb = self.prb
+        sigma = self.sigma
+
+        u = prb.fields(sigma)
+        Ahu = prb._AhVec(sigma, u)
+
+        V1 = Ahu[:,'b',1]
+        V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0]
+        self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6)
+
+        V1 = Ahu[:,'e',1]
+        return np.linalg.norm(V1) < 1.e-6
+
+        for i in range(2,prb.nT):
+
+            dt = prb.timeSteps[i]
+
+            V1 = Ahu[:,'b',i]
+            V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1]
+            # print np.linalg.norm(V1), np.linalg.norm(V2)
+            self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
+
+            V1 = Ahu[:,'e',i]
+            V2 = prb.MeSigma*u[:,'e',i]
+            # print np.linalg.norm(V1), np.linalg.norm(V2)
+            return np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6
+
+    def test_AhVecVSMat_OneTS(self):
+
+        prb = self.prb
+        prb.timeSteps = [1e-05]
+        sigma = self.sigma
+        prb.curModel = sigma
+
+        dt = prb.timeSteps[0]
+        a11 = 1/dt*prb.MfMui*sp.identity(prb.mesh.nF)
+        a12 = prb.MfMui*prb.mesh.edgeCurl
+        a21 = prb.mesh.edgeCurl.T*prb.MfMui
+        a22 = -prb.MeSigma
+        A = sp.bmat([[a11,a12],[a21,a22]])
+
+        f = prb.fields(sigma)
+        u1 = A*f.tovec()
+        u2 = prb._AhVec(sigma,f).tovec()
+
+        self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
+
+    def test_solveAhVSMat_OneTS(self):
+        prb = self.prb
+
+        prb.timeSteps = [1e-05]
+
+        sigma = self.sigma
+        prb.curModel = sigma
+
+        dt = prb.timeSteps[0]
+        a11 = 1.0/dt*prb.MfMui*sp.identity(prb.mesh.nF)
+        a12 = prb.MfMui*prb.mesh.edgeCurl
+        a21 = prb.mesh.edgeCurl.T*prb.MfMui
+        a22 = -prb.MeSigma
+        A = sp.bmat([[a11,a12],[a21,a22]])
+
+        f = prb.fields(sigma)
+        f[:,:,0] = {'b':0}
+        f[:,'b',1] = 0
+
+        self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec()))
+
+        u1 = prb.solveAh(sigma,f).tovec().flatten()
+        u2 = sp.linalg.spsolve(A.tocsr(),f.tovec())
+
+        self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
+
+    def test_solveAhVsAhVec(self):
+
+        prb = self.prb
+        mesh = self.prb.mesh
+        sigma = self.sigma
+        self.prb.curModel = sigma
+
+        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
+        f[:,'b',:] = 0.0
+        for i in range(prb.nT):
+            f[:,'e', i] = np.random.rand(mesh.nE, 1)
+
+        Ahf = prb._AhVec(sigma, f)
+        f_test = prb.solveAh(sigma, Ahf)
+
+        u1 = f.tovec()
+        u2 = f_test.tovec()
+        self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
+
+    def test_DerivG(self):
+        """
+            Test the derivative of c with respect to sigma
+        """
+
+        # Random model and perturbation
+        sigma = np.random.rand(self.prb.mapping.nP)
+
+        f = self.prb.fields(sigma)
+        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()]
+        print '\ntest_DerivG'
+        passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
+        return passed
+
+    def test_Deriv_dUdM(self):
+
+        prb = self.prb
+        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
+        mesh = self.mesh
+        sigma = self.sigma
+
+        dm = 10*np.random.rand(prb.mapping.nP)
+        f = prb.fields(sigma)
+
+        derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
+        print '\n'
+        print 'test_Deriv_dUdM'
+        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
+
+    def test_Deriv_J(self):
+
+        prb = self.prb
+        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
+        mesh = self.mesh
+        sigma = self.sigma
+
+        # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
+        d_sig = 10*np.random.rand(prb.mapping.nP)
 
 
-class TDEM_DerivTests(unittest.TestCase):
+        derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
+        print '\n'
+        print 'test_Deriv_J'
+        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
 
-# ====== TEST A ========== #
+    def test_projectAdjoint(self):
+        prb = self.prb
+        survey = prb.survey
+        mesh = self.mesh
 
-    def AderivTest(self, prbtype):
-        prb, m0, mesh = setUp(prbtype)
-        tInd = 2
-        if prbtype == 'b':
-            nu = mesh.nF
-        elif prbtype == 'e':
-            nu = mesh.nE
-        v = np.random.rand(nu)
+        # Generate random fields and data
+        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
+        for i in range(prb.nT):
+            f[:,'b',i] = np.random.rand(mesh.nF, 1)
+            f[:,'e',i] = np.random.rand(mesh.nE, 1)
+        d_vec = np.random.rand(survey.nD)
+        d = Survey.Data(survey,v=d_vec)
 
-        def AderivFun(m):
-            prb.curModel = m
-            A = prb.getAdiag(tInd)
-            Av = A*v
-            prb.curModel = m0
-            ADeriv_dm = lambda dm: prb.getAdiagDeriv(tInd, v, dm)
+        # Check that d.T*Q*f = f.T*Q.T*d
+        V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec())
+        V2 = f.tovec().dot(survey.evalDeriv(None, v=d, adjoint=True).tovec())
 
-            return Av, ADeriv_dm
+        self.assertTrue((V1-V2)/np.abs(V1) < tol)
 
-        print '\n Testing ADeriv %s'%(prbtype)
-        Tests.checkDerivative(AderivFun, m0, plotIt=False, num=4, eps=1e-20)
+    def test_adjointAhVsAht(self):
+        prb = self.prb
+        mesh = self.mesh
+        sigma = self.sigma
 
-    def A_adjointTest(self,prbtype):
-        prb, m0, mesh = setUp(prbtype)
-        tInd = 2
+        f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
+        for i in range(1,prb.nT+1):
+            f1[:,'b',i] = np.random.rand(mesh.nF, 1)
+            f1[:,'e',i] = np.random.rand(mesh.nE, 1)
 
-        print '\n Testing A_adjoint'
-        m = np.random.rand(prb.mapping.nP)
-        if prbtype == 'b':
-            nu = prb.mesh.nF
-        elif prbtype == 'e':
-            nu = prb.mesh.nE
+        f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
+        for i in range(1,prb.nT+1):
+            f2[:,'b',i] = np.random.rand(mesh.nF, 1)
+            f2[:,'e',i] = np.random.rand(mesh.nE, 1)
 
-        v = np.random.rand(nu)
-        u = np.random.rand(nu)
-        prb.curModel = m0
+        V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec())
+        V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec())
+        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
 
-        tInd = 2 # not actually used
-        V1 = v.dot(prb.getAdiagDeriv(tInd, u, m))
-        V2 = m.dot(prb.getAdiagDeriv(tInd, u, v, adjoint=True))
-        passed = np.abs(V1-V2) < TOL * (np.abs(V1) + np.abs(V2))/2.
-        print 'AdjointTest %s'%(prbtype), V1, V2, passed
-        self.assertTrue(passed)
+    # def test_solveAhtVsAhtVec(self):
+    #     prb = self.prb
+    #     mesh = self.mesh
+    #     sigma = np.random.rand(prb.mapping.nP)
 
-    def test_Aderiv_b(self):
-        self.AderivTest('b')
-    def test_Aderiv_e(self):
-        self.AderivTest('e')
+    #     f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey)
+    #     for i in range(1,prb.nT+1):
+    #         f1[:,'b',i] = np.random.rand(mesh.nF, 1)
+    #         f1[:,'e',i] = np.random.rand(mesh.nE, 1)
 
-    def test_Aadjoint_b(self):
-        self.A_adjointTest('b')
-    def test_Aadjoint_e(self):
-        self.A_adjointTest('e')
+    #     f2 = prb.solveAht(sigma, f1)
+    #     f3 = prb._AhtVec(sigma, f2)
 
-# ====== TEST Fields Deriv Pieces ========== #
+    #     if True:
+    #         import matplotlib.pyplot as plt
+    #         plt.plot(f3.tovec(),'b')
+    #         plt.plot(f1.tovec(),'r')
+    #         plt.show()
+    #     V1 = np.linalg.norm(f3.tovec()-f1.tovec())
+    #     V2 = np.linalg.norm(f1.tovec())
+    #     print 'AhtVsAhtVec', V1, V2, f1.tovec()
+    #     print 'I am gunna fail this one: boo. :('
+    #     self.assertLess(V1/V2, 1e-6)
 
-    def test_eDeriv_m_adjoint(self):
-        prb, m0, mesh = setUp()
-        tInd = 0
+    # def test_adjointsolveAhVssolveAht(self):
+    #     prb = self.prb
+    #     mesh = self.mesh
+    #     sigma = self.sigma
 
-        v = np.random.rand(mesh.nF)
+    #     f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
+    #     for i in range(1,prb.nT+1):
+    #         f1[:,'b',i] = np.random.rand(mesh.nF, 1)
+    #         f1[:,'e',i] = np.random.rand(mesh.nE, 1)
 
-        print '\n Testing eDeriv_m Adjoint'
+    #     f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
+    #     for i in range(1,prb.nT+1):
+    #         f2[:,'b',i] = np.random.rand(mesh.nF, 1)
+    #         f2[:,'e',i] = np.random.rand(mesh.nE, 1)
 
-        prb, m0, mesh = setUp()
-        f = prb.fields(m0)
+    #     V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
+    #     V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
+    #     print V1, V2
+    #     self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
+
+    def test_adjointGvecVsGtvec(self):
+        mesh = self.mesh
+        prb = self.prb
 
         m = np.random.rand(prb.mapping.nP)
-        e = np.random.randn(prb.mesh.nE)
-        V1 = e.dot(f._eDeriv_m(1, prb.survey.srcList[0], m))
-        V2 = m.dot(f._eDeriv_m(1, prb.survey.srcList[0], e, adjoint=True))
-        tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
-        passed = np.abs(V1-V2) < tol
+        sigma = np.random.rand(prb.mapping.nP)
 
-        print '    ', V1, V2, np.abs(V1-V2), tol, passed
+        u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
+        for i in range(1,prb.nT+1):
+            u[:,'b',i] = np.random.rand(mesh.nF, 1)
+            u[:,'e',i] = np.random.rand(mesh.nE, 1)
+
+        v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
+        for i in range(1,prb.nT+1):
+            v[:,'b',i] = np.random.rand(mesh.nF, 1)
+            v[:,'e',i] = np.random.rand(mesh.nE, 1)
+
+        V1 = m.dot(prb.Gtvec(sigma, v, u))
+        V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
+        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
+
+    def test_adjointJvecVsJtvec(self):
+        mesh = self.mesh
+        prb = self.prb
+        sigma = self.sigma
+
+        m = np.random.rand(prb.mapping.nP)
+        d = np.random.rand(prb.survey.nD)
+
+        V1 = d.dot(prb.Jvec(sigma, m))
+        V2 = m.dot(prb.Jtvec(sigma, d))
+        passed = np.abs(V1-V2)/np.abs(V1) < tol
+        print 'AdjointTest', V1, V2, passed
         self.assertTrue(passed)
 
-    def test_eDeriv_u_adjoint(self):
-        print '\n Testing eDeriv_u Adjoint'
-
-        prb, m0, mesh = setUp()
-        f = prb.fields(m0)
-
-        b = np.random.rand(prb.mesh.nF)
-        e = np.random.randn(prb.mesh.nE)
-        V1 = e.dot(f._eDeriv_u(1, prb.survey.srcList[0], b))
-        V2 = b.dot(f._eDeriv_u(1, prb.survey.srcList[0], e, adjoint=True))
-        tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
-        passed = np.abs(V1-V2) < tol
-
-        print '    ', V1, V2, np.abs(V1-V2), tol, passed
-        self.assertTrue(passed)
-
-
-# ====== TEST Jvec ========== #
-
-    if testDeriv:
-
-        def JvecTest(self, prbtype, rxcomp):
-            prb, m, mesh = setUp(prbtype, rxcomp)
-
-            derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)]
-            print '\n'
-            print 'test_Jvec_%s_%s' %(prbtype, rxcomp)
-            Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20)
-
-        def test_Jvec_b_bx(self):
-            self.JvecTest('b','bx')
-
-        def test_Jvec_b_bz(self):
-            self.JvecTest('b','bz')
-
-        def test_Jvec_b_dbxdt(self):
-            self.JvecTest('b','dbxdt')
-
-        def test_Jvec_b_dbzdt(self):
-            self.JvecTest('b','dbzdt')
-
-        def test_Jvec_b_ey(self):
-            self.JvecTest('b','ey')
-
-        def test_Jvec_e_ey(self):
-            self.JvecTest('e','ey')
-
-
-# ====== TEST Jtvec ========== #
-
-    if testAdjoint:
-
-        def JvecVsJtvecTest(self, prbtype='b', rxcomp='bz'):
-
-            print '\nAdjoint Testing Jvec, Jtvec %s' %(rxcomp)
-
-            prb, m0, mesh = setUp(prbtype, rxcomp)
-            m = np.random.rand(prb.mapping.nP)
-            d = np.random.randn(prb.survey.nD)
-            V1 = d.dot(prb.Jvec(m0, m))
-            V2 = m.dot(prb.Jtvec(m0, d))
-            tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
-            passed = np.abs(V1-V2) < tol
-
-            print '    ', V1, V2, np.abs(V1-V2), tol, passed
-            self.assertTrue(passed)
-
-        def test_Jvec_adjoint_b_bx(self):
-            self.JvecVsJtvecTest('b', 'bx')
-
-        def test_Jvec_adjoint_b_bz(self):
-            self.JvecVsJtvecTest('b', 'bz')
-
-        def test_Jvec_adjoint_b_dbxdt(self):
-            self.JvecVsJtvecTest('b', 'bx')
-
-        def test_Jvec_adjoint_b_dbzdt(self):
-            self.JvecVsJtvecTest('b', 'bz')
-
-        def test_Jvec_adjoint_b_ey(self):
-            self.JvecVsJtvecTest('b', 'ey')
-
-        # This is not working because Problem_e has not done
-        # def test_Jvec_adjoint_e_ey(self):
-        #     self.JvecVsJtvecTest('e', 'ey')
 
 
 

tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py

@@ -3,12 +3,10 @@ from SimPEG import *
 from SimPEG import EM
 
 plotIt = False
-testDeriv = True
-testAdjoint = True
 
-TOL = 1e-5
+class TDEM_bDerivTests(unittest.TestCase):
 
-def setUp(self, rxcomp='bz'):
+    def setUp(self):
 
         cs = 5.
         ncx = 20
@@ -23,78 +21,131 @@ def setUp(self, rxcomp='bz'):
         mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
 
         rxOffset = 40.
-        rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), rxcomp)
-        src = EM.TDEM.Src.MagDipole( [rx], loc=np.array([0., 0., 0.]))
-        rx2 = EM.TDEM.Rx(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), rxcomp)
-        src2 = EM.TDEM.Src.MagDipole( [rx2], loc=np.array([0., 0., 0.]))
+        rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
+        src = EM.TDEM.SrcTDEM_VMD_MVP( [rx], loc=np.array([0., 0., 0.]))
+        rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz')
+        src2 = EM.TDEM.SrcTDEM_VMD_MVP( [rx2], loc=np.array([0., 0., 0.]))
 
-        survey = EM.TDEM.Survey([src,src2])
+        survey = EM.TDEM.SurveyTDEM([src,src2])
 
-        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
-        # prb.timeSteps = [1e-5]
-        prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-        # prb.timeSteps = [(1e-05, 100)]
+        self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
+        # self.prb.timeSteps = [1e-5]
+        self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+        # self.prb.timeSteps = [(1e-05, 100)]
 
         try:
             from pymatsolver import MumpsSolver
-            prb.Solver = MumpsSolver
+            self.prb.Solver = MumpsSolver
         except ImportError, e:
-            prb.Solver  = SolverLU
+            self.prb.Solver  = SolverLU
 
-        m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP)
+        self.sigma = np.ones(mesh.nCz)*1e-8
+        self.sigma[mesh.vectorCCz<0] = 1e-1
+        self.sigma = np.log(self.sigma[active])
 
-        prb.pair(survey)
+        self.prb.pair(survey)
+        self.mesh = mesh
 
-        return mesh, prb, m
+    def test_DerivG(self):
+        """
+            Test the derivative of c with respect to sigma
+        """
 
-class TDEM_bDerivTests(unittest.TestCase):
+        # Random model and perturbation
+        sigma = np.random.rand(self.prb.mapping.nP)
+
+        f = self.prb.fields(sigma)
+        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()]
+        print '\ntest_DerivG'
+        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
+
+    def test_Deriv_dUdM(self):
+
+        prb = self.prb
+        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
+        mesh = self.mesh
+        sigma = self.sigma
+
+        dm = 10*np.random.rand(prb.mapping.nP)
+        f = prb.fields(sigma)
+
+        derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
+        print '\n'
+        print 'test_Deriv_dUdM'
+        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
+
+    def test_Deriv_J(self):
+
+        prb = self.prb
+        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
+        mesh = self.mesh
+        sigma = self.sigma
+
+        # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
+        d_sig = 10*np.random.rand(prb.mapping.nP)
 
 
-    if testDeriv:
-        def Deriv_J(self, rxcomp='bz'):
+        derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
+        print '\n'
+        print 'test_Deriv_J'
+        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
 
-            mesh, prb, m0 = setUp(rxcomp)
+    def test_projectAdjoint(self):
+        prb = self.prb
+        survey = prb.survey
+        nSrc = survey.nSrc
+        mesh = self.mesh
 
-            prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
+        # Generate random fields and data
+        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
+        for i in range(prb.nT):
+            f[:,'b',i] = np.random.rand(mesh.nF, nSrc)
+            f[:,'e',i] = np.random.rand(mesh.nE, nSrc)
+        d_vec = np.random.rand(survey.nD)
+        d = Survey.Data(survey,v=d_vec)
 
-            derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m0, mx)]
-            print '\n'
-            print 'test_Deriv_J %s'%rxcomp
-            Tests.checkDerivative(derChk, m0, plotIt=False, num=3, eps=1e-20)
+        # Check that d.T*Q*f = f.T*Q.T*d
+        V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec())
+        V2 = np.sum((f.tovec())*(survey.evalDeriv(None, v=d, adjoint=True).tovec()))
 
-        def test_Jvec_bx(self):
-            self.Deriv_J('bx')
+        self.assertTrue((V1-V2)/np.abs(V1) < 1e-6)
 
-        def test_Jvec_bz(self):
-            self.Deriv_J('bz')
+    def test_adjointGvecVsGtvec(self):
+        mesh = self.mesh
+        prb = self.prb
 
-        def test_Jvec_ey(self):
-            self.Deriv_J('ey')
+        m = np.random.rand(prb.mapping.nP)
+        sigma = np.random.rand(prb.mapping.nP)
 
-    if testAdjoint:
-        def adjointJvecVsJtvec(self, rxcomp='bz'):
-            print ' \n Testing Adjoint %s' %rxcomp
-            mesh, prb, m0 = setUp(rxcomp)
+        u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
+        for i in range(1,prb.nT+1):
+            u[:,'b',i] = np.random.rand(mesh.nF, 2)
+            u[:,'e',i] = np.random.rand(mesh.nE, 2)
 
-            m = np.random.rand(prb.mapping.nP)
-            d = np.random.rand(prb.survey.nD)
+        v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
+        for i in range(1,prb.nT+1):
+            v[:,'b',i] = np.random.rand(mesh.nF, 2)
+            v[:,'e',i] = np.random.rand(mesh.nE, 2)
 
-            V1 = d.dot(prb.Jvec(m0, m))
-            V2 = m.dot(prb.Jtvec(m0, d))
+        V1 = m.dot(prb.Gtvec(sigma, v, u))
+        V2 = np.sum(v.tovec()*prb.Gvec(sigma, m, u).tovec())
+        self.assertTrue(np.abs(V1-V2)/np.abs(V1) <1e-6)
 
-            tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
-            passed = np.abs(V1-V2) < tol
-            print '    ', V1, V2, np.abs(V1-V2), tol, passed
-            self.assertTrue(passed)
+    def test_adjointJvecVsJtvec(self):
+        mesh = self.mesh
+        prb = self.prb
+        sigma = self.sigma
 
-        def test_JvecVsJtvec_bx(self):
-            self.adjointJvecVsJtvec('bx')
+        m = np.random.rand(prb.mapping.nP)
+        d = np.random.rand(prb.survey.nD)
 
-        def test_JvecVsJtvec_bz(self):
-            self.adjointJvecVsJtvec('bz')
-
-        def test_JvecVsJtvec_ey(self):
-            self.adjointJvecVsJtvec('ey')
+        V1 = d.dot(prb.Jvec(sigma, m))
+        V2 = m.dot(prb.Jtvec(sigma, d))
+        print 'AdjointTest', V1, V2
+        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < 1e-6)
 
 
 

tests/em/tdem/test_TDEM_combos.py

@@ -0,0 +1,94 @@
+import unittest
+from SimPEG import *
+from SimPEG import EM
+
+plotIt = False
+
+def getProb(meshType='CYL',rxTypes='bx,bz',nSrc=1):
+    cs = 5.
+    ncx = 20
+    ncy = 6
+    npad = 20
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hy], '00C')
+
+    active = mesh.vectorCCz<0.
+    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
+
+    rxOffset = 40.
+
+    srcs = []
+    for ii in range(nSrc):
+        rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')]
+        srcs += [EM.TDEM.SrcTDEM_VMD_MVP(rxs,np.array([0., 0., 0.]))]
+
+    survey = EM.TDEM.SurveyTDEM(srcs)
+
+    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
+    # prb.timeSteps = [1e-5]
+    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+    # prb.timeSteps = [(1e-05, 100)]
+
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+        prb.Solver  = SolverLU
+
+    sigma = np.ones(mesh.nCz)*1e-8
+    sigma[mesh.vectorCCz<0] = 1e-1
+    sigma = np.log(sigma[active])
+
+    prb.pair(survey)
+    return prb, mesh, sigma
+
+def dotestJvec(prb, mesh, sigma):
+    prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
+    # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
+    d_sig = 10*np.random.rand(prb.mapping.nP)
+    derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
+    return Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20)
+
+def dotestAdjoint(prb, mesh, sigma):
+    m = np.random.rand(prb.mapping.nP)
+    d = np.random.rand(prb.survey.nD)
+
+    V1 = d.dot(prb.Jvec(sigma, m))
+    V2 = m.dot(prb.Jtvec(sigma, d))
+    print 'AdjointTest', V1, V2
+    return np.abs(V1-V2)/np.abs(V1), 1e-6
+
+class TDEM_bDerivTests(unittest.TestCase):
+
+    def test_Jvec_bx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx')))
+    def test_Adjoint_bx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx')))
+
+    def test_Jvec_bxbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz')))
+    def test_Adjoint_bxbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz')))
+
+    def test_Jvec_bxbz_2src(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nSrc=2)))
+    def test_Adjoint_bxbz_2src(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nSrc=2)))
+
+    def test_Jvec_bxbzbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz,bz')))
+    def test_Adjoint_bxbzbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz,bz')))
+
+    def test_Jvec_dbxdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbxdt')))
+    def test_Adjoint_dbxdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbxdt')))
+
+    def test_Jvec_dbzdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbzdt')))
+    def test_Adjoint_dbzdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbzdt')))
+
+    def test_Jvec_dbxdtbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbxdt,bz')))
+    def test_Adjoint_dbxdtbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbxdt,bz')))
+
+    def test_Jvec_ey(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey')))
+    def test_Adjoint_ey(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey')))
+
+    def test_Jvec_eybzdbxdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey,bz,dbxdt')))
+    def test_Adjoint_eybzdbxdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey,bz,dbxdt')))
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/em/tdem/test_TDEM_crosscheck.py

@@ -1,76 +0,0 @@
-import unittest
-from SimPEG import *
-from SimPEG import EM
-
-TOL = 1e-5
-FLR = 1e-20
-
-np.random.seed(seed=25) # set a seed so that the same conductivity model is used for all runs
-
-def setUp(prbtype = 'b', rxcomp='bz'):
-    cs = 5.
-    ncx = 20
-    ncy = 15
-    npad = 20
-    hx = [(cs,ncx), (cs,npad,1.3)]
-    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
-    mesh = Mesh.CylMesh([hx,1,hy], '00C')
-#
-    active = mesh.vectorCCz<0.
-    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
-    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
-
-    rxOffset = 10.
-    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,]
-    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 0.]))
-
-    survey = EM.TDEM.Survey([src])
-
-    if prbtype == 'b':
-        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
-    elif prbtype == 'e':
-        prb = EM.TDEM.Problem_e(mesh, mapping=mapping)
-
-    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-    # prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)]
-
-    try:
-        from pymatsolver import MumpsSolver
-        prb.Solver = MumpsSolver
-    except ImportError, e:
-        prb.Solver = SolverLU
-
-    m = np.log(1e-1)*np.ones(prb.mapping.nP) #+ 1e-2*np.random.randn(prb.mapping.nP)
-
-    prb.pair(survey)
-    mesh = mesh
-
-    return prb, m, mesh
-
-def CrossCheck(prbtype1='b', prbtype2='e', rxcomp='bz'):
-
-    prb1,m1,mesh1 = setUp(prbtype1, rxcomp)
-    prb2,m2,mesh2 = setUp(prbtype2, rxcomp)
-
-    assert (m1 == m2).all(), 'Models for two formulations are different'
-
-    d1 = prb1.survey.dpred(m1)
-    d2 = prb2.survey.dpred(m2)
-
-
-    check = np.linalg.norm(d1 - d2)
-    tol   = 0.5 * (np.linalg.norm(d1) + np.linalg.norm(d2)) * TOL
-    passed = check < tol
-
-    print 'Checking %s, %s for %s data'%(prbtype1, prbtype2, rxcomp)
-    print '  ', np.linalg.norm(d1), np.linalg.norm(d2), np.linalg.norm(check), tol, passed
-
-    assert passed
-
-class TDEM_cross_check_EB(unittest.TestCase):
-    def test_EB_ey(self):
-        CrossCheck('b','e','ey')
-
-if __name__ == '__main__':
-    unittest.main()
-

tests/em/tdem/test_TDEM_forward_Analytic.py

@@ -29,12 +29,12 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
     actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
     mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
 
-    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
-    src = EM.TDEM.Src.MagDipole([rx], waveform= EM.TDEM.Src.StepOffWaveform(), loc=np.array([0., 0., 0.]))
+    rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
+    src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
     # src = EM.TDEM.SrcTDEM([rx], loc=np.array([0., 0., 0.]))
 
-    survey = EM.TDEM.Survey([src])
-    prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+    survey = EM.TDEM.SurveyTDEM([src])
+    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
     prb.Solver = MumpsSolver
 
     prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40), (0.0001, 40), (0.0005, 40)]
@@ -50,8 +50,6 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
 
     ind = np.logical_and(rx.times > bounds[0],rx.times < bounds[1])
     log10diff = np.linalg.norm(np.log10(np.abs(bz_calc[ind])) - np.log10(np.abs(bz_ana[ind])))/np.linalg.norm(np.log10(np.abs(bz_ana[ind])))
-
-    print ' |bz_ana| = ',np.linalg.norm(bz_ana), ' |bz_num| = ', np.linalg.norm(bz_calc), ' |bz_ana - bz_num| =', np.linalg.norm(bz_ana-bz_calc)
     print 'Difference: ', log10diff
 
     if showIt == True:
@@ -63,12 +61,6 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
     return log10diff
 
 
-class TDEM_SimpleSrcTests(unittest.TestCase):
-    def test_source(self):
-        waveform = EM.TDEM.Src.StepOffWaveform()
-        assert waveform.eval(0.) == 0.
-
-
 class TDEM_bTests(unittest.TestCase):
 
     def test_analytic_p2_CYL_50m(self):