Updates to documentation and DC utilities.

Conflicts:
	SimPEG/DCIP/BaseDC.py
	SimPEG/DCIP/BaseIP.py

.bumpversion.cfg

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

.gitignore

@@ -40,4 +40,3 @@ 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/dcip
+    - TEST_DIR=tests/em/static
     - TEST_DIR=tests/flow
     - TEST_DIR=tests/mt
     - TEST_DIR=tests/examples
@@ -58,7 +58,8 @@ script:
 
 # Calculate coverage
 after_success:
-  - bash <(curl -s https://codecov.io/bash)
+  - coveralls --config_file .coveragerc
+
   - if [ "$TRAVIS_BRANCH" = "master" -a "$TRAVIS_PULL_REQUEST" = "false" ]; then
       if [ ${TEST_DIR} == "tests/docs" ]; then
         python scripts/fetch_gae_sdk.py $(dirname "${GAE_PYTHONPATH}");

README.rst

@@ -1,4 +1,4 @@
-.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png
+.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
     :alt: SimPEG Logo
 
 ======
@@ -28,11 +28,7 @@ SimPEG
 .. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
     :alt: gitter chat room at https://gitter.im/simpeg/simpeg
     :target: https://gitter.im/simpeg/simpeg
-
-.. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg
-  :target: https://codecov.io/gh/simpeg/simpeg
-
-
+    
 Simulation and Parameter Estimation in Geophysics  -  A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
 
 The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:

SimPEG/DCIP/BaseDC.py

@@ -1,292 +0,0 @@
-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

@@ -1,182 +0,0 @@
-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

@@ -1,38 +0,0 @@
-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

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

SimPEG/Directives.py

@@ -253,7 +253,8 @@ class SaveOutputDictEveryIteration(SaveEveryIteration):
 class Update_IRLS(InversionDirective):
 
     eps_min = None
-    eps = None
+    eps_p = None
+    eps_q = None
     norms = [2.,2.,2.,2.]
     factor = None
     gamma = None
@@ -262,7 +263,6 @@ 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,22 +297,9 @@ class Update_IRLS(InversionDirective):
             print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
 
             self.mode = 2
-
-            # 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)
-            
+            print self.eps_p, self.eps_q, self.norms
+            self.reg.eps_p = self.eps_p
+            self.reg.eps_q = self.eps_q
             self.reg.norms = self.norms
             self.coolingFactor = 1.
             self.coolingRate = 1
@@ -356,14 +343,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*omega(f)*epsilon
+    sig_hat = sig + 1j*omeg*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=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Ex, Ey, Ez = E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
     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=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Ex_galvanic, Ey_galvanic, Ez_galvanic = E_galvanic_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
     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=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Ex_inductive, Ey_inductive, Ez_inductive = E_inductive_from_ElectricDipoleWholeSpaced(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
     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=current, length=length, orientation=orientation, kappa=kappa, epsr=epsr)
+    Hx, Hy, Hz = H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.)
     Bx = mu*Hx
     By = mu*Hy
     Bz = mu*Hz

SimPEG/EM/Static/Utils/StaticUtils.py

@@ -219,10 +219,12 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
         for ii in range(0, int(nstn)-1):
 
 
-            if stype == 'dpdp':
+            if stype == 'dipole-dipole':
                 tx = np.c_[M[ii,:],N[ii,:]]
-            elif stype == 'pdp':
+            elif stype == 'pole-dipole':
                 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,:]])
 
@@ -256,10 +258,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 == 'dpdp':
-                    srcClass = DC.Src.Dipole([rxClass], M[ii,:],N[ii,:])
-            elif stype == 'pdp':
-                    srcClass = DC.Src.Pole([rxClass], M[ii,:])
+            if stype == 'dipole-dipole':
+                srcClass = DC.Src.Dipole([rxClass], M[ii,:],N[ii,:])
+            elif stype == 'pole-dipole':
+                srcClass = DC.Src.Pole([rxClass], M[ii,:])
             SrcList.append(srcClass)
 
     elif stype == 'gradient':
@@ -310,7 +312,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 'pdp', 'dpdp' or 'gradient'. """
+        print """stype must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """
 
 
     return SrcList

SimPEG/EM/Static/__init__.py

@@ -1,3 +1,4 @@
 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=True):
+def run(plotIt=False):
     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=True):
 
 
 if __name__ == '__main__':
-    print run()
+    print run(plotIt=True)

SimPEG/Examples/DC_Forward_PseudoSection.py

@@ -1,5 +1,5 @@
 from SimPEG import Mesh, Utils, np, sp
-import SimPEG.DCIP as DC
+from SimPEG.EM.Static import DC, Utils as StaticUtils
 import time
 
 def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', unitType='appConductivity', plotIt=True):
@@ -74,7 +74,13 @@ 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])
-    survey, 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
 
     # Define some global geometry
     dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )

SimPEG/Examples/Inversion_IRLS.py

@@ -42,33 +42,55 @@ 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()
-    
-    # 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)
+    IRLS = Directives.Update_IRLS( norms=norms,  eps_p=eps_p, eps_q=eps_q)
 
     inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
 

SimPEG/Examples/Maps_ComboMaps.py

@@ -1,62 +0,0 @@
-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

@@ -1,41 +0,0 @@
-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=True, nx=5, ny=5):
+def run(plotIt=False, nx=5, ny=5):
     """
 
         Utils: surface2ind_topo

SimPEG/Examples/__init__.py

@@ -10,8 +10,6 @@ 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
@@ -24,7 +22,7 @@ import MT_1D_ForwardAndInversion
 import MT_3D_Foward
 import Utils_surface2ind_topo
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "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"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Maps.py

@@ -502,9 +502,7 @@ class InjectActiveCells(IdentityMap):
         if Utils.isScalar(valInactive):
             self.valInactive = np.ones(self.nC)*float(valInactive)
         else:
-            self.valInactive = np.ones(self.nC)
-            self.valInactive[self.indInactive] = valInactive.copy()
-            
+            self.valInactive = valInactive.copy()
         self.valInactive[self.indActive] = 0
 
         inds = np.nonzero(self.indActive)[0]

SimPEG/__init__.py

@@ -15,7 +15,7 @@ import Directives
 import Inversion
 import Tests
 
-__version__   = '0.1.12'
+__version__   = '0.1.11'
 __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.12'
+version = '0.1.11'
 # The full version, including alpha/beta/rc tags.
-release = '0.1.12'
+release = '0.1.11'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

docs/content/api_core/api_Maps.rst

@@ -63,8 +63,26 @@ done by the :class:`SimPEG.Maps.ExpMap` described above.
 
 .. plot::
 
-    from SimPEG import Examples
-    Examples.Maps_ComboMaps.run()
+    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()
 
 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).
@@ -149,10 +167,31 @@ Map 2D Cross-Section to 3D Model
 Mesh to Mesh Map
 ----------------
 
+
 .. plot::
 
-    from SimPEG import Examples
-    Examples.Maps_Mesh2Mesh.run()
+    from SimPEG import *
+    import matplotlib.pyplot as plt
+    M = Mesh.TensorMesh([100,100])
+    h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
+    h1 = h1/h1.sum()
+    M2 = Mesh.TensorMesh([h1,h1])
+    V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
+    v = Utils.mkvc(V)
+    modh = Maps.Mesh2Mesh([M,M2])
+    modH = Maps.Mesh2Mesh([M2,M])
+    H = modH * v
+    h = modh * H
+    ax = plt.subplot(131)
+    M.plotImage(v, ax=ax)
+    ax.set_title('Fine Mesh (Original)')
+    ax = plt.subplot(132)
+    M2.plotImage(H,clim=[0,1],ax=ax)
+    ax.set_title('Course Mesh')
+    ax = plt.subplot(133)
+    M.plotImage(h,clim=[0,1],ax=ax)
+    ax.set_title('Fine Mesh (Interpolated)')
+    plt.show()
 
 
 .. autoclass:: SimPEG.Maps.Mesh2Mesh

docs/content/dc/index.rst

@@ -40,7 +40,7 @@
 Direct Current Resistivity
 **************************
 
-`SimPEG.DCIP` uses SimPEG as the framework for the forward and inverse
+`SimPEG.EM.Static.DC` and `SimPEG.EM.Static.IP` 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.DCIP.BaseDC
+.. automodule:: SimPEG.EM.Static.DC
     :show-inheritance:
     :members:
     :undoc-members:

docs/content/examples/Inversion_IRLS.rst

@@ -1,26 +0,0 @@
-.. _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

@@ -1,48 +0,0 @@
-.. _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

@@ -1,27 +0,0 @@
-.. _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.DCIP.BaseIP
+.. automodule:: SimPEG.EM.Static.IP
     :show-inheritance:
     :members:
     :undoc-members:

docs/index.rst

@@ -1,4 +1,4 @@
-.. image:: https://raw.github.com/simpeg/simpeg/master/docs/images/simpeg-logo.png
+.. image:: https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png
     :alt: SimPEG Logo
 
 SimPEG Documentation
@@ -24,10 +24,6 @@ SimPEG Documentation
     :target: https://coveralls.io/r/simpeg/simpeg?branch=master
     :alt: Coverage status
 
-.. image:: https://codecov.io/gh/simpeg/simpeg/branch/master/graph/badge.svg
-    :target: https://codecov.io/gh/simpeg/simpeg
-
-
 Simulation and Parameter Estimation in Geophysics  -  A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
 
 Our vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:

setup.py

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

tests/dcip/__init__.py

@@ -1,12 +0,0 @@
-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

@@ -1,77 +0,0 @@
-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

@@ -1,65 +0,0 @@
-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

@@ -1,82 +0,0 @@
-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()