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

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

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -257,7 +257,7 @@ class Fields3D_e(Fields):
         """
 
         # assuming primary does not depend on the model
-        return Zero()
+        return src.ePrimaryDeriv(self.prob, v, adjoint) #Zero()
 
     def _bPrimary(self, eSolution, srcList):
         """
@@ -600,8 +600,8 @@ class Fields3D_b(Fields):
 
 
         if adjoint:
-            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv
-        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv
+            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv + src.ePrimaryDeriv(self.prob, v, adjoint)
+        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv + src.ePrimaryDeriv(self.prob, v, adjoint)
 
     def _j(self, bSolution, srcList):
         """

SimPEG/EM/FDEM/ProblemFDEM.py

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

SimPEG/EM/FDEM/SrcFDEM.py

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

SimPEG/EM/Static/DC/RxDC.py

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

SimPEG/Maps.py

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

SimPEG/Regularization.py

@@ -39,7 +39,7 @@ class RegularizationMesh(object):
             if self.indActive is None:
                 self._nC = self.mesh.nC
             else:
-                self._nC = sum(self.indActive)
+                self._nC = int(sum(self.indActive))
         return self._nC
 
     @property
@@ -304,7 +304,7 @@ class BaseRegularization(object):
     mesh    = None    #: A SimPEG.Mesh instance.
     mref    = None    #: Reference model.
 
-    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+    def __init__(self, mesh=None, nP=None, mapping=None, indActive=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
         assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
         if indActive is not None and indActive.dtype != 'bool':
@@ -314,11 +314,19 @@ class BaseRegularization(object):
         if indActive is not None and mapping is None:
             mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
 
+        if mesh is None and nP is None:
+            raise Exception, 'either Mesh or number of parameters must be provided to the BaseRegularization'
+
         self.regmesh = RegularizationMesh(mesh,indActive)
-        self.mapping = mapping or self.mapPair(mesh)
-        self.mapping._assertMatchesPair(self.mapPair)
         self.indActive = indActive
 
+        if mesh is not None and nP is None:
+            nP = self.regmesh.nC
+        self.nP = nP
+
+        self.mapping = mapping or self.mapPair(nP=self.nP)
+        self.mapping._assertMatchesPair(self.mapPair)
+
     @property
     def parent(self):
         """This is the parent of the regularization."""
@@ -346,7 +354,7 @@ class BaseRegularization(object):
     @property
     def W(self):
         """Full regularization weighting matrix W."""
-        return sp.identity(self.regmesh.nC)
+        return sp.identity(self.nP)
 
     @Utils.timeIt
     def eval(self, m):

docs/api_Maps.rst

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

tests/base/test_maps.py

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