Merge pull request #339 from simpeg/em/secondary_Rx

Implementation for Inverting Secondary B field

SimPEG/Directives.py

@@ -144,6 +144,7 @@ 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.
@@ -242,12 +243,6 @@ 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'
@@ -258,56 +253,138 @@ 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):
 
-        # 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
+        if self.mode == 1:
+            self.reg.norms = [2., 2., 2., 2.]
 
     def endIter(self):
-        # Cool the threshold parameter if required
-        if getattr(self, 'factor', None) is not None:
-            eps = self.reg.eps / self.factor
 
-            if getattr(self, 'eps_min', None) is not None:
-                self.reg.eps = np.max([self.eps_min,eps])
+        # 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
             else:
-                self.reg.eps = eps
+                self.f_old = phim_new
 
-        # Get phi_m at the end of current iteration
-        self.phi_m_last = self.invProb.phi_m_last
+            # Cool the threshold parameter if required
+            if getattr(self, 'factor', None) is not None:
+                eps = self.reg.eps / self.factor
 
-         # Update the model used for the IRLS weights
-        self.reg.curModel = self.invProb.curModel
+                if getattr(self, 'eps_min', None) is not None:
+                    self.reg.eps = np.max([self.eps_min,eps])
+                else:
+                    self.reg.eps = eps
 
-        # Temporarely set gamma to 1. to get raw phi_m
-        self.reg.gamma = 1.
+            # Get phi_m at the end of current iteration
+            self.phi_m_last = self.invProb.phi_m_last
 
-        # Compute new model objective function value
-        phim_new = self.reg.eval(self.invProb.curModel)
+            # 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
 
-        # Update gamma to scale the regularization between IRLS iterations
-        self.reg.gamma = self.phi_m_last / phim_new
+             # Update the model used for the IRLS weights
+            self.reg.curModel = self.invProb.curModel
 
-        # Set the weighting matrix to None so that it is recomputed next time
-        # it is called in the inversion
-        self.reg._W = None
+            # 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
 
 class Update_lin_PreCond(InversionDirective):
     """
@@ -360,19 +437,3 @@ 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,6 +60,20 @@ 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
@@ -124,6 +138,21 @@ 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
@@ -257,7 +286,7 @@ class Fields3D_e(Fields):
         """
 
         # assuming primary does not depend on the model
-        return src.ePrimaryDeriv(self.prob, v, adjoint) #Zero()
+        return Zero()
 
     def _bPrimary(self, eSolution, srcList):
         """
@@ -471,6 +500,8 @@ 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:
@@ -600,8 +631,8 @@ class Fields3D_b(Fields):
 
 
         if adjoint:
-            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)
+            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv
+        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv
 
     def _j(self, bSolution, srcList):
         """

SimPEG/EM/FDEM/ProblemFDEM.py

@@ -74,8 +74,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         self.curModel = m
 
-        # Jv = self.dataPair(self.survey)
-        Jv = []
+        Jv = self.dataPair(self.survey)
 
         for freq in self.survey.freqs:
             A = self.getA(freq)
@@ -90,9 +89,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.append(rx.evalDeriv(src, self.mesh, f, df_dm_v))
+                    Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
             Ainv.clean()
-        return np.hstack(Jv)
+        return Utils.mkvc(Jv)
 
     def Jtvec(self, m, v, f=None):
         """
@@ -167,6 +166,7 @@ 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,6 +97,19 @@ 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,18 +60,6 @@ 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
@@ -84,18 +72,6 @@ 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
@@ -108,18 +84,6 @@ 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
@@ -132,18 +96,6 @@ 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
@@ -603,7 +555,7 @@ class CircularLoop(BaseSrc):
             a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
-            srcfct = MagneticLoopVectorPotential
+            srcfct = MagneticDipoleVectorPotential
             ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
             ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
             az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)
@@ -662,155 +614,5 @@ 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,14 +43,7 @@ class BaseRx(SimPEG.Survey.BaseRx):
         elif adjoint:
             return P.T*v
 
-# DC.Rx.Pole(locs)
-class Pole(BaseRx):
-
-    def __init__(self, locs, rxType = 'phi', **kwargs):
-        BaseRx.__init__(self, locs, rxType)
-
-
-# DC.Rx.Dipole(locsM, locsN)
+# DC.Rx.Dipole(locs)
 class Dipole(BaseRx):
 
     def __init__(self, locsM, locsN, rxType = 'phi', **kwargs):

SimPEG/Examples/Inversion_IRLS.py

@@ -1,7 +1,7 @@
 from SimPEG import *
 
 
-def run(N=200, plotIt=True):
+def run(N=100, plotIt=True):
     """
         Inversion: Linear Problem
         =========================
@@ -18,6 +18,8 @@ def run(N=200, 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.
@@ -50,57 +52,47 @@ def run(N=200, plotIt=True):
     wr = np.sum(prob.G**2.,axis=0)**0.5
     wr = ( wr/np.max(wr) )
 
-    reg = Regularization.Simple(mesh)
-    reg.wght = wr
-
+#    reg = Regularization.Simple(mesh)
+#    reg.mref = mref
+#    reg.cell_weights = wr
+#
     dmis = DataMisfit.l2_DataMisfit(survey)
     dmis.Wd = 1./wd
-
-    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()
-
+#
+#    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()
+#
     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)
-    reg.eps_p = 5e-2
-    reg.eps_q = 1e-2
-    reg.norms   = [0., 0., 2., 2.]
-    reg.wght = wr
+    eps_p = 5e-2
+    eps_q = 5e-2
+    norms   = [0., 0., 2., 2.]
 
-    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 )
+    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)
 
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
-
-    m0 = mrec
+    inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
 
     # Run inversion
     mrec = inv.run(m0)
@@ -117,7 +109,7 @@ def run(N=200, plotIt=True):
         axes[0].set_title('Columns of matrix G')
 
         axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
-        axes[1].plot(mesh.vectorCCx, ml2, 'r-')
+        axes[1].plot(mesh.vectorCCx, reg.l2model, 'r-')
         #axes[1].legend(('True Model', 'Recovered Model'))
         axes[1].set_ylim(-1.0,1.25)
 

SimPEG/Examples/Mesh_Basic_ForwardDC.py

@@ -1,22 +1,25 @@
 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(1./sigma)
-        MsigI = Utils.sdInv(Msigi)
-        A = D*MsigI*G
+        MsigI = mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True)
+        A = -D*MsigI*D.T
         A[-1,-1] /= mesh.vol[-1] # Remove null space
         rhs = np.zeros(mesh.nC)
         txind = Utils.meshutils.closestPoints(mesh, pts)
@@ -37,39 +40,17 @@ 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)
-
-    #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)
-
+    dat = rM.plotImage(phirM, ax=axes[1], clim=(vmin, vmax), grid=True)
     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/__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", "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"]
+__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"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Maps.py

@@ -1,4 +1,3 @@
-from __future__ import division
 import Utils, numpy as np, scipy.sparse as sp
 from scipy.sparse.linalg import LinearOperator
 from Tests import checkDerivative
@@ -6,7 +5,6 @@ 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):
     """
@@ -19,7 +17,7 @@ class IdentityMap(object):
         Utils.setKwargs(self, **kwargs)
 
         if nP is not None:
-            assert type(nP) in [int, long, np.int64], ' Number of parameters must be an integer.'
+            assert type(nP) in [int, long], ' Number of parameters must be an integer.'
 
         self.mesh = mesh
         self._nP  = nP
@@ -131,15 +129,7 @@ class IdentityMap(object):
 
 
 class ComboMap(IdentityMap):
-    """
-        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.
-
-    """
+    """Combination of various maps."""
 
     def __init__(self, maps, **kwargs):
         IdentityMap.__init__(self, None, **kwargs)
@@ -188,12 +178,6 @@ 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
@@ -465,32 +449,6 @@ 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):
@@ -543,19 +501,11 @@ 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:
-            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.')
+            self.valInactive = valInactive.copy()
+        self.valInactive[self.indActive] = 0
 
         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
@@ -624,37 +574,6 @@ 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
@@ -697,13 +616,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]
 
@@ -1056,697 +975,7 @@ 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,6 +2,7 @@ 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
@@ -10,7 +11,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):
+class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView):
     """
     CurvilinearMesh is a mesh class that deals with curvilinear meshes.
 
@@ -330,102 +331,6 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
 
 
 
-    #############################################
-    #            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,7 +552,8 @@ class CurvView(object):
     def __init__(self):
         pass
 
-    def plotGrid(self, length=0.05, showIt=False):
+
+    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.
 
 
@@ -560,60 +561,63 @@ class CurvView(object):
             :include-source:
 
             from SimPEG import Mesh, Utils
-            X, Y = Utils.exampleCurvGird([3,3],'rotate')
+            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
+
+        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()
 
-            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()
+            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()
 
-            X = np.r_[X1, X2]
-            Y = np.r_[Y1, Y2]
+                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()
 
-            plt.plot(X, Y)
+                X = np.r_[X1, X2]
+                Y = np.r_[Y1, Y2]
 
-            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')
+                ax.plot(X, Y, 'b-')
+            if centers:
+                ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
 
-            plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
+            # 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 = 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-')
+            # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
 
-            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.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-')
 
-            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-')
+            # 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-')
 
-            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')
+            # 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:
-            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()
@@ -630,16 +634,50 @@ class CurvView(object):
             Y = np.r_[Y1, Y2, Y3]
             Z = np.r_[Z1, Z2, Z3]
 
-            plt.plot(X, Y, 'b', zs=Z)
+            ax.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,4 +1,6 @@
-import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
+import Utils, Maps, Mesh
+import numpy as np
+import scipy.sparse as sp
 
 class RegularizationMesh(object):
     """
@@ -39,7 +41,7 @@ class RegularizationMesh(object):
             if self.indActive is None:
                 self._nC = self.mesh.nC
             else:
-                self._nC = int(sum(self.indActive))
+                self._nC = sum(self.indActive)
         return self._nC
 
     @property
@@ -304,7 +306,7 @@ class BaseRegularization(object):
     mesh    = None    #: A SimPEG.Mesh instance.
     mref    = None    #: Reference model.
 
-    def __init__(self, mesh=None, nP=None, mapping=None, indActive=None, **kwargs):
+    def __init__(self, mesh, 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':
@@ -314,18 +316,10 @@ 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.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 = mapping or self.mapPair(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
+        self.indActive = indActive
 
     @property
     def parent(self):
@@ -354,7 +348,7 @@ class BaseRegularization(object):
     @property
     def W(self):
         """Full regularization weighting matrix W."""
-        return sp.identity(self.nP)
+        return sp.identity(self.regmesh.nC)
 
     @Utils.timeIt
     def eval(self, m):
@@ -411,7 +405,238 @@ class BaseRegularization(object):
 
         return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
 
-class Tikhonov(BaseRegularization):
+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):
     """
     L2 Tikhonov regularization with both smallness and smoothness (first order
     derivative) contributions.
@@ -501,56 +726,131 @@ class Tikhonov(BaseRegularization):
             self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
         return self._Wzz
 
+
     @property
-    def Wsmooth(self):
+    def Wsmooth2(self):
         """Full smoothness regularization matrix W"""
         if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx, self.Wxx)
+            wlist = (self.Wxx)
             if self.regmesh.dim > 1:
-                wlist += (self.Wy, self.Wyy)
+                wlist += (self.Wyy)
             if self.regmesh.dim > 2:
-                wlist += (self.Wz, self.Wzz)
+                wlist += (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 _evalSmall(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
+    def _evalSmoothxx(self, m):
         if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
+            r = self.Wxx * ( self.mapping * (m - self.mref) )
         elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * (m) )
+            r = self.Wxx * ( self.mapping * (m) )
         return 0.5 * r.dot(r)
 
+    @Utils.timeIt
+    def _evalSmoothyy(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wyy * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wyy * ( 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)
-
-    @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) )
+        return self._evalSmall(m) + self._evalSmooth(m) + self._evalSmooth2(m)
 
     @Utils.timeIt
     def evalDeriv(self, m):
@@ -568,184 +868,134 @@ class Tikhonov(BaseRegularization):
             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
-        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+        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)
 
 
-class Simple(Tikhonov):
+
+class Sparse(Simple):
     """
-    Simple regularization that does not include length scales in the derivatives.
+        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.
     """
-
-    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.
+        
+    # 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
 
     def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
 
-        if isinstance(self.wght,float):
-            self.wght = np.ones(self.regmesh.nC) * self.wght
+        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.regmesh.vol*self.alpha_s*self.wght)**0.5)
+            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
+
         return self._Wsmall
 
     @property
     def Wx(self):
         """Regularization matrix Wx"""
-        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
+        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
+
         return self._Wx
 
     @property
     def Wy(self):
         """Regularization matrix Wy"""
-        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
+        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
+
         return self._Wy
 
     @property
     def Wz(self):
         """Regularization matrix Wz"""
-        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
+        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
+
         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.)
 

docs/api_Maps.rst

@@ -122,10 +122,92 @@ When these are used in the inverse problem, this is extremely important!!
 The API
 =======
 
-.. automodule:: SimPEG.Maps
+.. autoclass:: SimPEG.Maps.IdentityMap
     :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/examples/DC_Forward_PseudoSection.rst

@@ -22,7 +22,6 @@ 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/Mesh_Basic_ForwardDC.rst

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

tests/base/test_maps.py

@@ -5,10 +5,8 @@ from scipy.sparse.linalg import dsolve
 
 TOL = 1e-14
 
-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"]
-
+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"]
 
 class MapTests(unittest.TestCase):
 
@@ -19,8 +17,6 @@ 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:
@@ -32,15 +28,6 @@ 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.