cleanup in mapping tests

Analytics

.bumpversion.cfg

@@ -1,4 +1,4 @@
 [bumpversion]
-current_version = 0.1.11
+current_version = 0.1.12
 files = setup.py SimPEG/__init__.py docs/conf.py
 

.gitignore

@@ -40,3 +40,4 @@ nosetests.xml
 docs/_build/
 Makefile
 docs/warnings.txt
+.DS_Store

.travis.yml

@@ -18,7 +18,7 @@ env:
     - TEST_DIR="tests/mesh tests/base tests/utils"
     - TEST_DIR=tests/em/fdem/inverse/derivs
     - TEST_DIR=tests/em/tdem
-    - TEST_DIR=tests/em/static
+    - TEST_DIR=tests/dcip
     - TEST_DIR=tests/flow
     - TEST_DIR=tests/mt
     - TEST_DIR=tests/examples

README.rst

@@ -1,4 +1,4 @@
-.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
+.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png
     :alt: SimPEG Logo
 
 ======

SimPEG/DCIP/BaseDC.py

@@ -0,0 +1,292 @@
+from SimPEG import *
+
+class FieldsDC_CC(Problem.Fields):
+    knownFields = {'phi_sol':'CC'}
+    aliasFields = {
+                    'phi' : ['phi_sol','CC','_phi'],
+                    'e' : ['phi_sol','F','_e'],
+                    'j' : ['phi_sol','F','_j']
+                  }
+
+    def __init__(self,mesh,survey,**kwargs):
+        super(FieldsDC_CC, self).__init__(mesh, survey, **kwargs)
+
+    def startup(self):
+        self._cellGrad = self.survey.prob.mesh.cellGrad
+        self._Mfinv = self.survey.prob.mesh.getFaceInnerProduct(invMat=True)
+
+    def _phi(self, phi_sol, srcList):
+        phi = phi_sol
+        # for i, src in enumerate(srcList):
+        #     phi_p = src.phi_p(self.survey.prob)
+        #     if phi_p is not None:
+        #         phi[:,i] += phi_p
+        return phi
+
+    def _e(self, phi_sol, srcList):
+        e = -self._cellGrad*phi_sol
+        # for i, src in enumerate(srcList):
+        #     e_p = src.e_p(self.survey.prob)
+        #     if e_p is not None:
+        #         e[:,i] += e_p
+        return e
+
+    def _j(self, phi_sol, srcList):
+
+        j = -self._Mfinv*self.survey.prob.Msig*self._cellGrad*phi_sol
+        # for i, src in enumerate(srcList):
+        #     j_p = src.j_p(self.survey.prob)
+        #     if j_p is not None:
+        #         j[:,i] += j_p
+        return j
+
+
+
+class SrcDipole(Survey.BaseSrc):
+    """A dipole source, locA and locB are moved to the closest cell-centers"""
+
+    current = 1
+    loc  = None
+    # _rhsDict = None
+
+    def __init__(self, rxList, locA, locB, **kwargs):
+        self.loc = (locA, locB)
+        super(SrcDipole, self).__init__(rxList, **kwargs)
+
+    def eval(self, prob):
+        # Recompute rhs
+        # if getattr(self, '_rhsDict', None) is None:
+        #     self._rhsDict = {}
+        # if mesh not in self._rhsDict:
+        pts = [self.loc[0], self.loc[1]]
+        inds = Utils.closestPoints(prob.mesh, pts)
+        q = np.zeros(prob.mesh.nC)
+        q[inds] = - self.current * ( np.r_[1., -1.] / prob.mesh.vol[inds] )
+        # self._rhsDict[mesh] = q
+        # return self._rhsDict[mesh]
+        return q
+
+
+class RxDipole(Survey.BaseRx):
+    """A dipole source, locA and locB are moved to the closest cell-centers"""
+    def __init__(self, locsM, locsN, **kwargs):
+        locs = (locsM, locsN)
+        assert locsM.shape == locsN.shape, 'locs must be the same shape.'
+        super(RxDipole, self).__init__(locs, 'dipole', storeProjections=False, **kwargs)
+
+    @property
+    def nD(self):
+        """Number of data in the receiver."""
+        return self.locs[0].shape[0]
+
+    def getP(self, mesh):
+        P0 = mesh.getInterpolationMat(self.locs[0], self.projGLoc)
+        P1 = mesh.getInterpolationMat(self.locs[1], self.projGLoc)
+        return P0 - P1
+
+
+class SurveyDC(Survey.BaseSurvey):
+    """
+        **SurveyDC**
+
+        Geophysical DC resistivity data.
+
+    """
+    uncert = None
+    def __init__(self, srcList, **kwargs):
+        self.srcList = srcList
+        Survey.BaseSurvey.__init__(self, **kwargs)
+        # self._rhsDict = {}
+        self._Ps = {}
+
+    def eval(self, u):
+        """
+            Predicted data.
+
+            .. math::
+                d_\\text{pred} = Pu(m)
+        """
+        P = self.getP(self.prob.mesh)
+        return P*mkvc(u[self.srcList, 'phi_sol'])
+
+    def getP(self, mesh):
+        if mesh in self._Ps:
+            return self._Ps[mesh]
+
+        P_src = [sp.vstack([rx.getP(mesh) for rx in src.rxList]) for src in self.srcList]
+
+        self._Ps[mesh] = sp.block_diag(P_src)
+        return self._Ps[mesh]
+
+
+class ProblemDC_CC(Problem.BaseProblem):
+    """
+        **ProblemDC**
+
+        Geophysical DC resistivity problem.
+
+    """
+
+    surveyPair = SurveyDC
+    Solver     = Solver
+    fieldsPair = FieldsDC_CC
+    Ainv = None
+
+    def __init__(self, mesh, **kwargs):
+        Problem.BaseProblem.__init__(self, mesh)
+        self.mesh.setCellGradBC('neumann')
+        Utils.setKwargs(self, **kwargs)
+
+
+    deleteTheseOnModelUpdate = ['_A', '_Msig', '_dMdsig']
+
+    @property
+    def Msig(self):
+        if getattr(self, '_Msig', None) is None:
+            sigma = self.curModel.transform
+            Av = self.mesh.aveF2CC
+            self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
+        return self._Msig
+
+    @property
+    def dMdsig(self):
+        if getattr(self, '_dMdsig', None) is None:
+            sigma = self.curModel.transform
+            Av = self.mesh.aveF2CC
+            dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2)
+            self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop
+        return self._dMdsig
+
+    @property
+    def A(self):
+        """
+            Makes the matrix A(m) for the DC resistivity problem.
+
+            :param numpy.ndarray m: model
+            :rtype: scipy.sparse.csc_matrix
+            :return: A(m)
+
+            .. math::
+                c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
+
+            Where M() is the mass matrix and mT is the model transform.
+        """
+        if getattr(self, '_A', None) is None:
+            D = self.mesh.faceDiv
+            G = self.mesh.cellGrad
+            self._A = D*self.Msig*G
+            # Remove the null space from the matrix.
+            self._A[0,0] /= self.mesh.vol[0]
+            self._A = self._A.tocsc()
+        return self._A
+
+    def getRHS(self):
+        # if self.mesh not in self._rhsDict:
+        RHS = np.array([src.eval(self) for src in self.survey.srcList]).T
+        # self._rhsDict[mesh] = RHS
+        # return self._rhsDict[mesh]
+        return RHS
+
+    def fields(self, m):
+
+        F = self.fieldsPair(self.mesh, self.survey)
+        self.curModel = m
+        A    = self.A
+        self.Ainv = self.Solver(A, **self.solverOpts)
+        RHS    = self.getRHS()
+        Phi  = self.Ainv * RHS
+        Srcs = self.survey.srcList
+        F[Srcs, 'phi_sol'] = Phi
+
+        return F
+
+    def Jvec(self, m, v, f=None):
+        """
+            :param numpy.array m: model
+            :param numpy.array v: vector to multiply
+            :param Fields f: fields
+            :rtype: numpy.array
+            :return: Jv
+
+            .. math::
+                c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
+
+                \\nabla_u (A(m)u - q) = A(m)
+
+                \\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m(M(mT(m)))
+
+            Where M() is the mass matrix and mT is the model transform.
+
+            .. math::
+                J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u)
+
+                J(v) = - P ( A(m)^{-1} ( G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) v ) )
+        """
+        # Set current model; clear dependent property $\mathbf{A(m)}$
+        self.curModel = m
+        sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
+        if f is None:
+            # Run forward simulation if $u$ not provided
+            f = self.fields(self.curModel)
+        u = f[self.survey.srcList, 'phi_sol']
+
+        D = self.mesh.faceDiv
+        G = self.mesh.cellGrad
+        # Derivative of model transform, $\deriv{\sigma}{\m}$
+        dsigdm_x_v = self.curModel.transformDeriv * v
+
+        # Take derivative of $C(m,u)$ w.r.t. $m$
+        dCdm_x_v = np.empty_like(u)
+        # loop over fields for each source
+        for i in range(self.survey.nSrc):
+            # Derivative of inner product, $\left(\mathbf{M}_{1/\sigma}^f\right)^{-1}$
+            dAdsig         = D * self.dMdsig( G * u[:,i] )
+            dCdm_x_v[:, i] = dAdsig *  dsigdm_x_v
+
+        # Take derivative of $C(m,u)$ w.r.t. $u$
+        dA_du = self.A
+        # Solve for $\deriv{u}{m}$
+        # dCdu_inv = self.Solver(dCdu, **self.solverOpts)
+        if self.Ainv is None:
+            self.Ainv = self.Solver(dA_du, **self.solverOpts)
+
+        P     = self.survey.getP(self.mesh)
+        Jv    = - P * mkvc( self.Ainv * dCdm_x_v )
+        return Jv
+
+    def Jtvec(self, m, v, f=None):
+
+        self.curModel = m
+        sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
+        if f is None:
+            # Run forward simulation if $f$ not provided
+            f = self.fields(self.curModel)
+        u = f[self.survey.srcList, 'phi_sol']
+
+        shp = u.shape
+        P = self.survey.getP(self.mesh)
+        PT_x_v = (P.T*v).reshape(shp, order='F')
+
+        D = self.mesh.faceDiv
+        G = self.mesh.cellGrad
+        dA_du = self.A
+        mT_dm = self.mapping.deriv(m)
+
+        # We probably always need this due to the linesearch .. (?)
+        self.Ainv = self.Solver(dA_du.T, **self.solverOpts)
+        # if self.Ainv is None:
+        #     self.Ainv = self.Solver(dCdu, **self.solverOpts)
+
+        w = self.Ainv * PT_x_v
+
+        Jtv = 0
+        for i, ui in enumerate(u.T):  # loop over each column
+            Jtv += self.dMdsig( G * ui ).T * ( D.T * w[:,i] )
+
+        Jtv = - mT_dm.T * ( Jtv )
+        return Jtv
+
+
+
+
+

SimPEG/DCIP/BaseIP.py

@@ -0,0 +1,182 @@
+from SimPEG import *
+from BaseDC import SurveyDC, FieldsDC_CC
+
+class SurveyIP(SurveyDC):
+    """
+        **SurveyDC**
+
+        Geophysical DC resistivity data.
+
+    """
+
+    def __init__(self, srcList, **kwargs):
+        self.srcList = srcList
+        Survey.BaseSurvey.__init__(self, **kwargs)
+        self._Ps = {}
+
+    def dpred(self, m, f=None):
+        """
+            Predicted data.
+
+            .. math::
+                d_\\text{pred} = Pf(m)
+        """
+
+        return self.prob.forward(m)
+
+
+class ProblemIP(Problem.BaseProblem):
+    """
+        **ProblemIP**
+
+        Geophysical IP resistivity problem.
+
+    """
+
+    surveyPair = SurveyDC
+    Solver     = Solver
+    sigma = None
+    Ainv = None
+    u = None
+
+    def __init__(self, mesh, **kwargs):
+        Problem.BaseProblem.__init__(self, mesh)
+        self.mesh.setCellGradBC('neumann')
+        Utils.setKwargs(self, **kwargs)
+
+    # deleteTheseOnModelUpdate = ['_A', '_Msig', '_dMdsig']
+
+    @property
+    def Msig(self):
+        if getattr(self, '_Msig', None) is None:
+            # sigma = self.curModel.transform
+            sigma = self.sigma
+            Av = self.mesh.aveF2CC
+            self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
+        return self._Msig
+
+    @property
+    def dMdsig(self):
+        if getattr(self, '_dMdsig', None) is None:
+            # sigma = self.curModel.transform
+            sigma = self.sigma
+            Av = self.mesh.aveF2CC
+            dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2)
+            self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop
+        return self._dMdsig
+
+    @property
+    def A(self):
+        """
+            Makes the matrix A(m) for the DC resistivity problem.
+
+            :param numpy.array m: model
+            :rtype: scipy.sparse.csc_matrix
+            :return: A(m)
+
+            .. math::
+                c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
+
+            Where M() is the mass matrix and mT is the model transform.
+        """
+        if getattr(self, '_A', None) is None:
+            D = self.mesh.faceDiv
+            G = self.mesh.cellGrad
+            self._A = D*self.Msig*G
+            # Remove the null space from the matrix.
+            self._A[-1,-1] /= self.mesh.vol[-1]
+            self._A = self._A.tocsc()
+        return self._A
+
+    def getRHS(self):
+        # if self.mesh not in self._rhsDict:
+        RHS = np.array([src.eval(self) for src in self.survey.srcList]).T
+        # self._rhsDict[mesh] = RHS
+        # return self._rhsDict[mesh]
+        return RHS
+
+    def fields(self, m):
+        if self.u is None:
+	        A    = self.A
+	        if self.Ainv == None:
+	        	self.Ainv = self.Solver(A, **self.solverOpts)
+	        Q    = self.getRHS()
+	        self.u  = self.Ainv * Q
+        return self.u
+
+    def forward(self, m, u=None):
+        # Set current model; clear dependent property $\mathbf{A(m)}$
+        self.curModel = m
+        # sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
+        sigma = self.sigma
+        if self.u is None:
+            # Run forward simulation if $u$ not provided
+            u = self.fields(sigma)
+
+        shp = (self.mesh.nC, self.survey.nSrc)
+        u = self.u.reshape(shp, order='F')
+
+        D = self.mesh.faceDiv
+        G = self.mesh.cellGrad
+        # Derivative of model transform, $\deriv{\sigma}{\m}$
+        # dsigdm_x_v = self.curModel.transformDeriv * v
+
+        dsigdm_x_v = Utils.sdiag(sigma) * self.curModel.transformDeriv *  m
+
+        # Take derivative of $C(m,u)$ w.r.t. $m$
+        dCdm_x_v = np.empty_like(u)
+        # loop over fields for each source
+        for i in range(self.survey.nSrc):
+            # Derivative of inner product, $\left(\mathbf{M}_{1/\sigma}^f\right)^{-1}$
+            dAdsig         = D * self.dMdsig( G * u[:,i] )
+            dCdm_x_v[:, i] = dAdsig *  dsigdm_x_v
+
+        # Take derivative of $C(m,u)$ w.r.t. $u$
+
+        if self.Ainv == None:
+        	self.Ainv = self.Solver(A, **self.solverOpts)
+
+        # dCdu = self.A
+        # Solve for $\deriv{u}{m}$
+        # dCdu_inv = self.Solver(dCdu, **self.solverOpts)
+        P        = self.survey.getP(self.mesh)
+        J_x_v    = - P * mkvc( self.Ainv * dCdm_x_v )
+        return -J_x_v
+
+    def Jvec(self, m, v, f=None):
+        return self.forward(v)
+
+    def Jtvec(self, m, v, f=None):
+
+        self.curModel = m
+        # sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
+        sigma = self.sigma
+        if self.u is None:
+            u = self.fields(sigma)
+        else:
+        	u = self.u
+        shp = (self.mesh.nC, self.survey.nSrc)
+        u = u.reshape(shp, order='F')
+        P = self.survey.getP(self.mesh)
+        PT_x_v = (P.T*v).reshape(shp, order='F')
+
+        D = self.mesh.faceDiv
+        G = self.mesh.cellGrad
+        A = self.A
+        mT_dm = Utils.sdiag(sigma)*self.mapping.deriv(m)
+        # mT_dm = self.mapping.deriv(m)
+
+        # dCdu = A.T
+        # Ainv = self.Solver(dCdu, **self.solverOpts)
+        # if self.Ainv == None:
+    	self.Ainv = self.Solver(A.T, **self.solverOpts)
+
+        w = self.Ainv * PT_x_v
+
+        Jtv = 0
+        for i, ui in enumerate(u.T):  # loop over each column
+            Jtv += self.dMdsig( G * ui ).T * ( D.T * w[:,i] )
+
+        Jtv = - mT_dm.T * ( Jtv )
+        return -Jtv
+

SimPEG/DCIP/Utils.py

@@ -0,0 +1,38 @@
+import numpy as np
+
+def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
+
+    import SimPEG.DCIP as DC
+
+    elocs = np.arange(0,aSpacing*nElecs,aSpacing)
+    elocs -= (nElecs*aSpacing - aSpacing)/2
+    space = 1
+    WENNER = np.zeros((0,),dtype=int)
+    for ii in range(nElecs):
+        for jj in range(nElecs):
+            test = np.r_[jj,jj+space,jj+space*2,jj+space*3]
+            if np.any(test >= nElecs):
+                break
+            WENNER = np.r_[WENNER, test]
+        space += 1
+    WENNER = WENNER.reshape((-1,4))
+
+
+    if plotIt:
+        for i, s in enumerate('rbkg'):
+            plt.plot(elocs[WENNER[:,i]],s+'.')
+        plt.show()
+
+    # Create sources and receivers
+    i = 0
+    if in2D:
+        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0]
+    else:
+        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0, 0]
+    srcList = []
+    for i in range(WENNER.shape[0]):
+        rx = DC.RxDipole(getLoc(i,1),getLoc(i,2))
+        src = DC.SrcDipole([rx], getLoc(i,0),getLoc(i,3))
+        srcList += [src]
+
+    return srcList

SimPEG/DCIP/__init__.py

@@ -0,0 +1,4 @@
+from BaseDC import *
+from BaseIP import *
+from DCIPUtils import *
+import Utils

SimPEG/Directives.py

@@ -253,8 +253,7 @@ class SaveOutputDictEveryIteration(SaveEveryIteration):
 class Update_IRLS(InversionDirective):
 
     eps_min = None
-    eps_p = None
-    eps_q = None
+    eps = None
     norms = [2.,2.,2.,2.]
     factor = None
     gamma = None
@@ -263,6 +262,7 @@ class Update_IRLS(InversionDirective):
     f_old = None
     f_min_change = 1e-2
     beta_tol = 5e-2
+    prctile = 95
 
     # Solving parameter for IRLS (mode:2)
     IRLSiter   = 0
@@ -297,9 +297,22 @@ class Update_IRLS(InversionDirective):
             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
+
+            # Either use the supplied epsilon, or fix base on distribution of
+            # model values
+            if getattr(self, 'reg.eps', None) is None:
+                self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile)
+            else:                 
+                self.reg.eps_p = self.eps[0]
+                
+            if getattr(self, 'reg.eps', None) is None:
+                self.reg.eps_q = np.percentile(np.abs(self.reg.regmesh.cellDiffxStencil*(self.reg.mapping * self.invProb.curModel)),self.prctile)
+            else:                 
+                self.reg.eps_q = self.eps[1]
+            
+            print "L[p qx qy qz]-norm : " + str(self.reg.norms)
+            print "eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q)
+            
             self.reg.norms = self.norms
             self.coolingFactor = 1.
             self.coolingRate = 1
@@ -343,14 +356,14 @@ class Update_IRLS(InversionDirective):
             else:
                 self.f_old = phim_new
 
-            # 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])
-                else:
-                    self.reg.eps = eps
+#            # 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])
+#                else:
+#                    self.reg.eps = eps
 
             # Get phi_m at the end of current iteration
             self.phi_m_last = self.invProb.phi_m_last

SimPEG/EM/Analytics/FDEMDipolarfields.py

@@ -114,7 +114,7 @@ def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l
     """
     mu = mu_0*(1+kappa)
     epsilon = epsilon_0*epsr
-    sig_hat = sig + 1j*omeg*epsilon
+    sig_hat = sig + 1j*omega(f)*epsilon
 
     XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
     # Check
@@ -160,7 +160,7 @@ def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1.,
             Add description of parameters
     """
 
-    Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
+    Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
     Jx = sig*Ex
     Jy = sig*Ey
     Jz = sig*Ez
@@ -175,7 +175,7 @@ def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., le
             Add description of parameters
     """
 
-    Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
+    Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
     Jx_galvanic = sig*Ex_galvanic
     Jy_galvanic = sig*Ey_galvanic
     Jz_galvanic = sig*Ez_galvanic
@@ -190,7 +190,7 @@ def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., l
             Add description of parameters
     """
 
-    Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
+    Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
     Jx_inductive = sig*Ex_inductive
     Jy_inductive = sig*Ey_inductive
     Jz_inductive = sig*Ez_inductive
@@ -248,7 +248,7 @@ def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1.,
             Add description of parameters
     """
 
-    Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
+    Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
     Bx = mu*Hx
     By = mu*Hy
     Bz = mu*Hz

SimPEG/EM/Static/Utils/StaticUtils.py

@@ -219,12 +219,10 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
         for ii in range(0, int(nstn)-1):
 
 
-            if stype == 'dipole-dipole':
+            if stype == 'dpdp':
                 tx = np.c_[M[ii,:],N[ii,:]]
-            elif stype == 'pole-dipole':
+            elif stype == 'pdp':
                 tx = np.c_[M[ii,:],M[ii,:]]
-            else:
-                raise Exception('The stype must be "dipole-dipole" or "pole-dipole"')
 
             # Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]])
 
@@ -258,10 +256,10 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
                 P2 = np.c_[stn_x+a*dl_x, np.ones(nstn).T*ztop]
                 rxClass = DC.Rx.Dipole_ky(P1, P2)
 
-            if stype == 'dipole-dipole':
-                srcClass = DC.Src.Dipole([rxClass], M[ii,:],N[ii,:])
-            elif stype == 'pole-dipole':
-                srcClass = DC.Src.Pole([rxClass], M[ii,:])
+            if stype == 'dpdp':
+                    srcClass = DC.Src.Dipole([rxClass], M[ii,:],N[ii,:])
+            elif stype == 'pdp':
+                    srcClass = DC.Src.Pole([rxClass], M[ii,:])
             SrcList.append(srcClass)
 
     elif stype == 'gradient':
@@ -312,7 +310,7 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
             srcClass = DC.Src.Dipole([rxClass], M[0,:], N[-1,:])
         SrcList.append(srcClass)
     else:
-        print """stype must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """
+        print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
 
 
     return SrcList

SimPEG/EM/Static/__init__.py

@@ -1,4 +1,3 @@
 import DC
 import IP
 import SIP
-import Utils

SimPEG/Examples/DC_Analytic_Dipole.py

@@ -1,7 +1,7 @@
 from SimPEG import *
 import SimPEG.EM.Static.DC as DC
 
-def run(plotIt=False):
+def run(plotIt=True):
     cs = 25.
     hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
     hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
@@ -65,4 +65,4 @@ def run(plotIt=False):
 
 
 if __name__ == '__main__':
-    print run(plotIt=True)
+    print run()

SimPEG/Examples/DC_Forward_PseudoSection.py

@@ -1,5 +1,5 @@
 from SimPEG import Mesh, Utils, np, sp
-from SimPEG.EM.Static import DC, Utils as StaticUtils
+import SimPEG.DCIP as DC
 import time
 
 def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', unitType='appConductivity', plotIt=True):
@@ -74,13 +74,7 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
 
     # We will handle the geometry of the survey for you and create all the combination of tx-rx along line
     # [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])
-    srcList = StaticUtils.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])
-    print(srcList)
-
-    raise Exception('The output of the function changed. We need to create the survey and get Rx to continue.')
-
-    # Something like this!
-    survey, Tx, Rx
+    survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])
 
     # Define some global geometry
     dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )

SimPEG/Examples/Inversion_IRLS.py

@@ -42,55 +42,33 @@ def run(N=100, plotIt=True):
     survey = Survey.LinearSurvey()
     survey.pair(prob)
     survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk)
-    #survey.makeSyntheticData(mtrue, std=std_noise)
 
     wd = np.ones(nk) * std_noise
 
-    #print survey.std[0]
-    #M = prob.mesh
     # Distance weighting
     wr = np.sum(prob.G**2.,axis=0)**0.5
     wr = ( wr/np.max(wr) )
 
-#    reg = Regularization.Simple(mesh)
-#    reg.mref = mref
-#    reg.cell_weights = wr
-#
     dmis = DataMisfit.l2_DataMisfit(survey)
     dmis.Wd = 1./wd
-#
-#    opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4)
-#    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
-#    invProb.curModel = m0
-#
-#    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
-#    target = Directives.TargetMisfit()
-#
+
     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
 
     reg = Regularization.Sparse(mesh)
     reg.mref = mref
     reg.cell_weights = wr
 
     reg.mref = np.zeros(mesh.nC)
-    eps_p = 5e-2
-    eps_q = 5e-2
-    norms   = [0., 0., 2., 2.]
+    
 
     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)
+    
+    # Set the IRLS directive, penalize the lowest 25 percentile of model values
+    # Start with an l2-l2, then switch to lp-norms
+    norms   = [0., 0., 2., 2.]    
+    IRLS = Directives.Update_IRLS( norms=norms, prctile = 25, maxIRLSiter = 15, minGNiter=3)
 
     inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
 

SimPEG/Examples/Maps_ComboMaps.py

@@ -0,0 +1,62 @@
+from SimPEG import Mesh, Maps, np
+
+def run(plotIt=True):
+    """
+
+        Maps: ComboMaps
+        ===============
+
+        We will use an example where we want a 1D layered earth as
+        our model, but we want to map this to a 2D discretization to do our forward
+        modeling. We will also assume that we are working in log conductivity still,
+        so after the transformation we want to map to conductivity space.
+        To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`),
+        which does the first part of what we just described. The second part will be
+        done by the :class:`SimPEG.Maps.ExpMap` described above.
+
+        .. code-block:: python
+            :linenos:
+
+            M = Mesh.TensorMesh([7,5])
+            v1dMap = Maps.SurjectVertical1D(M)
+            expMap = Maps.ExpMap(M)
+            myMap = expMap * v1dMap
+            m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
+            sig = myMap * m
+
+        If you noticed, it was pretty easy to combine maps. What is even cooler is
+        that the derivatives also are made for you (if everything goes right).
+        Just to be sure that the derivative is correct, you should always run the test
+        on the mapping that you create.
+
+    """
+
+
+    M = Mesh.TensorMesh([7,5])
+    v1dMap = Maps.SurjectVertical1D(M)
+    expMap = Maps.ExpMap(M)
+    myMap = expMap * v1dMap
+    m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
+    sig = myMap * m
+
+    if not plotIt: return
+
+    import matplotlib.pyplot as plt
+    figs, axs = plt.subplots(1,2)
+    axs[0].plot(m, M.vectorCCy, 'b-o')
+    axs[0].set_title('Model')
+    axs[0].set_ylabel('Depth, y')
+    axs[0].set_xlabel('Value, $m_i$')
+    axs[0].set_xlim(0,3)
+    axs[0].set_ylim(0,1)
+    clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0])
+    axs[1].set_title('Physical Property')
+    axs[1].set_ylabel('Depth, y')
+    clbar.set_label('$\sigma = \exp(\mathbf{P}m)$')
+    plt.tight_layout()
+    plt.show()
+
+
+if __name__ == '__main__':
+    run()
+

SimPEG/Examples/Maps_Mesh2Mesh.py

@@ -0,0 +1,41 @@
+from SimPEG import Mesh, Maps, Utils
+
+def run(plotIt=True):
+    """
+
+        Maps: Mesh2Mesh
+        ===============
+
+        This mapping allows you to go from one mesh to another.
+
+    """
+
+    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
+
+    if not plotIt: return
+
+    import matplotlib.pyplot as plt
+    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()
+
+
+if __name__ == '__main__':
+    run()
+

SimPEG/Examples/Utils_surface2ind_topo.py

@@ -2,7 +2,7 @@ from SimPEG import *
 from SimPEG.Utils import surface2ind_topo
 
 
-def run(plotIt=False, nx=5, ny=5):
+def run(plotIt=True, nx=5, ny=5):
     """
 
         Utils: surface2ind_topo

SimPEG/Examples/__init__.py

@@ -10,6 +10,8 @@ import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
 import Inversion_IRLS
 import Inversion_Linear
+import Maps_ComboMaps
+import Maps_Mesh2Mesh
 import Mesh_Basic_ForwardDC
 import Mesh_Basic_PlotImage
 import Mesh_Basic_Types
@@ -22,7 +24,7 @@ import MT_1D_ForwardAndInversion
 import MT_3D_Foward
 import Utils_surface2ind_topo
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Maps_ComboMaps", "Maps_Mesh2Mesh", "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,6 @@
-import Utils, numpy as np, scipy.sparse as sp
+import Utils
+import numpy as np
+import scipy.sparse as sp
 from scipy.sparse.linalg import LinearOperator
 from Tests import checkDerivative
 from PropMaps import PropMap, Property
@@ -6,6 +8,7 @@ from numpy.polynomial import polynomial
 from scipy.interpolate import UnivariateSpline
 import warnings
 
+
 class IdentityMap(object):
     """
     SimPEG Map
@@ -17,10 +20,11 @@ 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], 'Number of parameters '
+            'must be an integer.'
 
         self.mesh = mesh
-        self._nP  = nP
+        self._nP = nP
 
     @property
     def nP(self):
@@ -50,14 +54,14 @@ class IdentityMap(object):
             return ('*', self.nP)
         return (self.mesh.nC, self.nP)
 
-
     def _transform(self, m):
         """
             Changes the model into the physical property.
 
             .. note::
 
-                This can be called by the __mul__ property against a numpy.ndarray.
+                This can be called by the __mul__ property against a
+                :meth:numpy.ndarray.
 
             :param numpy.array m: model
             :rtype: numpy.array
@@ -81,7 +85,7 @@ class IdentityMap(object):
         """
         raise NotImplementedError('The transformInverse is not implemented.')
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
         """
             The derivative of the transformation.
 
@@ -90,13 +94,16 @@ class IdentityMap(object):
             :return: derivative of transformed model
 
         """
+        if v is not None:
+            return v
         return sp.identity(self.nP)
 
     def test(self, m=None, **kwargs):
         """Test the derivative of the mapping.
 
             :param numpy.array m: model
-            :param kwargs: key word arguments of :meth:`SimPEG.Tests.checkDerivative`
+            :param kwargs: key word arguments of
+                           :meth:`SimPEG.Tests.checkDerivative`
             :rtype: bool
             :return: passed the test?
 
@@ -106,26 +113,52 @@ class IdentityMap(object):
             m = abs(np.random.rand(self.nP))
         if 'plotIt' not in kwargs:
             kwargs['plotIt'] = False
-        return checkDerivative(lambda m : [self * m, self.deriv(m)], m, num=4, **kwargs)
+        return checkDerivative(lambda m : [self * m, self.deriv(m)], m, num=4,
+                               **kwargs)
+
+    def testVec(self, m=None, **kwargs):
+        """Test the derivative of the mapping times a vector.
+
+            :param numpy.array m: model
+            :param kwargs: key word arguments of
+                           :meth:`SimPEG.Tests.checkDerivative`
+            :rtype: bool
+            :return: passed the test?
+
+        """
+        print 'Testing %s' % str(self)
+        if m is None:
+            m = abs(np.random.rand(self.nP))
+        if 'plotIt' not in kwargs:
+            kwargs['plotIt'] = False
+        return checkDerivative(lambda m: [self*m, lambda x: self.deriv(m, x)],
+                               m, num=4, **kwargs)
 
     def _assertMatchesPair(self, pair):
         assert (isinstance(self, pair) or
-            isinstance(self, ComboMap) and isinstance(self.maps[0], pair)
-            ), "Mapping object must be an instance of a %s class."%(pair.__name__)
+                isinstance(self, ComboMap) and isinstance(self.maps[0], pair)
+                ), ("Mapping object must be an instance of a %s"
+                    " class."% (pair.__name__))
 
     def __mul__(self, val):
         if isinstance(val, IdentityMap):
-            if not (self.shape[1] == '*' or val.shape[0] == '*') and not self.shape[1] == val.shape[0]:
-                raise ValueError('Dimension mismatch in %s and %s.' % (str(self), str(val)))
+            if (not (self.shape[1] == '*' or val.shape[0] == '*') and not
+                self.shape[1] == val.shape[0]):
+                raise ValueError('Dimension mismatch in %s and %s.'
+                                 % (str(self), str(val)))
             return ComboMap([self, val])
+
         elif isinstance(val, np.ndarray):
             if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]:
-                raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape)))
+                raise ValueError('Dimension mismatch in %s and np.ndarray%s.'
+                                 % (str(self), str(val.shape)))
             return self._transform(val)
-        raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!')
+        raise Exception('Unrecognized data type to multiply. '
+                        'Try a map or a numpy.ndarray!')
 
     def __str__(self):
-        return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
+        return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0],
+                              self.shape[1])
 
 
 class ComboMap(IdentityMap):
@@ -136,11 +169,17 @@ class ComboMap(IdentityMap):
 
         self.maps = []
         for ii, m in enumerate(maps):
-            assert isinstance(m, IdentityMap), 'Unrecognized data type, inherit from an IdentityMap or ComboMap!'
-            if ii > 0 and not (self.shape[1] == '*' or m.shape[0] == '*') and not self.shape[1] == m.shape[0]:
+            assert isinstance(m, IdentityMap), "Unrecognized data type, "
+            "inherit from an IdentityMap or ComboMap!"
+
+            if (ii > 0 and not (self.shape[1] == '*' or m.shape[0] == '*') and
+                not self.shape[1] == m.shape[0]):
                 prev = self.maps[-1]
-                errArgs = (prev.__class__.__name__, prev.shape[0], prev.shape[1], m.__class__.__name__, m.shape[0], m.shape[1])
-                raise ValueError('Dimension mismatch in map[%s] (%s, %s) and map[%s] (%s, %s).' % errArgs)
+                errArgs = (prev.__class__.__name__, prev.shape[0],
+                           prev.shape[1], m.__class__.__name__, m.shape[0],
+                           m.shape[1])
+                raise ValueError('Dimension mismatch in map[%s] (%s, %s) '
+                                 'and map[%s] (%s, %s).' % errArgs)
 
             if isinstance(m, ComboMap):
                 self.maps += m.maps
@@ -164,8 +203,13 @@ class ComboMap(IdentityMap):
             m = map_i * m
         return m
 
-    def deriv(self, m):
-        deriv = 1
+    def deriv(self, m, v=None):
+
+        if v is not None:
+            deriv = v
+        else:
+            deriv = 1
+
         mi = m
         for map_i in reversed(self.maps):
             deriv = map_i.deriv(mi) * deriv
@@ -212,8 +256,7 @@ class ExpMap(IdentityMap):
         """
         return np.log(Utils.mkvc(D))
 
-
-    def deriv(self, m):
+    def deriv(self, m, v=None):
         """
             :param numpy.array m: model
             :rtype: scipy.sparse.csr_matrix
@@ -236,7 +279,11 @@ class ExpMap(IdentityMap):
 
                 \\frac{\partial \exp{m}}{\partial m} = \\text{sdiag}(\exp{m})
         """
-        return Utils.sdiag(np.exp(Utils.mkvc(m)))
+        deriv = Utils.sdiag(np.exp(Utils.mkvc(m)))
+        if v is not None:
+            return deriv * v
+        return deriv
+
 
 class ReciprocalMap(IdentityMap):
     """
@@ -253,10 +300,12 @@ class ReciprocalMap(IdentityMap):
     def inverse(self, D):
         return 1.0 / Utils.mkvc(m)
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
         # TODO: if this is a tensor, you might have a problem.
-        return Utils.sdiag( - Utils.mkvc(m)**(-2) )
-
+        deriv = Utils.sdiag( - Utils.mkvc(m)**(-2) )
+        if v is not None:
+            return deriv * v
+        return deriv
 
 
 class LogMap(IdentityMap):
@@ -265,13 +314,13 @@ class LogMap(IdentityMap):
 
         If \\(p\\) is the physical property and \\(m\\) is the model, then
 
-        ..math::
+        .. math::
 
             p = \\log(m)
 
         and
 
-        ..math::
+        .. math::
 
             m = \\exp(p)
 
@@ -286,17 +335,20 @@ class LogMap(IdentityMap):
     def _transform(self, m):
         return np.log(Utils.mkvc(m))
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
         mod = Utils.mkvc(m)
         deriv = np.zeros(mod.shape)
         tol = 1e-16 # zero
-        ind = np.greater_equal(np.abs(mod),tol)
+        ind = np.greater_equal(np.abs(mod), tol)
         deriv[ind] = 1.0/mod[ind]
+        if v is not None:
+            return Utils.sdiag(deriv)*v
         return Utils.sdiag(deriv)
 
     def inverse(self, m):
         return np.exp(Utils.mkvc(m))
 
+
 class SurjectFull(IdentityMap):
     """
     SurjectFull
@@ -305,8 +357,8 @@ class SurjectFull(IdentityMap):
     full model space.
     """
 
-    def __init__(self,mesh,**kwargs):
-        IdentityMap.__init__(self, mesh,**kwargs)
+    def __init__(self, mesh, **kwargs):
+        IdentityMap.__init__(self, mesh, **kwargs)
 
     @property
     def nP(self):
@@ -318,22 +370,28 @@ class SurjectFull(IdentityMap):
             :rtype: numpy.array
             :return: transformed model
         """
-        return np.ones(self.mesh.nC)*m
+        return np.ones(self.mesh.nC) * m
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
         """
             :param numpy.array m: model
             :rtype: numpy.array
             :return: derivative of transformed model
         """
-        return np.ones([self.mesh.nC,1])
+        deriv = np.ones([self.mesh.nC,1])
+        if v is not None:
+            return deriv * v
+        return deriv
 
 class FullMap(SurjectFull):
-    def __init__(self,mesh,**kwargs):
+    """FullMap is depreciated. Use SurjectVertical1DMap instead.
+    """
+    def __init__(self, mesh, **kwargs):
         warnings.warn(
             "`FullMap` is deprecated and will be removed in future versions. Use `SurjectFull` instead",
             FutureWarning)
-        SurjectFull.__init__(self,mesh,**kwargs)
+        SurjectFull.__init__(self, mesh, **kwargs)
+
 
 class SurjectVertical1D(IdentityMap):
     """SurjectVertical1DMap
@@ -363,7 +421,7 @@ class SurjectVertical1D(IdentityMap):
         repNum = self.mesh.vnC[:self.mesh.dim-1].prod()
         return Utils.mkvc(m).repeat(repNum)
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
         """
             :param numpy.array m: model
             :rtype: scipy.sparse.csr_matrix
@@ -374,14 +432,22 @@ class SurjectVertical1D(IdentityMap):
                     (np.ones(repNum),
                     (range(repNum), np.zeros(repNum))
                     ), shape=(repNum, 1))
-        return sp.kron(sp.identity(self.nP), repVec)
+        deriv = sp.kron(sp.identity(self.nP), repVec)
+        if v is not None:
+            return deriv * v
+        return deriv
+
 
 class Vertical1DMap(SurjectVertical1D):
-    def __init__(self,mesh,**kwargs):
+    """
+        Vertical1DMap is depreciated. Use SurjectVertical1D instead.
+    """
+    def __init__(self, mesh, **kwargs):
         warnings.warn(
             "`Vertical1DMap` is deprecated and will be removed in future versions. Use `SurjectVertical1D` instead",
             FutureWarning)
-        SurjectVertical1D.__init__(self,mesh,**kwargs)
+        SurjectVertical1D.__init__(self, mesh, **kwargs)
+
 
 class Surject2Dto3D(IdentityMap):
     """Map2Dto3D
@@ -390,12 +456,12 @@ class Surject2Dto3D(IdentityMap):
         3D model space.
     """
 
-    normal = 'Y' #: The normal
+    normal = 'Y'  #: The normal
 
     def __init__(self, mesh, **kwargs):
         assert mesh.dim == 3, 'Only works for a 3D Mesh'
         IdentityMap.__init__(self, mesh, **kwargs)
-        assert self.normal in ['X','Y','Z'], 'For now, only "Y" normal is supported'
+        assert self.normal in ['X', 'Y', 'Z'], 'For now, only "Y" normal is supported'
 
     @property
     def nP(self):
@@ -424,7 +490,7 @@ class Surject2Dto3D(IdentityMap):
         elif self.normal == 'X':
             return Utils.mkvc(m.reshape(self.mesh.vnC[[1,2]], order='F')[np.newaxis,:,:].repeat(self.mesh.nCx,axis=0))
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
         """
             :param numpy.array m: model
             :rtype: scipy.sparse.csr_matrix
@@ -436,19 +502,25 @@ class Surject2Dto3D(IdentityMap):
                     (np.ones(nC),
                     (range(nC), inds)
                 ), shape=(nC, nP))
+        if v is not None:
+            return P * v
         return P
 
+
 class Map2Dto3D(Surject2Dto3D):
-    def __init__(self,mesh,**kwargs):
+    """Map2Dto3D is depreciated. Use Surject2Dto3D instead
+    """
+
+    def __init__(self, mesh, **kwargs):
         warnings.warn(
             "`Map2Dto3D` is deprecated and will be removed in future versions. Use `Surject2Dto3D` instead",
             FutureWarning)
-        Surject2Dto3D.__init__(self,mesh,**kwargs)
+        Surject2Dto3D.__init__(self, mesh, **kwargs)
+
 
 class Mesh2Mesh(IdentityMap):
     """
         Takes a model on one mesh are translates it to another mesh.
-
     """
 
     def __init__(self, meshes, **kwargs):
@@ -461,7 +533,7 @@ class Mesh2Mesh(IdentityMap):
         self.mesh  = meshes[0]
         self.mesh2 = meshes[1]
 
-        self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC,'CC',zerosOutside=True)
+        self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC, 'CC', zerosOutside=True)
 
     @property
     def shape(self):
@@ -472,9 +544,13 @@ class Mesh2Mesh(IdentityMap):
     def nP(self):
         """Number of parameters in the model."""
         return self.mesh2.nC
+
     def _transform(self, m):
-        return self.P*m
-    def deriv(self, m):
+        return self.P * m
+
+    def deriv(self, m, v=None):
+        if v is not None:
+            return self.P * v
         return self.P
 
 
@@ -484,9 +560,9 @@ class InjectActiveCells(IdentityMap):
 
     """
 
-    indActive   = None #: Active Cells
-    valInactive = None #: Values of inactive Cells
-    nC          = None #: Number of cells in the full model
+    indActive   = None  #: Active Cells
+    valInactive = None  #: Values of inactive Cells
+    nC          = None  #: Number of cells in the full model
 
     def __init__(self, mesh, indActive, valInactive, nC=None):
         self.mesh  = mesh
@@ -494,7 +570,7 @@ class InjectActiveCells(IdentityMap):
         self.nC = nC or mesh.nC
 
         if indActive.dtype is not bool:
-            z = np.zeros(self.nC,dtype=bool)
+            z = np.zeros(self.nC, dtype=bool)
             z[indActive] = True
             indActive = z
         self.indActive = indActive
@@ -502,11 +578,15 @@ class InjectActiveCells(IdentityMap):
         if Utils.isScalar(valInactive):
             self.valInactive = np.ones(self.nC)*float(valInactive)
         else:
-            self.valInactive = valInactive.copy()
+            self.valInactive = np.ones(self.nC)
+            self.valInactive[self.indInactive] = valInactive.copy()
+
         self.valInactive[self.indActive] = 0
 
         inds = np.nonzero(self.indActive)[0]
-        self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
+        self.P = sp.csr_matrix((np.ones(inds.size), (inds, range(inds.size))),
+                                shape=(self.nC, self.nP)
+                               )
 
     @property
     def shape(self):
@@ -518,18 +598,25 @@ class InjectActiveCells(IdentityMap):
         return self.indActive.sum()
 
     def _transform(self, m):
-        return self.P*m + self.valInactive
+        return self.P * m + self.valInactive
 
     def inverse(self, D):
         return self.P.T*D
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
+        if v is not None:
+            return self.P * v
         return self.P
 
+
 class ActiveCells(InjectActiveCells):
+    """ActiveCells is depreciated. Use InjectActiveCells instead.
+    """
+
     def __init__(self, mesh, indActive, valInactive, nC=None):
         warnings.warn(
-            "`ActiveCells` is deprecated and will be removed in future versions. Use `InjectActiveCells` instead",
+            "`ActiveCells` is deprecated and will be removed in future "
+            "versions. Use `InjectActiveCells` instead",
             FutureWarning)
         InjectActiveCells.__init__(self, mesh, indActive, valInactive, nC)
 
@@ -537,11 +624,10 @@ class ActiveCells(InjectActiveCells):
 class Weighting(IdentityMap):
     """
         Model weight parameters.
-
     """
 
-    weights     = None #: Active Cells
-    nC          = None #: Number of cells in the full model
+    weights = None  #: Active Cells
+    nC = None  #: Number of cells in the full model
 
     def __init__(self, mesh, weights=None, nC=None):
         self.mesh  = mesh
@@ -571,7 +657,9 @@ class Weighting(IdentityMap):
         Pinv = Utils.sdiag(self.weights**(-1.))
         return Pinv*D
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
+        if v is not None:
+            return self.P * v
         return self.P
 
 
@@ -599,26 +687,32 @@ class ComplexMap(IdentityMap):
         nC = self.mesh.nC
         return m[:nC] + m[nC:]*1j
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
         nC = self.nP/2
         shp = (nC, nC*2)
+
         def fwd(v):
             return v[:nC] + v[nC:]*1j
+
         def adj(v):
-            return np.r_[v.real,v.imag]
-        return LinearOperator(shp,matvec=fwd,rmatvec=adj)
+            return np.r_[v.real, v.imag]
+        if v is not None:
+            return LinearOperator(shp, matvec=fwd, rmatvec=adj) * v
+        return LinearOperator(shp, matvec=fwd, rmatvec=adj)
 
     inverse = deriv
 
 
-class CircleMap(IdentityMap):
-    """CircleMap
+class ParametricCircleMap(IdentityMap):
+    """ParametricCircleMap
 
         Parameterize the model space using a circle in a wholespace.
 
         ..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)
+            \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:
 
@@ -627,46 +721,69 @@ class CircleMap(IdentityMap):
             m = [\sigma_1, \sigma_2, x_0, y_0, r]
 
     """
-    def __init__(self, mesh, logSigma=True):
-        assert mesh.dim == 2, "Working for a 2D mesh only right now. But it isn't that hard to change.. :)"
-        IdentityMap.__init__(self, mesh)
-        self.logSigma = logSigma
 
     slope = 1e-1
 
+    def __init__(self, mesh, logSigma=True):
+        assert mesh.dim == 2, "Working for a 2D mesh only right now. "
+        "But it isn't that hard to change.. :)"
+        IdentityMap.__init__(self, mesh)
+        # TODO: this should be done through a composition with and ExpMap
+        self.logSigma = logSigma
+
     @property
     def nP(self):
         return 5
 
     def _transform(self, m):
         a = self.slope
-        sig1,sig2,x,y,r = m[0],m[1],m[2],m[3],m[4]
+        sig1, sig2, x, y, r = m[0], m[1], m[2], m[3], m[4]
         if self.logSigma:
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
-        X = self.mesh.gridCC[:,0]
-        Y = self.mesh.gridCC[:,1]
-        return sig1 + (sig2 - sig1)*(np.arctan(a*(np.sqrt((X-x)**2 + (Y-y)**2) - r))/np.pi + 0.5)
+        X = self.mesh.gridCC[:, 0]
+        Y = self.mesh.gridCC[:, 1]
+        return sig1 + (sig2 - sig1)*(np.arctan(a*(np.sqrt((X-x)**2 +
+                                     (Y-y)**2) - r))/np.pi + 0.5)
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
         a = self.slope
-        sig1,sig2,x,y,r = m[0],m[1],m[2],m[3],m[4]
+        sig1, sig2, x, y, r = m[0], m[1], m[2], m[3], m[4]
         if self.logSigma:
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
-        X = self.mesh.gridCC[:,0]
-        Y = self.mesh.gridCC[:,1]
+        X = self.mesh.gridCC[:, 0]
+        Y = self.mesh.gridCC[:, 1]
         if self.logSigma:
-            g1 = -(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5)*sig1 + sig1
-            g2 = (np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5)*sig2
+            g1 = -(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi +
+                   0.5)*sig1 + sig1
+            g2 = (np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi +
+                  0.5)*sig2
         else:
-            g1 = -(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5) + 1.0
-            g2 = (np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi + 0.5)
+            g1 = -(np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi +
+                   0.5) + 1.0
+            g2 = (np.arctan(a*(-r + np.sqrt((X - x)**2 + (Y - y)**2)))/np.pi +
+                  0.5)
         g3 = a*(-X + x)*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1)*np.sqrt((X - x)**2 + (Y - y)**2))
         g4 = a*(-Y + y)*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1)*np.sqrt((X - x)**2 + (Y - y)**2))
         g5 = -a*(-sig1 + sig2)/(np.pi*(a**2*(-r + np.sqrt((X - x)**2 + (Y - y)**2))**2 + 1))
-        return sp.csr_matrix(np.c_[g1,g2,g3,g4,g5])
+
+        if v is not None:
+            return sp.csr_matrix(np.c_[g1, g2, g3, g4, g5]) * v
+        return sp.csr_matrix(np.c_[g1, g2, g3, g4, g5])
 
 
-class PolyMap(IdentityMap):
+class CircleMap(ParametricCircleMap):
+    """CircleMap is depreciated. Use ParametricCircleMap instead.
+    """
+
+    def __init__(self, mesh, logSigma=True):
+        warnings.warn(
+            "`CircleMap` is deprecated and will be removed in future "
+            "versions. Use `ParametricCircleMap` instead",
+            FutureWarning)
+        ParametricCircleMap.__init__(self, mesh, logSigma)
+
+
+class ParametricPolyMap(IdentityMap):
 
     """PolyMap
 
@@ -685,7 +802,8 @@ class PolyMap(IdentityMap):
         Can take in an actInd vector to account for topography.
 
     """
-    def __init__(self, mesh, order, logSigma=True, normal='X', actInd = None):
+
+    def __init__(self, mesh, order, logSigma=True, normal='X', actInd=None):
         IdentityMap.__init__(self, mesh)
         self.logSigma = logSigma
         self.order = order
@@ -710,78 +828,88 @@ class PolyMap(IdentityMap):
         if np.isscalar(self.order):
             nP = self.order+3
         else:
-            nP =(self.order[0]+1)*(self.order[1]+1)+2
+            nP = (self.order[0]+1)*(self.order[1]+1)+2
         return nP
 
     def _transform(self, m):
         # Set model parameters
         alpha = self.slope
-        sig1,sig2 = m[0],m[1]
+        sig1, sig2 = m[0], m[1]
         c = m[2:]
         if self.logSigma:
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
-        #2D
+
+        # 2D
         if self.mesh.dim == 2:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
+            X = self.mesh.gridCC[self.actInd, 0]
+            Y = self.mesh.gridCC[self.actInd, 1]
             if self.normal =='X':
                 f = polynomial.polyval(Y, c) - X
             elif self.normal =='Y':
                 f = polynomial.polyval(X, c) - Y
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
-        #3D
+
+        # 3D
         elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
-            Z = self.mesh.gridCC[self.actInd,2]
-            if self.normal =='X':
-                f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
-            elif self.normal =='Y':
-                f = polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - Y
-            elif self.normal =='Z':
-                f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
+            X = self.mesh.gridCC[self.actInd, 0]
+            Y = self.mesh.gridCC[self.actInd, 1]
+            Z = self.mesh.gridCC[self.actInd, 2]
+
+            if self.normal == 'X':
+                f = (polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,
+                     self.order[1]+1))) - X)
+            elif self.normal == 'Y':
+                f = (polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,
+                     self.order[1]+1))) - Y)
+            elif self.normal == 'Z':
+                f = (polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,
+                     self.order[1]+1))) - Z)
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
 
         else:
             raise(Exception("Only supports 2D"))
 
-
         return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
         alpha = self.slope
-        sig1,sig2, c = m[0],m[1],m[2:]
+        sig1, sig2, c = m[0], m[1], m[2:]
         if self.logSigma:
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
-        #2D
-        if self.mesh.dim == 2:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
 
-            if self.normal =='X':
+        # 2D
+        if self.mesh.dim == 2:
+            X = self.mesh.gridCC[self.actInd, 0]
+            Y = self.mesh.gridCC[self.actInd, 1]
+
+            if self.normal == 'X':
                 f = polynomial.polyval(Y, c) - X
                 V = polynomial.polyvander(Y, len(c)-1)
-            elif self.normal =='Y':
+            elif self.normal == 'Y':
                 f = polynomial.polyval(X, c) - Y
                 V = polynomial.polyvander(X, len(c)-1)
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
-        #3D
-        elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
-            Z = self.mesh.gridCC[self.actInd,2]
 
-            if self.normal =='X':
-                f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
+        # 3D
+        elif self.mesh.dim == 3:
+            X = self.mesh.gridCC[self.actInd, 0]
+            Y = self.mesh.gridCC[self.actInd, 1]
+            Z = self.mesh.gridCC[self.actInd, 2]
+
+            if self.normal == 'X':
+                f = (polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,
+                     self.order[1]+1))) - X)
                 V = polynomial.polyvander2d(Y, Z, self.order)
-            elif self.normal =='Y':
-                f = polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - Y
+            elif self.normal == 'Y':
+                f = (polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,
+                     self.order[1]+1))) - Y)
                 V = polynomial.polyvander2d(X, Z, self.order)
-            elif self.normal =='Z':
-                f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
+            elif self.normal == 'Z':
+                f = (polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,
+                     self.order[1]+1))) - Z)
                 V = polynomial.polyvander2d(X, Y, self.order)
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
@@ -795,13 +923,17 @@ class PolyMap(IdentityMap):
 
         g3 = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*V
 
-        return sp.csr_matrix(np.c_[g1,g2,g3])
+        if v is not None:
+            return sp.csr_matrix(np.c_[g1, g2, g3]) * v
+        return sp.csr_matrix(np.c_[g1, g2, g3])
 
-class SplineMap(IdentityMap):
+
+class ParametricSplineMap(IdentityMap):
 
     """SplineMap
 
-        Parameterize the boundary of two geological units using a spline interpolation
+        Parameterize the boundary of two geological units using
+        a spline interpolation
 
         ..math::
 
@@ -814,7 +946,10 @@ class SplineMap(IdentityMap):
             m = [\sigma_1, \sigma_2, y]
 
     """
-    def __init__(self, mesh, pts, ptsv=None,order=3, logSigma=True, normal='X'):
+
+    slope = 1e4
+
+    def __init__(self, mesh, pts, ptsv=None, order=3, logSigma=True, normal='X'):
         IdentityMap.__init__(self, mesh)
         self.logSigma = logSigma
         self.order = order
@@ -824,7 +959,6 @@ class SplineMap(IdentityMap):
         self.ptsv = ptsv
         self.spl = None
 
-    slope = 1e4
     @property
     def nP(self):
         if self.mesh.dim == 2:
@@ -837,18 +971,18 @@ class SplineMap(IdentityMap):
     def _transform(self, m):
         # Set model parameters
         alpha = self.slope
-        sig1,sig2 = m[0],m[1]
+        sig1, sig2 = m[0], m[1]
         c = m[2:]
         if self.logSigma:
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
-        #2D
+        # 2D
         if self.mesh.dim == 2:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
+            X = self.mesh.gridCC[:, 0]
+            Y = self.mesh.gridCC[:, 1]
             self.spl = UnivariateSpline(self.pts, c, k=self.order, s=0)
-            if self.normal =='X':
+            if self.normal == 'X':
                 f = self.spl(Y) - X
-            elif self.normal =='Y':
+            elif self.normal == 'Y':
                 f = self.spl(X) - Y
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
@@ -860,18 +994,18 @@ class SplineMap(IdentityMap):
         # Using 2D interpolation  is possible
 
         elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
-            Z = self.mesh.gridCC[:,2]
+            X = self.mesh.gridCC[:, 0]
+            Y = self.mesh.gridCC[:, 1]
+            Z = self.mesh.gridCC[:, 2]
 
             npts = np.size(self.pts)
             if np.mod(c.size, 2):
                 raise(Exception("Put even points!"))
 
-            self.spl = {"splb":UnivariateSpline(self.pts, c[:npts], k=self.order, s=0),
-                        "splt":UnivariateSpline(self.pts, c[npts:], k=self.order, s=0)}
+            self.spl = {"splb": UnivariateSpline(self.pts, c[:npts], k=self.order, s=0),
+                        "splt": UnivariateSpline(self.pts, c[npts:], k=self.order, s=0)}
 
-            if self.normal =='X':
+            if self.normal == 'X':
                 zb = self.ptsv[0]
                 zt = self.ptsv[1]
                 flines = (self.spl["splt"](Y)-self.spl["splb"](Y))*(Z-zb)/(zt-zb) + self.spl["splb"](Y)
@@ -883,30 +1017,29 @@ class SplineMap(IdentityMap):
         else:
             raise(Exception("Only supports 2D and 3D"))
 
-
         return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
 
-    def deriv(self, m):
+    def deriv(self, m, v=None):
         alpha = self.slope
-        sig1,sig2, c = m[0],m[1],m[2:]
+        sig1, sig2,  c = m[0], m[1], m[2:]
         if self.logSigma:
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
         #2D
         if self.mesh.dim == 2:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
+            X = self.mesh.gridCC[:, 0]
+            Y = self.mesh.gridCC[:, 1]
 
-            if self.normal =='X':
+            if self.normal == 'X':
                 f = self.spl(Y) - X
-            elif self.normal =='Y':
+            elif self.normal == 'Y':
                 f = self.spl(X) - Y
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
         #3D
         elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
-            Z = self.mesh.gridCC[:,2]
+            X = self.mesh.gridCC[:, 0]
+            Y = self.mesh.gridCC[:, 1]
+            Z = self.mesh.gridCC[:, 2]
             if self.normal =='X':
                 zb = self.ptsv[0]
                 zt = self.ptsv[1]
@@ -924,10 +1057,9 @@ class SplineMap(IdentityMap):
             g1 = -(np.arctan(alpha*f)/np.pi + 0.5) + 1.0
             g2 = (np.arctan(alpha*f)/np.pi + 0.5)
 
-
-        if self.mesh.dim ==2:
+        if self.mesh.dim == 2:
             g3 = np.zeros((self.mesh.nC, self.npts))
-            if self.normal =='Y':
+            if self.normal == 'Y':
                 # Here we use perturbation to compute sensitivity
                 # TODO: bit more generalization of this ...
                 # Modfications for X and Z directions ...
@@ -942,11 +1074,11 @@ class SplineMap(IdentityMap):
                     spla = UnivariateSpline(self.pts, ca, k=self.order, s=0)
                     splb = UnivariateSpline(self.pts, cb, k=self.order, s=0)
                     fderiv = (spla(X)-splb(X))/(2*dy)
-                    g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
+                    g3[:, i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
 
-        elif self.mesh.dim==3:
+        elif self.mesh.dim == 3:
             g3 = np.zeros((self.mesh.nC, self.npts*2))
-            if self.normal =='X':
+            if self.normal == 'X':
                 # Here we use perturbation to compute sensitivity
                 for i in range(self.npts*2):
                     ctemp = c[i]
@@ -956,26 +1088,40 @@ class SplineMap(IdentityMap):
                     dy = self.mesh.hy[ind]*1.5
                     ca[i] = ctemp+dy
                     cb[i] = ctemp-dy
-                    #treat bottom boundary
-                    if i< self.npts:
+
+                    # treat bottom boundary
+                    if i < self.npts:
                         splba = UnivariateSpline(self.pts, ca[:self.npts], k=self.order, s=0)
                         splbb = UnivariateSpline(self.pts, cb[:self.npts], k=self.order, s=0)
                         flinesa = (self.spl["splt"](Y)-splba(Y))*(Z-zb)/(zt-zb) + splba(Y) - X
                         flinesb = (self.spl["splt"](Y)-splbb(Y))*(Z-zb)/(zt-zb) + splbb(Y) - X
-                    #treat top boundary
+
+                    # treat top boundary
                     else:
                         splta = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
                         spltb = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
                         flinesa = (self.spl["splt"](Y)-splta(Y))*(Z-zb)/(zt-zb) + splta(Y) - X
                         flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X
                     fderiv = (flinesa-flinesb)/(2*dy)
-                    g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
+                    g3[:, i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
         else :
             raise(Exception("Not Implemented for Y and Z, your turn :)"))
-        return sp.csr_matrix(np.c_[g1,g2,g3])
 
-
+        if v is not None:
+            return sp.csr_matrix(np.c_[g1, g2, g3]) * v
+        return sp.csr_matrix(np.c_[g1, g2, g3])
 
 
+class SplineMap(ParametricSplineMap):
+    """SplineMap is depreciated. Use ParametricSplineMap instead.
+    """
 
+    def __init__(self, mesh, pts, ptsv=None, order=3, logSigma=True,
+                 normal='X'):
+        warnings.warn(
+            "`SplineMap` is deprecated and will be removed in future "
+            "versions. Use `ParametricSplineMap` instead",
+            FutureWarning)
+        ParametricSplineMap.__init__(self, mesh, pts, ptsv, order, logSigma,
+                                     normal)
 

SimPEG/__init__.py

@@ -15,7 +15,7 @@ import Directives
 import Inversion
 import Tests
 
-__version__   = '0.1.11'
+__version__   = '0.1.12'
 __author__    = 'Rowan Cockett'
 __license__   = 'MIT'
 __copyright__ = 'Copyright 2014 Rowan Cockett'

docs/conf.py

@@ -51,9 +51,9 @@ copyright = u'2013 - 2016, SimPEG Developers'
 # built documents.
 #
 # The short X.Y version.
-version = '0.1.11'
+version = '0.1.12'
 # The full version, including alpha/beta/rc tags.
-release = '0.1.11'
+release = '0.1.12'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

docs/content/api_core/api_Maps.rst

@@ -63,26 +63,8 @@ done by the :class:`SimPEG.Maps.ExpMap` described above.
 
 .. plot::
 
-    from SimPEG import *
-    import matplotlib.pyplot as plt
-    M = Mesh.TensorMesh([7,5])
-    v1dMap = Maps.SurjectVertical1D(M)
-    expMap = Maps.ExpMap(M)
-    myMap = expMap * v1dMap
-    m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
-    sig = myMap * m
-    figs, axs = plt.subplots(1,2)
-    axs[0].plot(m, M.vectorCCy, 'b-o')
-    axs[0].set_title('Model')
-    axs[0].set_ylabel('Depth, y')
-    axs[0].set_xlabel('Value, $m_i$')
-    axs[0].set_xlim(0,3)
-    axs[0].set_ylim(0,1)
-    clbar = plt.colorbar(M.plotImage(sig,ax=axs[1],grid=True,gridOpts=dict(color='grey'))[0])
-    axs[1].set_title('Physical Property')
-    axs[1].set_ylabel('Depth, y')
-    clbar.set_label('$\sigma = \exp(\mathbf{P}m)$')
-    plt.tight_layout()
+    from SimPEG import Examples
+    Examples.Maps_ComboMaps.run()
 
 If you noticed, it was pretty easy to combine maps. What is even cooler is
 that the derivatives also are made for you (if everything goes right).
@@ -167,31 +149,10 @@ Map 2D Cross-Section to 3D Model
 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()
+    from SimPEG import Examples
+    Examples.Maps_Mesh2Mesh.run()
 
 
 .. autoclass:: SimPEG.Maps.Mesh2Mesh

docs/content/dc/index.rst

@@ -40,7 +40,7 @@
 Direct Current Resistivity
 **************************
 
-`SimPEG.EM.Static.DC` and `SimPEG.EM.Static.IP` uses SimPEG as the framework for the forward and inverse
+`SimPEG.DCIP` uses SimPEG as the framework for the forward and inverse
 direct current (DC) resistivity and induced polarization (IP) geophysical problems.
 
 
@@ -150,7 +150,7 @@ Comparing to the analytic function:
 API for DC codes
 ================
 
-.. automodule:: SimPEG.EM.Static.DC
+.. automodule:: SimPEG.DCIP.BaseDC
     :show-inheritance:
     :members:
     :undoc-members:

docs/content/examples/Inversion_IRLS.rst

@@ -0,0 +1,26 @@
+.. _examples_Inversion_IRLS:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+Inversion: Linear Problem
+=========================
+
+Here we go over the basics of creating a linear problem and inversion.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.Inversion_IRLS.run()
+
+.. literalinclude:: ../../../SimPEG/Examples/Inversion_IRLS.py
+    :language: python
+    :linenos:

docs/content/examples/Maps_ComboMaps.rst

@@ -0,0 +1,48 @@
+.. _examples_Maps_ComboMaps:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+
+Maps: ComboMaps
+===============
+
+We will use an example where we want a 1D layered earth as
+our model, but we want to map this to a 2D discretization to do our forward
+modeling. We will also assume that we are working in log conductivity still,
+so after the transformation we want to map to conductivity space.
+To do this we will introduce the vertical 1D map (:class:`SimPEG.Maps.SurjectVertical1D`),
+which does the first part of what we just described. The second part will be
+done by the :class:`SimPEG.Maps.ExpMap` described above.
+
+.. code-block:: python
+    :linenos:
+
+    M = Mesh.TensorMesh([7,5])
+    v1dMap = Maps.SurjectVertical1D(M)
+    expMap = Maps.ExpMap(M)
+    myMap = expMap * v1dMap
+    m = np.r_[0.2,1,0.1,2,2.9] # only 5 model parameters!
+    sig = myMap * m
+
+If you noticed, it was pretty easy to combine maps. What is even cooler is
+that the derivatives also are made for you (if everything goes right).
+Just to be sure that the derivative is correct, you should always run the test
+on the mapping that you create.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.Maps_ComboMaps.run()
+
+.. literalinclude:: ../../../SimPEG/Examples/Maps_ComboMaps.py
+    :language: python
+    :linenos:

docs/content/examples/Maps_Mesh2Mesh.rst

@@ -0,0 +1,27 @@
+.. _examples_Maps_Mesh2Mesh:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+
+Maps: Mesh2Mesh
+===============
+
+This mapping allows you to go from one mesh to another.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.Maps_Mesh2Mesh.run()
+
+.. literalinclude:: ../../../SimPEG/Examples/Maps_Mesh2Mesh.py
+    :language: python
+    :linenos:

docs/content/ip/index.rst

@@ -7,7 +7,7 @@ Todo: docs for IP!
 API for IP codes
 ================
 
-.. automodule:: SimPEG.EM.Static.IP
+.. automodule:: SimPEG.DCIP.BaseIP
     :show-inheritance:
     :members:
     :undoc-members:

docs/index.rst

@@ -1,4 +1,4 @@
-.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
+.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png
     :alt: SimPEG Logo
 
 SimPEG Documentation

setup.py

@@ -83,7 +83,7 @@ with open("README.rst") as f:
 
 setup(
     name = "SimPEG",
-    version = "0.1.11",
+    version = "0.1.12",
     packages = find_packages(),
     install_requires = ['numpy>=1.7',
                         'scipy>=0.13',

tests/base/test_maps.py

@@ -2,39 +2,71 @@ import numpy as np
 import unittest
 from SimPEG import *
 from scipy.sparse.linalg import dsolve
+import inspect
 
 TOL = 1e-14
 
-MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull","FullMap"]
-MAPS_TO_TEST_3D = [             "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull","FullMap"]
+MAPS_TO_EXCLUDE_2D = ["ComboMap", "ActiveCells", "InjectActiveCells"]
+MAPS_TO_EXCLUDE_3D = ["ComboMap", "ActiveCells", "InjectActiveCells",
+                      "CircleMap"]
+
+MAPS_TO_TEST_INVERSE = ["ExpMap"]
+
 
 class MapTests(unittest.TestCase):
 
     def setUp(self):
 
+        maps2test2D = [M for M in dir(Maps) if M not in MAPS_TO_EXCLUDE_2D]
+        maps2test3D = [M for M in dir(Maps) if M not in MAPS_TO_EXCLUDE_3D]
+
+        self.maps2test2D = [getattr(Maps, m) for m in maps2test2D if
+                            inspect.isclass(getattr(Maps, M)) and
+                            issubclass(getattr(Maps, M), Maps.IdentityMap)]
+
+        self.maps2test3D = [getattr(Maps, m) for m in maps2test3D if
+                            inspect.isclass(getattr(Maps, M)) and
+                            issubclass(getattr(Maps, M), Maps.IdentityMap)]
+
         a = np.array([1, 1, 1])
         b = np.array([1, 2])
+
         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.mesh3 = Mesh.TensorMesh([a, b, [3, 4]], x0=np.array([3, 5, 2]))
         self.mesh22 = Mesh.TensorMesh([b, a], x0=np.array([3, 5]))
 
     def test_transforms2D(self):
-        for M in MAPS_TO_TEST_2D:
-            maps = getattr(Maps, M)(self.mesh2)
-            self.assertTrue(maps.test())
+        for M in self.maps2test2D:
+            self.assertTrue(M.test())
+
+    def test_transforms2Dvec(self):
+        for M in self.maps2test2D:
+            self.assertTrue(M.testVec())
 
     def test_transforms3D(self):
-        for M in MAPS_TO_TEST_3D:
-            maps = getattr(Maps, M)(self.mesh3)
-            self.assertTrue(maps.test())
+        for M in self.maps2test3D:
+            self.assertTrue(M.test())
 
+    def test_transforms3Dvec(self):
+        for M in self.maps2test3D:
+            self.assertTrue(M.test())
 
     def test_transforms_logMap_reciprocalMap(self):
-        # Note that log/reciprocal maps can be kinda finicky, so we are being explicit about the random seed.
-        v2 = np.r_[ 0.40077291, 0.14410044, 0.58452314, 0.96323738, 0.01198519, 0.79754415]
-        dv2 = np.r_[ 0.80653921, 0.13132446, 0.4901117, 0.03358737, 0.65473762, 0.44252488]
-        v3 = np.r_[ 0.96084865, 0.34385186, 0.39430044, 0.81671285, 0.65929109, 0.2235217, 0.87897526, 0.5784033, 0.96876393, 0.63535864, 0.84130763, 0.22123854]
-        dv3 = np.r_[ 0.96827838, 0.26072111, 0.45090749, 0.10573893, 0.65276365, 0.15646586, 0.51679682, 0.23071984, 0.95106218, 0.14201845, 0.25093564, 0.3732866 ]
+
+        # Note that log/reciprocal maps can be kinda finicky, so we are being
+        # explicit about the random seed.
+
+        v2 = np.r_[0.40077291, 0.14410044, 0.58452314, 0.96323738, 0.01198519,
+                   0.79754415]
+        dv2 = np.r_[0.80653921, 0.13132446, 0.4901117, 0.03358737, 0.65473762,
+                    0.44252488]
+        v3 = np.r_[0.96084865, 0.34385186, 0.39430044, 0.81671285, 0.65929109,
+                   0.2235217, 0.87897526, 0.5784033, 0.96876393, 0.63535864,
+                   0.84130763, 0.22123854]
+        dv3 = np.r_[0.96827838, 0.26072111, 0.45090749, 0.10573893,
+                    0.65276365, 0.15646586, 0.51679682, 0.23071984,
+                    0.95106218, 0.14201845, 0.25093564, 0.3732866 ]
+
         maps = Maps.LogMap(self.mesh2)
         self.assertTrue(maps.test(v2, dx=dv2))
         maps = Maps.LogMap(self.mesh3)
@@ -49,100 +81,126 @@ class MapTests(unittest.TestCase):
         maps = Maps.Mesh2Mesh([self.mesh22, self.mesh2])
         self.assertTrue(maps.test())
 
+    def test_Mesh2MeshMapVec(self):
+        maps = Maps.Mesh2Mesh([self.mesh22, self.mesh2])
+        self.assertTrue(maps.testVec())
+
     def test_mapMultiplication(self):
-        M = Mesh.TensorMesh([2,3])
+        M = Mesh.TensorMesh([2, 3])
         expMap = Maps.ExpMap(M)
         vertMap = Maps.SurjectVertical1D(M)
         combo = expMap*vertMap
         m = np.arange(3.0)
-        t_true = np.exp(np.r_[0,0,1,1,2,2.])
-        self.assertLess(np.linalg.norm((combo * m)-t_true,np.inf),TOL)
-        self.assertLess(np.linalg.norm((expMap * vertMap * m)-t_true,np.inf),TOL)
-        self.assertLess(np.linalg.norm(expMap * (vertMap * m)-t_true,np.inf),TOL)
-        self.assertLess(np.linalg.norm((expMap * vertMap) * m-t_true,np.inf),TOL)
-        #Try making a model
+        t_true = np.exp(np.r_[0, 0, 1, 1, 2, 2.])
+        self.assertLess(np.linalg.norm((combo * m) - t_true, np.inf), TOL)
+        self.assertLess(np.linalg.norm((expMap * vertMap * m)-t_true, np.inf),
+                        TOL)
+        self.assertLess(np.linalg.norm(expMap * (vertMap * m)-t_true, np.inf),
+                        TOL)
+        self.assertLess(np.linalg.norm((expMap * vertMap) * m-t_true, np.inf),
+                        TOL)
+        # Try making a model
         mod = Models.Model(m, mapping=combo)
         # print mod.transform
         # import matplotlib.pyplot as plt
         # plt.colorbar(M.plotImage(mod.transform)[0])
         # plt.show()
-        self.assertLess(np.linalg.norm(mod.transform-t_true,np.inf),TOL)
+        self.assertLess(np.linalg.norm(mod.transform - t_true, np.inf), TOL)
 
-        self.assertRaises(Exception,Models.Model,np.r_[1.0],mapping=combo)
+        self.assertRaises(Exception, Models.Model, np.r_[1.0], mapping=combo)
 
         self.assertRaises(ValueError, lambda: combo * (vertMap * expMap))
         self.assertRaises(ValueError, lambda: (combo * vertMap) * expMap)
         self.assertRaises(ValueError, lambda: vertMap * expMap)
         self.assertRaises(ValueError, lambda: expMap * np.ones(100))
-        self.assertRaises(ValueError, lambda: expMap * np.ones((100.0,1)))
-        self.assertRaises(ValueError, lambda: expMap * np.ones((100.0,5)))
+        self.assertRaises(ValueError, lambda: expMap * np.ones((100.0, 1)))
+        self.assertRaises(ValueError, lambda: expMap * np.ones((100.0, 5)))
         self.assertRaises(ValueError, lambda: combo * np.ones(100))
-        self.assertRaises(ValueError, lambda: combo * np.ones((100.0,1)))
-        self.assertRaises(ValueError, lambda: combo * np.ones((100.0,5)))
+        self.assertRaises(ValueError, lambda: combo * np.ones((100.0, 1)))
+        self.assertRaises(ValueError, lambda: combo * np.ones((100.0, 5)))
 
     def test_activeCells(self):
-        M = Mesh.TensorMesh([2,4],'0C')
+        M = Mesh.TensorMesh([2, 4], '0C')
         expMap = Maps.ExpMap(M)
-        for actMap in [Maps.InjectActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy), Maps.ActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy)]:
-        # actMap = Maps.InjectActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy)
+        for actMap in [Maps.InjectActiveCells(M, M.vectorCCy <= 0, 10,
+                       nC=M.nCy), Maps.ActiveCells(M, M.vectorCCy <= 0, 10,
+                       nC=M.nCy)]:
+
             vertMap = Maps.SurjectVertical1D(M)
             combo = vertMap * actMap
-            m = np.r_[1,2.]
-            mod = Models.Model(m,combo)
-            # import matplotlib.pyplot as plt
-            # plt.colorbar(M.plotImage(mod.transform)[0])
-            # plt.show()
-            self.assertLess(np.linalg.norm(mod.transform - np.r_[1,1,2,2,10,10,10,10.]), TOL)
-            self.assertLess((mod.transformDeriv - combo.deriv(m)).toarray().sum(), TOL)
+            m = np.r_[1., 2.]
+            mod = Models.Model(m, combo)
+
+            self.assertLess(np.linalg.norm(mod.transform -
+                            np.r_[1, 1, 2, 2, 10, 10, 10, 10.]), TOL)
+            self.assertLess((mod.transformDeriv -
+                            combo.deriv(m)).toarray().sum(), TOL)
 
     def test_tripleMultiply(self):
-        M = Mesh.TensorMesh([2,4],'0C')
+        M = Mesh.TensorMesh([2, 4], '0C')
         expMap = Maps.ExpMap(M)
         vertMap = Maps.SurjectVertical1D(M)
-        actMap = Maps.InjectActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy)
-        m = np.r_[1,2.]
-        t_true = np.exp(np.r_[1,1,2,2,10,10,10,10.])
-        self.assertLess(np.linalg.norm((expMap * vertMap * actMap * m)-t_true,np.inf),TOL)
-        self.assertLess(np.linalg.norm(((expMap * vertMap * actMap) * m)-t_true,np.inf),TOL)
-        self.assertLess(np.linalg.norm((expMap * vertMap * (actMap * m))-t_true,np.inf),TOL)
-        self.assertLess(np.linalg.norm((expMap * (vertMap * actMap) * m)-t_true,np.inf),TOL)
-        self.assertLess(np.linalg.norm(((expMap * vertMap) * actMap * m)-t_true,np.inf),TOL)
+        actMap = Maps.InjectActiveCells(M, M.vectorCCy <= 0, 10, nC=M.nCy)
+        m = np.r_[1., 2.]
+        t_true = np.exp(np.r_[1, 1, 2, 2, 10, 10, 10, 10.])
 
-        self.assertRaises(ValueError, lambda: expMap * actMap * vertMap )
-        self.assertRaises(ValueError, lambda: actMap * vertMap * expMap )
+        self.assertLess(np.linalg.norm((expMap * vertMap * actMap * m) -
+                        t_true, np.inf), TOL)
+        self.assertLess(np.linalg.norm(((expMap * vertMap * actMap) * m) -
+                        t_true, np.inf), TOL)
+        self.assertLess(np.linalg.norm((expMap * vertMap * (actMap * m)) -
+                        t_true, np.inf), TOL)
+        self.assertLess(np.linalg.norm((expMap * (vertMap * actMap) * m) -
+                        t_true, np.inf), TOL)
+        self.assertLess(np.linalg.norm(((expMap * vertMap) * actMap * m) -
+                        t_true, np.inf), TOL)
 
+        self.assertRaises(ValueError, lambda: expMap * actMap * vertMap)
+        self.assertRaises(ValueError, lambda: actMap * vertMap * expMap)
 
     def test_map2Dto3D_x(self):
-        M2 = Mesh.TensorMesh([2,4])
-        M3 = Mesh.TensorMesh([3,2,4])
+        M2 = Mesh.TensorMesh([2, 4])
+        M3 = Mesh.TensorMesh([3, 2, 4])
         m = np.random.rand(M2.nC)
 
-        for m2to3 in [Maps.Surject2Dto3D(M3, normal='X'), Maps.Map2Dto3D(M3, normal='X')]:
-        # m2to3 = Maps.Surject2Dto3D(M3, normal='X')
+        for m2to3 in [Maps.Surject2Dto3D(M3, normal='X'),
+                      Maps.Map2Dto3D(M3, normal='X')]:
+
+            # m2to3 = Maps.Surject2Dto3D(M3, normal='X')
             m = np.arange(m2to3.nP)
             self.assertTrue(m2to3.test())
-            self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[0,:,:] ) == m))
-
+            self.assertTrue(m2to3.testVec())
+            self.assertTrue(np.all(Utils.mkvc((m2to3 * m).reshape(M3.vnC,
+                            order='F')[0, :, :]) == m))
 
     def test_map2Dto3D_y(self):
-        M2 = Mesh.TensorMesh([3,4])
-        M3 = Mesh.TensorMesh([3,2,4])
+        M2 = Mesh.TensorMesh([3, 4])
+        M3 = Mesh.TensorMesh([3, 2, 4])
         m = np.random.rand(M2.nC)
-        for m2to3 in [Maps.Surject2Dto3D(M3, normal='Y'),Maps.Map2Dto3D(M3, normal='Y')]:
-        # m2to3 = Maps.Surject2Dto3D(M3, normal='Y')
+
+        for m2to3 in [Maps.Surject2Dto3D(M3, normal='Y'), Maps.Map2Dto3D(M3,
+                      normal='Y')]:
+            # m2to3 = Maps.Surject2Dto3D(M3, normal='Y')
             m = np.arange(m2to3.nP)
             self.assertTrue(m2to3.test())
-            self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[:,0,:] ) == m))
+            self.assertTrue(m2to3.testVec())
+            self.assertTrue(np.all(Utils.mkvc((m2to3 * m).reshape(M3.vnC,
+                            order='F')[:, 0, :]) == m))
 
     def test_map2Dto3D_z(self):
-        M2 = Mesh.TensorMesh([3,2])
-        M3 = Mesh.TensorMesh([3,2,4])
+        M2 = Mesh.TensorMesh([3, 2])
+        M3 = Mesh.TensorMesh([3, 2, 4])
         m = np.random.rand(M2.nC)
-        for m2to3 in [Maps.Surject2Dto3D(M3, normal='Z'),Maps.Map2Dto3D(M3, normal='Z')]:
-        # m2to3 = Maps.Surject2Dto3D(M3, normal='Z')
+
+        for m2to3 in [Maps.Surject2Dto3D(M3, normal='Z'), Maps.Map2Dto3D(M3,
+                      normal='Z')]:
+
+            # m2to3 = Maps.Surject2Dto3D(M3, normal='Z')
             m = np.arange(m2to3.nP)
             self.assertTrue(m2to3.test())
-            self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[:,:,0] ) == m))
+            self.assertTrue(m2to3.testVec())
+            self.assertTrue(np.all(Utils.mkvc((m2to3 * m).reshape(M3.vnC,
+                            order='F')[:, :, 0]) == m))
 
 
 if __name__ == '__main__':

tests/base/test_problem.py

@@ -5,12 +5,12 @@ from SimPEG import *
 class TestTimeProblem(unittest.TestCase):
 
     def setUp(self):
-        mesh = Mesh.TensorMesh([10,10])
+        mesh = Mesh.TensorMesh([10, 10])
         self.prob = Problem.BaseTimeProblem(mesh)
 
     def test_timeProblem_setTimeSteps(self):
         self.prob.timeSteps = [(1e-6, 3), 1e-5, (1e-4, 2)]
-        trueTS = np.r_[1e-6,1e-6,1e-6,1e-5,1e-4,1e-4]
+        trueTS = np.r_[1e-6, 1e-6, 1e-6, 1e-5, 1e-4, 1e-4]
         self.assertTrue(np.all(trueTS == self.prob.timeSteps))
 
         self.prob.timeSteps = trueTS
@@ -18,7 +18,7 @@ class TestTimeProblem(unittest.TestCase):
 
         self.assertTrue(self.prob.nT == 6)
 
-        self.assertTrue(np.all(self.prob.times == np.r_[0,trueTS].cumsum()))
+        self.assertTrue(np.all(self.prob.times == np.r_[0, trueTS].cumsum()))
 
 
 if __name__ == '__main__':

tests/base/test_regularization.py

@@ -8,6 +8,7 @@ TOL = 1e-20
 testReg = True
 testRegMesh = True
 
+
 class RegularizationTests(unittest.TestCase):
 
     def setUp(self):
@@ -16,41 +17,47 @@ class RegularizationTests(unittest.TestCase):
         mesh1 = Mesh.TensorMesh([hx])
         mesh2 = Mesh.TensorMesh([hx, hy])
         mesh3 = Mesh.TensorMesh([hx, hy, hz])
-        self.meshlist = [mesh1,mesh2, mesh3]
+        self.meshlist = [mesh1, mesh2, mesh3]
 
     if testReg:
         def test_regularization(self):
             for R in dir(Regularization):
                 r = getattr(Regularization, R)
-                if not inspect.isclass(r): continue
+                if not inspect.isclass(r):
+                    continue
                 if not issubclass(r, Regularization.BaseRegularization):
                     continue
 
                 for i, mesh in enumerate(self.meshlist):
 
-                    print 'Testing %iD'%mesh.dim
+                    print 'Testing %iD' % mesh.dim
 
                     mapping = r.mapPair(mesh)
                     reg = r(mesh, mapping=mapping)
                     m = np.random.rand(mapping.nP)
                     reg.mref = np.ones_like(m)*np.mean(m)
 
-                    print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                    print 'Check: phi_m (mref) = %f' % reg.eval(reg.mref)
                     passed = reg.eval(reg.mref) < TOL
                     self.assertTrue(passed)
 
                     print 'Check:', R
-                    passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                    passed = Tests.checkDerivative(lambda m: [reg.eval(m),
+                                                   reg.evalDeriv(m)], m,
+                                                   plotIt=False)
                     self.assertTrue(passed)
 
                     print 'Check 2 Deriv:', R
-                    passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                    passed = Tests.checkDerivative(lambda m: [reg.evalDeriv(m),
+                                                   reg.eval2Deriv(m)], m,
+                                                   plotIt=False)
                     self.assertTrue(passed)
 
         def test_regularization_ActiveCells(self):
             for R in dir(Regularization):
                 r = getattr(Regularization, R)
-                if not inspect.isclass(r): continue
+                if not inspect.isclass(r):
+                    continue
                 if not issubclass(r, Regularization.BaseRegularization):
                     continue
 

tests/dcip/__init__.py

@@ -0,0 +1,12 @@
+import os
+import glob
+import unittest
+
+if __name__ == '__main__':
+    test_file_strings = glob.glob('test_*.py')
+    module_strings = [str[0:len(str)-3] for str in test_file_strings]
+    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
+              in module_strings]
+    testSuite = unittest.TestSuite(suites)
+
+    unittest.TextTestRunner(verbosity=2).run(testSuite)

tests/dcip/test_forward_DCproblem.py

@@ -0,0 +1,77 @@
+import unittest
+from SimPEG import *
+import SimPEG.DCIP as DC
+
+
+class DCProblemTests(unittest.TestCase):
+
+    def setUp(self):
+
+        aSpacing=2.5
+        nElecs=10
+
+        surveySize = nElecs*aSpacing - aSpacing
+        cs = surveySize/nElecs/4
+
+        mesh = Mesh.TensorMesh([
+                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
+                [(cs,3, -1.3),(cs,3,1.3)],
+        #         [(cs,5, -1.3),(cs,10)]
+            ],'CN')
+
+        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
+        survey = DC.SurveyDC(srcList)
+        problem = DC.ProblemDC_CC(mesh)
+        problem.pair(survey)
+
+        mSynth = np.ones(mesh.nC)
+        survey.makeSyntheticData(mSynth)
+
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p = problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-10
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+
+    def test_massMatrices(self):
+        Gu = np.random.rand(self.mesh.nF)
+        def derChk(m):
+            self.p.curModel = m
+            return [self.p.Msig * Gu, self.p.dMdsig(Gu)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/dcip/test_forward_IPproblem.py

@@ -0,0 +1,65 @@
+import unittest
+import SimPEG.DCIP as DC
+from SimPEG import *
+
+class IPforwardTests(unittest.TestCase):
+
+    def test_IPforward(self):
+
+        cs = 12.5
+        nc = 200/cs+1
+        hx = [(cs,7, -1.3),(cs,nc),(cs,7, 1.3)]
+        hy = [(cs,7, -1.3),(cs,int(nc/2+1)),(cs,7, 1.3)]
+        hz = [(cs,7, -1.3),(cs,int(nc/2+1))]
+        mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
+        sighalf = 1e-2
+        sigma = np.ones(mesh.nC)*sighalf
+        p0 = np.r_[-50., 50., -50.]
+        p1 = np.r_[ 50.,-50., -150.]
+        blk_ind = Utils.ModelBuilder.getIndicesBlock(p0, p1, mesh.gridCC)
+        sigma[blk_ind] = 1e-3
+        eta = np.zeros_like(sigma)
+        eta[blk_ind] = 0.1
+        sigmaInf = sigma.copy()
+        sigma0 = sigma*(1-eta)
+
+        nElecs = 11
+        x_temp = np.linspace(-100, 100, nElecs)
+        aSpacing = x_temp[1]-x_temp[0]
+        y_temp = 0.
+        xyz = Utils.ndgrid(x_temp, np.r_[y_temp], np.r_[0.])
+        srcList = DC.Utils.WennerSrcList(nElecs,aSpacing)
+        survey = DC.SurveyDC(srcList)
+
+        imap   = Maps.IdentityMap(mesh)
+        problem = DC.ProblemDC_CC(mesh, mapping=imap)
+
+        try:
+            from pymatsolver import MumpsSolver
+            solver = MumpsSolver
+        except ImportError, e:
+            solver = SolverLU
+
+        problem.Solver = solver
+        problem.pair(survey)
+
+        phi0 = survey.dpred(sigma0)
+        phiInf = survey.dpred(sigmaInf)
+
+        phiIP_true = phi0-phiInf
+
+        surveyIP = DC.SurveyIP(srcList)
+        problemIP = DC.ProblemIP(mesh, sigma=sigma)
+        problemIP.pair(surveyIP)
+
+        problemIP.Solver = solver
+
+        phiIP_approx = surveyIP.dpred(eta)
+
+        err =  np.linalg.norm(phiIP_true-phiIP_approx) / np.linalg.norm(phiIP_true)
+
+        self.assertTrue(err < 0.02)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/dcip/test_sens_IPproblem.py

@@ -0,0 +1,82 @@
+import unittest
+from SimPEG import *
+import SimPEG.DCIP as DC
+
+class IPProblemTests(unittest.TestCase):
+
+    def setUp(self):
+
+        cs = 12.5
+        nc = 500/cs+1
+        hx = [(cs,0, -1.3),(cs,nc),(cs,0, 1.3)]
+        hy = [(cs,0, -1.3),(cs,int(nc/2+1)),(cs,0, 1.3)]
+        hz = [(cs,0, -1.3),(cs,int(nc/2+1))]
+        mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
+        sighalf = 1e-2
+        sigma = np.ones(mesh.nC)*sighalf
+        p0 = np.r_[-50., 50., -50.]
+        p1 = np.r_[ 50.,-50., -150.]
+        blk_ind = Utils.ModelBuilder.getIndicesBlock(p0, p1, mesh.gridCC)
+        sigma[blk_ind] = 1e-3
+        eta = np.zeros_like(sigma)
+        eta[blk_ind] = 0.1
+
+        nElecs = 5
+        x_temp = np.linspace(-250, 250, nElecs)
+        aSpacing = x_temp[1]-x_temp[0]
+        y_temp = 0.
+        xyz = Utils.ndgrid(x_temp, np.r_[y_temp], np.r_[0.])
+        srcList = DC.Utils.WennerSrcList(nElecs,aSpacing)
+        survey = DC.SurveyIP(srcList)
+        imap   = Maps.IdentityMap(mesh)
+        problem = DC.ProblemIP(mesh, sigma=sigma, mapping= imap)
+        problem.pair(survey)
+
+        try:
+            from pymatsolver import MumpsSolver
+            problem.Solver = MumpsSolver
+        except ImportError, e:
+            problem.Solver = SolverLU
+
+        mSynth = eta
+        survey.makeSyntheticData(mSynth)
+
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p = problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0*0, plotIt=False)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-10
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/mesh/test_operators.py

@@ -3,7 +3,7 @@ import unittest
 from SimPEG.Tests import OrderTest
 import matplotlib.pyplot as plt
 
-#TODO: 'randomTensorMesh'
+# TODO: 'randomTensorMesh'
 MESHTYPES = ['uniformTensorMesh', 'uniformCurv', 'rotateCurv']
 call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
 call3 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1], xyz[:, 2])