Merge pull request #362 from simpeg/patch/coveragedocs

Patch/coveragedocs

.bumpversion.cfg

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

.gitignore

@@ -39,3 +39,5 @@ nosetests.xml
 *.sublime-workspace
 docs/_build/
 Makefile
+docs/warnings.txt
+.DS_Store

.travis.yml

@@ -24,18 +24,25 @@ env:
     - TEST_DIR=tests/examples
     - TEST_DIR=tests/em/fdem/inverse/adjoint
     - TEST_DIR=tests/em/fdem/forward
+    - TEST_DIR=tests/docs;
+      GAE_PYTHONPATH=${HOME}/.cache/google_appengine;
+      PATH=$PATH:${HOME}/google-cloud-sdk/bin;
+      PYTHONPATH=${PYTHONPATH}:${GAE_PYTHONPATH};
+      CLOUDSDK_CORE_DISABLE_PROMPTS=1
 
 # Setup anaconda
 before_install:
-  - if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh -O miniconda.sh; fi
+# Install packages
+  - if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.8.3-Linux-x86_64.sh
+    -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.8.3-Linux-x86_64.sh
+    -O miniconda.sh; fi
   - chmod +x miniconda.sh
   - ./miniconda.sh -b
   - export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH
   - conda update --yes conda
 
-# Install packages
 install:
-  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython ipywidgets nose vtk
+  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk sphinx
   - pip install nose-cov python-coveralls
 
   - git clone https://github.com/rowanc1/pymatsolver.git
@@ -46,11 +53,26 @@ install:
 
 # Run test
 script:
+  # test docs
   - nosetests $TEST_DIR --with-cov --cov SimPEG --cov-config .coveragerc -v -s
 
 # Calculate coverage
 after_success:
-  - coveralls --config_file .coveragerc
+  - bash <(curl -s https://codecov.io/bash)
+  - 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}");
+        openssl aes-256-cbc -K $encrypted_93066031461c_key -iv $encrypted_93066031461c_iv
+          -in docs/credentials.tar.gz.enc -out credentials.tar.gz -d ;
+        if [ ! -d ${HOME}/google-cloud-sdk ]; then curl https://sdk.cloud.google.com | bash; fi ;
+        tar -xzf credentials.tar.gz ;
+        gcloud auth activate-service-account --key-file client-secret.json ;
+        gcloud config set project simpegdocs;
+        gcloud -q components update gae-python;
+        gcloud -q preview app deploy ./docs/app.yaml --version ${TRAVIS_COMMIT} --promote;
+      fi;
+    fi
+
 
 notifications:
   email:

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
 
 ======
@@ -15,7 +15,7 @@ SimPEG
 
 .. image:: https://img.shields.io/badge/license-MIT-blue.svg
     :target: https://github.com/simpeg/simpeg/blob/master/LICENSE
-    :alt: BSD 3 clause license.
+    :alt: MIT license
 
 .. image:: https://api.travis-ci.org/simpeg/simpeg.svg?branch=master
     :target: https://travis-ci.org/simpeg/simpeg
@@ -28,7 +28,12 @@ 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
+    :alt: Coverage status
+
+
 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

@@ -162,8 +162,8 @@ class ProblemDC_CC(Problem.BaseProblem):
         """
             Makes the matrix A(m) for the DC resistivity problem.
 
-            :param numpy.array m: model
-            :rtype: scipy.csc_matrix
+            :param numpy.ndarray m: model
+            :rtype: scipy.sparse.csc_matrix
             :return: A(m)
 
             .. math::

SimPEG/DCIP/BaseIP.py

@@ -71,7 +71,7 @@ class ProblemIP(Problem.BaseProblem):
             Makes the matrix A(m) for the DC resistivity problem.
 
             :param numpy.array m: model
-            :rtype: scipy.csc_matrix
+            :rtype: scipy.sparse.csc_matrix
             :return: A(m)
 
             .. math::

SimPEG/Directives.py

@@ -167,7 +167,7 @@ class TargetMisfit(InversionDirective):
 
 
 
-class _SaveEveryIteration(InversionDirective):
+class SaveEveryIteration(InversionDirective):
     @property
     def name(self):
         if getattr(self, '_name', None) is None:
@@ -188,7 +188,7 @@ class _SaveEveryIteration(InversionDirective):
         self._fileName = value
 
 
-class SaveModelEveryIteration(_SaveEveryIteration):
+class SaveModelEveryIteration(SaveEveryIteration):
     """SaveModelEveryIteration"""
 
     def initialize(self):
@@ -198,7 +198,7 @@ class SaveModelEveryIteration(_SaveEveryIteration):
         np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
 
 
-class SaveOutputEveryIteration(_SaveEveryIteration):
+class SaveOutputEveryIteration(SaveEveryIteration):
     """SaveModelEveryIteration"""
 
     def initialize(self):
@@ -212,40 +212,48 @@ class SaveOutputEveryIteration(_SaveEveryIteration):
         f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
         f.close()
 
-class SaveOutputDictEveryIteration(_SaveEveryIteration):
-    """
-    Saves inversion parameters at every iteraion.
-
-
-    """
+class SaveOutputDictEveryIteration(SaveEveryIteration):
+    """SaveOutputDictEveryIteration"""
 
     def initialize(self):
         print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
 
     def endIter(self):
-
-        # Initialize the output dict
-        outDict = {}
         # Save the data.
-        outDict['iter'] = self.opt.iter
-        outDict['beta'] = self.invProb.beta
-        outDict['phi_d'] = self.invProb.phi_d
-        outDict['phi_ms'] = self.reg._evalSmall(self.invProb.curModel)
-        outDict['phi_mx'] = self.reg._evalSmoothx(self.invProb.curModel)
-        outDict['phi_my'] = self.reg._evalSmoothy(self.invProb.curModel) if self.prob.mesh.dim >= 2 else 'NaN'
-        outDict['phi_mz'] = self.reg._evalSmoothz(self.invProb.curModel) if self.prob.mesh.dim==3 else 'NaN'
-        outDict['f'] = self.opt.f
-        outDict['m'] = self.invProb.curModel
-        outDict['dpred'] = self.invProb.dpred
+        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
+        phi_ms = 0.5*ms.dot(ms)
+        if self.reg.mrefInSmooth == True:
+            mref = self.reg.mref
+        else:
+            mref = 0
+        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
+        phi_mx = 0.5 * mx.dot(mx)
+        if self.prob.mesh.dim >= 2:
+            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_my = 0.5 * my.dot(my)
+        else:
+            phi_my = 'NaN'
+        if self.prob.mesh.dim==3:
+            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_mz = 0.5 * mz.dot(mz)
+        else:
+            phi_mz = 'NaN'
+
 
         # Save the file as a npz
-        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), outDict)
+        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
+
+#         mref = getattr(self, 'm_prev', None)
+#         if mref is None:
+#             if self.debug: print 'UpdateReferenceModel is using mref0'
+#             mref = self.mref0
+#         self.m_prev = self.invProb.m_current
+#         return mref
 
 class Update_IRLS(InversionDirective):
 
     eps_min = None
-    eps_p = None
-    eps_q = None
+    eps = None
     norms = [2.,2.,2.,2.]
     factor = None
     gamma = None
@@ -254,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
@@ -288,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
@@ -334,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

@@ -0,0 +1,302 @@
+from __future__ import division
+import numpy as np
+from scipy.constants import mu_0, pi, epsilon_0
+from scipy.special import erf
+from SimPEG import Utils
+
+omega = lambda f: 2.*np.pi*f
+# TODO:
+# r = lambda dx, dy, dz: np.sqrt( dx**2. + dy**2. + dz**2.)
+# k = lambda f, mu, epsilon, sig: np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+def E_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=0., epsr=1.):
+
+    """
+        Computing Analytic Electric fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    sig_hat = sig + 1j*omega(f)*epsilon
+
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
+    mid   = -k**2 * r**2 + 3*1j*k*r + 3
+
+    if orientation.upper() == 'X':
+        Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ey = front*(dx*dy  / r**2)*mid
+        Ez = front*(dx*dz  / r**2)*mid
+        return Ex, Ey, Ez
+
+    elif orientation.upper() == 'Y':
+        #  x--> y, y--> z, z-->x
+        Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ez = front*(dy*dz  / r**2)*mid
+        Ex = front*(dy*dx  / r**2)*mid
+        return Ex, Ey, Ez
+
+    elif orientation.upper() == 'Z':
+        # x --> z, y --> x, z --> y
+        Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
+        Ex = front*(dz*dx  / r**2)*mid
+        Ey = front*(dz*dy  / r**2)*mid
+        return Ex, Ey, Ez
+
+
+def E_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Galvanic portion of Electric fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    sig_hat = sig + 1j*omega(f)*epsilon
+
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
+    mid   = -k**2 * r**2 + 3*1j*k*r + 3
+
+    if orientation.upper() == 'X':
+        Ex_galvanic = front*((dx**2 / r**2)*mid + (-1j*k*r-1.))
+        Ey_galvanic = front*(dx*dy  / r**2)*mid
+        Ez_galvanic = front*(dx*dz  / r**2)*mid
+        return Ex_galvanic, Ey_galvanic, Ez_galvanic
+
+    elif orientation.upper() == 'Y':
+        #  x--> y, y--> z, z-->x
+        Ey_galvanic = front*((dy**2 / r**2)*mid + (-1j*k*r-1.))
+        Ez_galvanic = front*(dy*dz  / r**2)*mid
+        Ex_galvanic = front*(dy*dx  / r**2)*mid
+        return Ex_galvanic, Ey_galvanic, Ez_galvanic
+
+    elif orientation.upper() == 'Z':
+        # x --> z, y --> x, z --> y
+        Ez_galvanic = front*((dz**2 / r**2)*mid + (-1j*k*r-1.))
+        Ex_galvanic = front*(dz*dx  / r**2)*mid
+        Ey_galvanic = front*(dz*dy  / r**2)*mid
+        return Ex_galvanic, Ey_galvanic, Ez_galvanic
+
+
+def E_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Inductive portion of Electric fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    sig_hat = sig + 1j*omega(f)*epsilon
+
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*sig_hat* r**3) * np.exp(-1j*k*r)
+
+    if orientation.upper() == 'X':
+        Ex_inductive = front*(k**2 * r**2)
+        Ey_inductive = np.zeros_like(Ex_inductive)
+        Ez_inductive = np.zeros_like(Ex_inductive)
+        return Ex_inductive, Ey_inductive, Ez_inductive
+
+    elif orientation.upper() == 'Y':
+        #  x--> y, y--> z, z-->x
+        Ey_inductive = front*(k**2 * r**2)
+        Ez_inductive = np.zeros_like(Ey_inductive)
+        Ex_inductive = np.zeros_like(Ey_inductive)
+        return Ex_inductive, Ey_inductive, Ez_inductive
+
+    elif orientation.upper() == 'Z':
+        # x --> z, y --> x, z --> y
+        Ez_inductive = front*(k**2 * r**2)
+        Ex_inductive = np.zeros_like(Ez_inductive)
+        Ey_inductive = np.zeros_like(Ez_inductive)
+        return Ex_inductive, Ey_inductive, Ez_inductive
+
+
+def J_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Current densities from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+
+    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
+    return Jx, Jy, Jz
+
+
+def J_galvanic_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Galvanic portion of Current densities from Electrical Dipole in a Wholespace
+        TODO:
+            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)
+    Jx_galvanic = sig*Ex_galvanic
+    Jy_galvanic = sig*Ey_galvanic
+    Jz_galvanic = sig*Ez_galvanic
+    return Jx_galvanic, Jy_galvanic, Jz_galvanic
+
+
+def J_inductive_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Inductive portion of Current densities from Electrical Dipole in a Wholespace
+        TODO:
+            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)
+    Jx_inductive = sig*Ex_inductive
+    Jy_inductive = sig*Ey_inductive
+    Jz_inductive = sig*Ez_inductive
+    return Jx_inductive, Jy_inductive, Jz_inductive
+
+
+def H_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Magnetic fields from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    # k  = np.sqrt( -1j*2.*np.pi*f*mu*sig )
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi* r**2) * (-1j*k*r + 1) * np.exp(-1j*k*r)
+
+    if orientation.upper() == 'X':
+        Hy = front*(-dz  / r)
+        Hz = front*(dy  / r)
+        Hx = np.zeros_like(Hy)
+        return Hx, Hy, Hz
+
+    elif orientation.upper() == 'Y':
+        Hx = front*(dz  / r)
+        Hz = front*(-dx  / r)
+        Hy = np.zeros_like(Hx)
+        return Hx, Hy, Hz
+
+    elif orientation.upper() == 'Z':
+        Hx = front*(-dy / r)
+        Hy = front*(dx  / r)
+        Hz = np.zeros_like(Hx)
+        return Hx, Hy, Hz
+
+
+def B_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Magnetic flux densites from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+
+    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
+    return Bx, By, Bz
+
+
+def A_from_ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', kappa=1., epsr=1.):
+
+    """
+        Computing Electric vector potentials from Electrical Dipole in a Wholespace
+        TODO:
+            Add description of parameters
+    """
+    mu = mu_0*(1+kappa)
+    epsilon = epsilon_0*epsr
+    XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
+    # Check
+    if XYZ.shape[0] > 1 & f.shape[0] > 1:
+        raise Exception("I/O type error: For multiple field locations only a single frequency can be specified.")
+
+    dx = XYZ[:,0]-srcLoc[0]
+    dy = XYZ[:,1]-srcLoc[1]
+    dz = XYZ[:,2]-srcLoc[2]
+
+    r  = np.sqrt( dx**2. + dy**2. + dz**2.)
+    k  = np.sqrt( omega(f)**2. *mu*epsilon -1j*omega(f)*mu*sig )
+
+    front = current * length / (4.*np.pi*r)
+
+    if orientation.upper() == 'X':
+        Ax = front*np.exp(-1j*k*r)
+        Ay = np.zeros_like(Ax)
+        Az = np.zeros_like(Ax)
+        return Ax, Ay, Az
+
+    elif orientation.upper() == 'Y':
+        Ay = front*np.exp(-1j*k*r)
+        Ax = np.zeros_like(Ay)
+        Az = np.zeros_like(Ay)
+        return Ax, Ay, Az
+
+    elif orientation.upper() == 'Z':
+        Az = front*np.exp(-1j*k*r)
+        Ax = np.zeros_like(Ay)
+        Ay = np.zeros_like(Ay)
+        return Ax, Ay, Az
+
+
+
+
+

SimPEG/EM/Analytics/__init__.py

@@ -2,3 +2,4 @@ from TDEM import hzAnalyticDipoleT
 from FDEM import hzAnalyticDipoleF
 from FDEMcasing import *
 from DC import DCAnalyticHalf, DCAnalyticSphere
+from FDEMDipolarfields import *

SimPEG/EM/Base.py

@@ -20,10 +20,10 @@ class BaseEMProblem(Problem.BaseProblem):
         Problem.BaseProblem.__init__(self, mesh, **kwargs)
 
 
-    surveyPair = Survey.BaseSurvey
-    dataPair = Survey.Data
+    surveyPair = Survey.BaseSurvey #: The survey to pair with.
+    dataPair = Survey.Data #: The data to pair with.
 
-    PropMap = EMPropMap
+    PropMap = EMPropMap #: The property mapping
 
     Solver = SimpegSolver
     solverOpts = {}
@@ -217,6 +217,7 @@ class BaseEMSurvey(Survey.BaseSurvey):
     def eval(self, f):
         """
         Project fields to receiver locations
+
         :param Fields u: fields object
         :rtype: numpy.ndarray
         :return: data

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -6,11 +6,11 @@ from SimPEG.EM.Utils import omega
 from SimPEG.Utils import Zero, Identity, sdiag
 
 
-class Fields(SimPEG.Problem.Fields):
+class FieldsFDEM(SimPEG.Problem.Fields):
     """
 
     Fancy Field Storage for a FDEM survey. Only one field type is stored for
-    each problem, the rest are computed. The fields obejct acts like an array and is indexed by
+    each problem, the rest are computed. The fields object acts like an array and is indexed by
 
     .. code-block:: python
 
@@ -92,7 +92,7 @@ class Fields(SimPEG.Problem.Fields):
         """
         Total derivative of e with respect to the inversion model. Returns :math:`d\mathbf{e}/d\mathbf{m}` for forward and (:math:`d\mathbf{e}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
 
-        :param Src src: sorce
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
         :param numpy.ndarray v: vector to take sensitivity product with
         :param bool adjoint: adjoint?
@@ -110,7 +110,7 @@ class Fields(SimPEG.Problem.Fields):
         """
         Total derivative of b with respect to the inversion model. Returns :math:`d\mathbf{b}/d\mathbf{m}` for forward and (:math:`d\mathbf{b}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
 
-        :param Src src: sorce
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
         :param numpy.ndarray v: vector to take sensitivity product with
         :param bool adjoint: adjoint?
@@ -128,7 +128,7 @@ class Fields(SimPEG.Problem.Fields):
         """
         Total derivative of h with respect to the inversion model. Returns :math:`d\mathbf{h}/d\mathbf{m}` for forward and (:math:`d\mathbf{h}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
 
-        :param Src src: sorce
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
         :param numpy.ndarray v: vector to take sensitivity product with
         :param bool adjoint: adjoint?
@@ -146,7 +146,7 @@ class Fields(SimPEG.Problem.Fields):
         """
         Total derivative of j with respect to the inversion model. Returns :math:`d\mathbf{j}/d\mathbf{m}` for forward and (:math:`d\mathbf{j}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
 
-        :param Src src: sorce
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: source
         :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
         :param numpy.ndarray v: vector to take sensitivity product with
         :param bool adjoint: adjoint?
@@ -160,12 +160,12 @@ class Fields(SimPEG.Problem.Fields):
             return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
         return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex)
 
-class Fields3D_e(Fields):
+class Fields3D_e(FieldsFDEM):
     """
     Fields object for Problem3D_e.
 
-    :param Mesh mesh: mesh
-    :param Survey survey: survey
+    :param BaseMesh mesh: mesh
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'eSolution':'E'}
@@ -180,9 +180,6 @@ class Fields3D_e(Fields):
                     'h' : ['eSolution','CCV','_h'],
                   }
 
-    def __init__(self, mesh, survey, **kwargs):
-        Fields.__init__(self, mesh, survey, **kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -426,12 +423,12 @@ class Fields3D_e(Fields):
 
 
 
-class Fields3D_b(Fields):
+class Fields3D_b(FieldsFDEM):
     """
     Fields object for Problem3D_b.
 
-    :param Mesh mesh: mesh
-    :param Survey survey: survey
+    :param BaseMesh mesh: mesh
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'bSolution':'F'}
@@ -446,9 +443,6 @@ class Fields3D_b(Fields):
                     'h' : ['bSolution','CCV','_h'],
                   }
 
-    def __init__(self,mesh,survey,**kwargs):
-        Fields.__init__(self,mesh,survey,**kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -693,12 +687,12 @@ class Fields3D_b(Fields):
         return Zero()
 
 
-class Fields3D_j(Fields):
+class Fields3D_j(FieldsFDEM):
     """
     Fields object for Problem3D_j.
 
-    :param Mesh mesh: mesh
-    :param Survey survey: survey
+    :param BaseMesh mesh: mesh
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'jSolution':'F'}
@@ -713,9 +707,6 @@ class Fields3D_j(Fields):
                     'b' : ['jSolution','CCV','_b'],
                   }
 
-    def __init__(self,mesh,survey,**kwargs):
-        Fields.__init__(self,mesh,survey,**kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl
@@ -988,12 +979,12 @@ class Fields3D_j(Fields):
         return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
 
 
-class Fields3D_h(Fields):
+class Fields3D_h(FieldsFDEM):
     """
     Fields object for Problem3D_h.
 
-    :param Mesh mesh: mesh
-    :param Survey survey: survey
+    :param BaseMesh mesh: mesh
+    :param SimPEG.EM.FDEM.SurveyFDEM.Survey survey: survey
     """
 
     knownFields = {'hSolution':'E'}
@@ -1008,9 +999,6 @@ class Fields3D_h(Fields):
                     'b' : ['hSolution','CCV','_b'],
                   }
 
-    def __init__(self,mesh,survey,**kwargs):
-        Fields.__init__(self,mesh,survey,**kwargs)
-
     def startup(self):
         self.prob = self.survey.prob
         self._edgeCurl = self.survey.prob.mesh.edgeCurl

SimPEG/EM/FDEM/ProblemFDEM.py

@@ -1,7 +1,7 @@
 from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
 from SurveyFDEM import Survey as SurveyFDEM
-from FieldsFDEM import Fields, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
+from FieldsFDEM import FieldsFDEM, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
 from SimPEG.EM.Base import BaseEMProblem
 from SimPEG.EM.Utils import omega
 
@@ -31,10 +31,11 @@ class BaseFDEMProblem(BaseEMProblem):
         if using the H-J formulation (:code:`Problem3D_j` or :code:`Problem3D_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
 
         The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
+
     """
 
     surveyPair = SurveyFDEM
-    fieldsPair = Fields
+    fieldsPair = FieldsFDEM
 
     def fields(self, m):
         """
@@ -64,7 +65,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         :param numpy.array m: inversion model (nP,)
         :param numpy.array v: vector which we take sensitivity product with (nP,)
-        :param SimPEG.EM.FDEM.Fields u: fields object
+        :param SimPEG.EM.FDEM.FieldsFDEM.FieldsFDEM u: fields object
         :rtype numpy.array:
         :return: Jv (ndata,)
         """
@@ -99,7 +100,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         :param numpy.array m: inversion model (nP,)
         :param numpy.array v: vector which we take adjoint product with (nP,)
-        :param SimPEG.EM.FDEM.Fields u: fields object
+        :param SimPEG.EM.FDEM.FieldsFDEM.FieldsFDEM u: fields object
         :rtype numpy.array:
         :return: Jv (ndata,)
         """
@@ -153,8 +154,8 @@ class BaseFDEMProblem(BaseEMProblem):
         Evaluates the sources for a given frequency and puts them in matrix form
 
         :param float freq: Frequency
-        :rtype: (numpy.ndarray, numpy.ndarray)
-        :return: s_m, s_e (nE or nF, nSrc)
+        :rtype: tuple
+        :return: (s_m, s_e) (nE or nF, nSrc)
         """
         Srcs = self.survey.getSrcByFreq(freq)
         if self._formulation is 'EB':
@@ -194,7 +195,7 @@ class Problem3D_e(BaseFDEMProblem):
 
     which we solve for :math:`\mathbf{e}`.
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'eSolution'
@@ -269,7 +270,7 @@ class Problem3D_e(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -305,7 +306,7 @@ class Problem3D_b(BaseFDEMProblem):
     .. note ::
         The inverse problem will not work with full anisotropy
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'bSolution'
@@ -400,7 +401,7 @@ class Problem3D_b(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -444,6 +445,7 @@ class Problem3D_j(BaseFDEMProblem):
 
         \mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j}  + \mathbf{M^e} \mathbf{s_m} \\right)
 
+
     and solve for \\\(\\\mathbf{j}\\\) using
 
     .. math ::
@@ -453,7 +455,7 @@ class Problem3D_j(BaseFDEMProblem):
     .. note::
         This implementation does not yet work with full anisotropy!!
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'jSolution'
@@ -529,8 +531,8 @@ class Problem3D_j(BaseFDEMProblem):
             \mathbf{RHS} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega \mathbf{s_e}
 
         :param float freq: Frequency
-        :rtype: numpy.ndarray (nE, nSrc)
-        :return: RHS
+        :rtype: numpy.ndarray
+        :return: RHS (nE, nSrc)
         """
 
         s_m, s_e = self.getSourceTerm(freq)
@@ -549,7 +551,7 @@ class Problem3D_j(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -591,7 +593,7 @@ class Problem3D_h(BaseFDEMProblem):
 
         \\left(\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e}
 
-    :param SimPEG.Mesh mesh: mesh
+    :param SimPEG.Mesh.BaseMesh.BaseMesh mesh: mesh
     """
 
     _solutionType = 'hSolution'
@@ -608,9 +610,11 @@ class Problem3D_h(BaseFDEMProblem):
         .. math::
             \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}
 
+
         :param float freq: Frequency
         :rtype: scipy.sparse.csr_matrix
         :return: A
+
         """
 
         MeMu = self.MeMu
@@ -653,6 +657,7 @@ class Problem3D_h(BaseFDEMProblem):
         :param float freq: Frequency
         :rtype: numpy.ndarray
         :return: RHS (nE, nSrc)
+
         """
 
         s_m, s_e = self.getSourceTerm(freq)
@@ -666,7 +671,7 @@ class Problem3D_h(BaseFDEMProblem):
         Derivative of the right hand side with respect to the model
 
         :param float freq: frequency
-        :param SimPEG.EM.FDEM.Src src: FDEM source
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray

SimPEG/EM/FDEM/RxFDEM.py

@@ -25,10 +25,10 @@ class BaseRx(SimPEG.Survey.BaseRx):
 
     def eval(self, src, mesh, f):
         """
-        Project fields to recievers to get data.
+        Project fields to receivers to get data.
 
-        :param Source src: FDEM source
-        :param Mesh mesh: mesh used
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
+        :param BaseMesh mesh: mesh used
         :param Fields f: fields object
         :rtype: numpy.ndarray
         :return: fields projected to recievers
@@ -44,8 +44,8 @@ class BaseRx(SimPEG.Survey.BaseRx):
         """
         Derivative of projected fields with respect to the inversion model times a vector.
 
-        :param Source src: FDEM source
-        :param Mesh mesh: mesh used
+        :param SimPEG.EM.FDEM.SrcFDEM.BaseSrc src: FDEM source
+        :param BaseMesh mesh: mesh used
         :param Fields f: fields object
         :param numpy.ndarray v: vector to multiply
         :rtype: numpy.ndarray

SimPEG/EM/FDEM/SrcFDEM.py

@@ -23,8 +23,8 @@ class BaseSrc(Survey.BaseSrc):
         - :math:`s_m` : magnetic source term
         - :math:`s_e` : electric source term
 
-        :param Problem prob: FDEM Problem
-        :rtype: (numpy.ndarray, numpy.ndarray)
+        :param BaseFDEMProblem prob: FDEM Problem
+        :rtype: tuple
         :return: tuple with magnetic source term and electric source term
         """
         s_m = self.s_m(prob)
@@ -37,10 +37,10 @@ class BaseSrc(Survey.BaseSrc):
         - :code:`s_mDeriv` : derivative of the magnetic source term
         - :code:`s_eDeriv` : derivative of the electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
-        :rtype: (numpy.ndarray, numpy.ndarray)
+        :rtype: tuple
         :return: tuple with magnetic source term and electric source term derivatives times a vector
         """
         if v is not None:
@@ -52,7 +52,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary magnetic flux density
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary magnetic flux density
         """
@@ -64,7 +64,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary magnetic field
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -76,7 +76,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary electric field
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary electric field
         """
@@ -88,7 +88,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Primary current density
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: primary current density
         """
@@ -100,7 +100,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Magnetic source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: magnetic source term on mesh
         """
@@ -110,7 +110,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: electric source term on mesh
         """
@@ -120,7 +120,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Derivative of magnetic source term with respect to the inversion model
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -133,7 +133,7 @@ class BaseSrc(Survey.BaseSrc):
         """
         Derivative of electric source term with respect to the inversion model
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :param numpy.ndarray v: vector to take product with
         :param bool adjoint: adjoint?
         :rtype: numpy.ndarray
@@ -162,7 +162,7 @@ class RawVec_e(BaseSrc):
         """
         Electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: electric source term on mesh
         """
@@ -191,7 +191,7 @@ class RawVec_m(BaseSrc):
         """
         Magnetic source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: magnetic source term on mesh
         """
@@ -220,7 +220,7 @@ class RawVec(BaseSrc):
         """
         Magnetic source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: magnetic source term on mesh
         """
@@ -232,7 +232,7 @@ class RawVec(BaseSrc):
         """
         Electric source term
 
-        :param Problem prob: FDEM Problem
+        :param BaseFDEMProblem prob: FDEM Problem
         :rtype: numpy.ndarray
         :return: electric source term on mesh
         """
@@ -301,7 +301,7 @@ class MagDipole(BaseSrc):
         """
         The primary magnetic flux density from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -339,7 +339,7 @@ class MagDipole(BaseSrc):
         """
         The primary magnetic field from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -350,7 +350,7 @@ class MagDipole(BaseSrc):
         """
         The magnetic source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -364,7 +364,7 @@ class MagDipole(BaseSrc):
         """
         The electric source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -416,7 +416,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The primary magnetic flux density from the analytic solution for magnetic fields from a dipole
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -455,7 +455,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The primary magnetic field from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -466,7 +466,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The magnetic source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -479,7 +479,7 @@ class MagDipole_Bfield(BaseSrc):
         """
         The electric source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -530,7 +530,7 @@ class CircularLoop(BaseSrc):
         """
         The primary magnetic flux density from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -567,7 +567,7 @@ class CircularLoop(BaseSrc):
         """
         The primary magnetic field from a magnetic vector potential
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -578,7 +578,7 @@ class CircularLoop(BaseSrc):
         """
         The magnetic source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
@@ -591,7 +591,7 @@ class CircularLoop(BaseSrc):
         """
         The electric source term
 
-        :param Problem prob: FDEM problem
+        :param BaseFDEMProblem prob: FDEM problem
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """

SimPEG/EM/TDEM/BaseTDEM.py

@@ -112,7 +112,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.ndarray v: vector (model object)
-            :param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m
+            :param FieldsTDEM f: Fields resulting from m
             :rtype: numpy.ndarray
             :return: w (data object)
 
@@ -136,8 +136,8 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
     def Jtvec(self, m, v, f=None):
         """
             :param numpy.array m: Conductivity model
-            :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
+            :param numpy.ndarray v: vector (or a :class:`SimPEG.Survey.Data` object)
+            :param FieldsTDEM u: Fields resulting from m
             :rtype: numpy.ndarray
             :return: w (model object)
 

SimPEG/EM/TDEM/TDEM_b.py

@@ -87,8 +87,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.array vec: vector (like a model)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :param FieldsTDEM u: Fields resulting from m
+            :rtype: FieldsTDEM
             :return: f
 
             Multiply G by a vector
@@ -125,9 +125,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
         """
             :param numpy.array m: Conductivity model
             :param numpy.array vec: vector (like a fields)
-            :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m
-            :rtype: np.ndarray (like a model)
-            :return: p
+            :param FieldsTDEM u: Fields resulting from m
+            :rtype: numpy.ndarray
+            :return: p (like a model)
 
             Multiply G.T by a vector
         """
@@ -153,8 +153,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def solveAh(self, m, p):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM p: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :param FieldsTDEM p: Fields object
+            :rtype: FieldsTDEM
             :return: y
 
             Solve the block-matrix system \\\(\\\hat{A} \\\hat{y} = \\\hat{p}\\\):
@@ -200,8 +200,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def solveAht(self, m, p):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM p: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :param FieldsTDEM p: Fields object
+            :rtype: FieldsTDEM
             :return: y
 
             Solve the block-matrix system \\\(\\\hat{A}^\\\\top \\\hat{y} = \\\hat{p}\\\):
@@ -270,8 +270,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def _AhVec(self, m, vec):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM vec: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :param FieldsTDEM vec: Fields object
+            :rtype: FieldsTDEM
             :return: f
 
             Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where
@@ -315,8 +315,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
     def _AhtVec(self, m, vec):
         """
             :param numpy.array m: Conductivity model
-            :param simpegEM.TDEM.FieldsTDEM vec: Fields object
-            :rtype: simpegEM.TDEM.FieldsTDEM
+            :param FieldsTDEM vec: Fields object
+            :rtype: FieldsTDEM
             :return: f
 
             Multiply the matrix \\\(\\\hat{A}\\\) by a fields vector where

SimPEG/Examples/DC_Analytic_Dipole.py

@@ -1,7 +1,7 @@
 from SimPEG import *
-import SimPEG.DCIP as DC
+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)]
@@ -21,10 +21,10 @@ def run(plotIt=False):
     #     ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
     #     ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
 
-    rx = DC.RxDipole(xyz_rxP, xyz_rxN)
-    src = DC.SrcDipole([rx], [-200, 0, -12.5], [+200, 0, -12.5])
-    survey = DC.SurveyDC([src])
-    problem = DC.ProblemDC_CC(mesh)
+    rx = DC.Rx.Dipole(xyz_rxP, xyz_rxN)
+    src = DC.Src.Dipole([rx], np.r_[-200, 0, -12.5], np.r_[+200, 0, -12.5])
+    survey = DC.Survey([src])
+    problem = DC.Problem3D_CC(mesh)
     problem.pair(survey)
     try:
         from pymatsolver import MumpsSolver
@@ -65,4 +65,4 @@ def run(plotIt=False):
 
 
 if __name__ == '__main__':
-    print run(plotIt=True)
+    print run()

SimPEG/Examples/EM_Schenkel_Morrison_Casing.py

@@ -19,10 +19,13 @@ def run(plotIt=True):
         Morrison Casing Model, and the results are used in a 2016 SEG abstract by
         Yang et al.
 
-        - Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
+        .. code-block:: text
+
+            Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686.
 
 
         The model consists of:
+
         - Air: Conductivity 1e-8 S/m, above z = 0
         - Background: conductivity 1e-2 S/m, below z = 0
         - Casing: conductivity 1e6 S/m
@@ -215,7 +218,7 @@ def run(plotIt=True):
     # ------------ Problem and Survey ---------------
     survey = FDEM.Survey(sg_p + dg_p)
     mapping = [('sigma', Maps.IdentityMap(mesh))]
-    problem = FDEM.Problem3D_h(mesh, mapping=mapping)
+    problem = FDEM.Problem3D_h(mesh, mapping=mapping, Solver=solver)
     problem.pair(survey)
 
     # ------------- Solve ---------------------------

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/MT_1D_ForwardAndInversion.py

@@ -1,15 +1,13 @@
 import SimPEG as simpeg
 import numpy as np
-from SimPEG import NSEM
+import SimPEG.MT as MT
 from scipy.constants import mu_0
 import matplotlib.pyplot as plt
 
-np.random.seed(1983)
-
 def run(plotIt=True):
     """
         MT: 1D: Inversion
-        =======================
+        =================
 
         Forward model 1D MT data.
         Setup and run a MT 1D inversion.
@@ -19,13 +17,13 @@ def run(plotIt=True):
     ## Setup the forward modeling
     # Setting up 1D mesh and conductivity models to forward model data.
     # Frequency
-    nFreq = 26
-    freqs = np.logspace(2,-3,nFreq)
+    nFreq = 31
+    freqs = np.logspace(3,-3,nFreq)
     # Set mesh parameters
-    ct = 10
-    air = simpeg.Utils.meshTensor([(ct,25,1.4)])
+    ct = 20
+    air = simpeg.Utils.meshTensor([(ct,16,1.4)])
     core = np.concatenate( (  np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )
-    bot = simpeg.Utils.meshTensor([(core[0],25,-1.4)])
+    bot = simpeg.Utils.meshTensor([(core[0],10,-1.4)])
     x0 = -np.array([np.sum(np.concatenate((core,bot)))])
     # Make the model
     m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
@@ -35,7 +33,7 @@ def run(plotIt=True):
     layer1 = (m1d.vectorCCx<-500.) & (m1d.vectorCCx>=-800.)
     layer2 = (m1d.vectorCCx<-3500.) & (m1d.vectorCCx>=-5000.)
     # Set the conductivity values
-    sig_half = 1e-2
+    sig_half = 2e-3
     sig_air = 1e-8
     sig_layer1 = .2
     sig_layer2 = .2
@@ -59,31 +57,31 @@ def run(plotIt=True):
     # Receivers
     rxList = []
     for rxType in ['z1dr','z1di']:
-        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([-0.5]),2).T,rxType))
+        rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
     # Source list
     srcList =[]
     for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
+            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
     # Make the survey
-    survey = NSEM.Survey(srcList)
+    survey = MT.Survey(srcList)
     survey.mtrue = m_true
 
     ## Set the problem
-    problem = NSEM.Problem1D_ePrimSec(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
+    problem = MT.Problem1D.eForm_psField(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
     problem.pair(survey)
 
     ## Forward model data
     # Project the data
     survey.dtrue = survey.dpred(m_true)
-    survey.dobs = survey.dtrue + 0.01*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
+    survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
 
     if plotIt:
-        fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[])
+        fig = MT.Utils.dataUtils.plotMT1DModelData(problem, [m_0])
         fig.suptitle('Target - smooth true')
 
 
     # Assign uncertainties
-    std = 0.025 # 5% std
+    std = 0.05 # 5% std
     survey.std = np.abs(survey.dobs*std)
     # Assign the data weight
     Wd = 1./survey.std
@@ -92,33 +90,30 @@ def run(plotIt=True):
     # Define a counter
     C =  simpeg.Utils.Counter()
     # Set the optimization
-    opt = simpeg.Optimization.ProjectedGNCG(maxIter = 25)
+    opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)
     opt.counter = C
-    opt.lower = np.log(1e-4)
-    opt.upper = np.log(5)
-    opt.LSshorten = 0.1
+    opt.LSshorten = 0.5
     opt.remember('xc')
     # Data misfit
     dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
     dmis.Wd = Wd
     # Regularization - with a regularization mesh
-    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[active]],m1d.x0)
+    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
     reg = simpeg.Regularization.Tikhonov(regMesh)
     reg.mrefInSmooth = True
-    reg.alpha_s = 1e-1
+    reg.alpha_s = 1e-7
     reg.alpha_x = 1.
-
     # Inversion problem
     invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)
     invProb.counter = C
     # Beta cooling
     beta = simpeg.Directives.BetaSchedule()
-    beta.coolingRate = 4.
-    beta.coolingFactor = 4.
-    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=1.)
-    betaest.beta0 = 1.
+    beta.coolingRate = 4
+    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)
     targmis = simpeg.Directives.TargetMisfit()
     targmis.target = survey.nD
+    saveModel = simpeg.Directives.SaveModelEveryIteration()
+    saveModel.fileName = 'Inversion_TargMisEqnD_smoothTrue'
     # Create an inversion object
     inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis])
 
@@ -126,9 +121,8 @@ def run(plotIt=True):
     mopt = inv.run(m_0)
 
     if plotIt:
-        fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
+        fig = MT.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
         fig.suptitle('Target - smooth true')
-        fig.axes[0].set_ylim([-10000,500])
         plt.show()
 
 if __name__ == '__main__':

SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py

@@ -1,428 +0,0 @@
-from scipy.constants import epsilon_0, mu_0
-import matplotlib.pyplot as plt
-import numpy as np
-from ipywidgets import *
-from SimPEG.EM.Utils import k, omega
-
-"""
-MT1D: n layered earth problem
-*****************************
-
-Author: Thibaut Astic
-Contact: thast@eos.ubc.ca
-Date: January 2016
-
-This code compute the analytic response of a n-layered Earth to a plane wave (Magneto-Tellurics).
-
-We start by looking at Maxwell's equations in the electric
-field \\\(\\\mathbf{E}\\) and the magnetic flux
-\\\(\\\mathbf{H}\\) to write the wave equations
-\\(\\ \nabla ^2  \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) &
-\\(\\ \nabla ^2  \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\)
-
-Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of
-\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\)
-\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\)
-
-With U and D the Up and Down components of the E-field.
-
-The iteration from one layer to another is ensure by:
-
-\\(\\  \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) =
-    P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\)
-
-And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0)
-and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability)
-
-The layer 0 is assumed to be the air layer.
-
-"""
-
-#Define a frquency range for a survey
-frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.)
-
-#Functions to create random physical Perties for a n-layered earth
-thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]),
-                                                     np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1)
-                                                            ,decimals =1))
-
-sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]),
-                                               np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1)
-                                                      ,decimals=3))
-
-mu  = lambda minmu, maxmu, nlayer: np.append(np.array([1.]),
-                                             np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1)
-                                                    ,decimals=1))
-
-eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]),
-                                               np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1)
-                                                                ,decimals=1))
-
-#Evaluate Impedance Z of a layer
-ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k
-
-#Complex Cole-Cole Conductivity - EM utils
-PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c)))
-
-#Converted thickness array into top of layer array
-top = lambda thick: np.cumsum(thick)
-
-#Propagation Matrix and theirs inverses
-
-#matrix T for transition of Up and Down components accross a layer
-T = lambda h,k: np.matrix([[np.exp(1j*k*h),0.],[0.,np.exp(-1j*k*h)]],dtype='complex_')
-
-Tinv = lambda h,k: np.matrix([[np.exp(-1j*k*h),0.],[0.,np.exp(1j*k*h)]],dtype='complex_')
-
-#transition of Up and Down components accross a layer
-UD_Z = lambda UD,z,zj,k : T((z-zj),k)*UD
-
-
-#matrix P relating Up and Down components with E and H fields
-P = lambda z: np.matrix([[1.,1,],[-1./z,1./z]],dtype='complex_')
-
-Pinv = lambda z: np.matrix([[1.,-z],[1.,z]],dtype='complex_')/2.
-
-
-#Time Variation of E and H
-E_ZT = lambda U,D,f,t : np.exp(1j*omega(f)*t)*(U+D)
-H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U)
-
-#Plot the configuration of the problem
-def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
-
-    topn = top(thick)
-    widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.)
-
-    ax.set_ylim([z.min(),z.max()])
-    ax.set_xlim([-widthg,widthg])
-
-    ax.set_ylabel("Depth (m)", fontsize=16.)
-    ax.yaxis.tick_right()
-    ax.yaxis.set_label_position("right")
-
-    #define filling for the different layers
-    hatches=['/'  , '+', 'x', '|' , '\\', '-'  , 'o' , 'O' , '.' , '*' ]
-
-    #Write the physical properties of air
-    ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)),
-            xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data',
-            xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data',
-            fontsize=14.)
-
-    ax.annotate(("$\epsilon_r$= %1i")%(eps[0]),
-            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
-            xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data',
-            fontsize=14.)
-
-    ax.annotate(("$\mu_r$= %1i")%(mu[0]),
-            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
-            xytext=(0, -np.abs(z.max())/3.), textcoords='data',
-            fontsize=14.)
-
-    #Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern
-    for i in range(1,len(topn)-1,1):
-        if topn[i] == topn[i+1]:
-            pass
-        else:
-            ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[i]*10**(3)),
-                xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data',
-                xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data',
-                fontsize=14.)
-
-            ax.annotate(("$\epsilon_r$= %1i")%(eps[i]),
-                xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data',
-                xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data',
-                fontsize=14.)
-
-            ax.annotate(("$\mu_r$= %1.2f")%(mu[i]),
-                xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data',
-                xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data',
-                fontsize=14.)
-
-            ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black')
-            ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10])
-
-    #Write the physical properties of the n-th layer and fill it with pattern
-    ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black')
-    ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10])
-
-    ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)),
-            xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data',
-            xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data',
-            fontsize=14.)
-
-    ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]),
-            xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data',
-            xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data',
-            fontsize=14.)
-
-    ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]),
-            xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data',
-            xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data',
-            fontsize=14.)
-
-    #plot Trees!
-    ax.annotate("",
-            xy=(widthg/2., -1.*z.max()/5.), xycoords='data',
-            xytext=(widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(widthg/2., -3./4.*z.max()/5.), xycoords='data',
-            xytext=(widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(widthg/2., -1./2.*z.max()/5.), xycoords='data',
-            xytext=(widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data',
-            xytext=(1.2*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.2*widthg/2., -3./4.*z.max()/5.), xycoords='data',
-            xytext=(1.2*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.2*widthg/2., -1./2.*z.max()/5.), xycoords='data',
-            xytext=(1.2*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.5*widthg/2., -1.*z.max()/5.), xycoords='data',
-            xytext=(1.5*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.5*widthg/2., -3./4.*z.max()/5.), xycoords='data',
-            xytext=(1.5*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.5*widthg/2., -1./2.*z.max()/5.), xycoords='data',
-            xytext=(1.5*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
-            )
-
-
-    ax.invert_yaxis()
-
-    return ax
-
-#Propagate Up and Down component for a certain frequency & evaluate E and H field
-
-def Propagate(f,H,sig,chg,taux,c,mu,eps,n):
-
-    sigcm = np.zeros_like(sig,dtype='complex_')
-
-    for j in range(1,len(sig)):
-        sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f)
-
-    K = k(f, sigcm, mu, eps)
-    Z = ImpZ(f,mu,K)
-
-    EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
-    UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
-
-    UD[1,-1] = 1.
-
-    for i in range(-2,-(n+2),-1):
-
-        UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1]
-        UD = UD/((np.abs(UD[0,:]+UD[1,:])).max())
-
-    for j in range(0,n+1):
-        EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j]
-
-    return UD, EH, Z ,K
-
-
-#Evaluate the apparent resistivity and phase for a frequency range
-def appres(F,H,sig,chg,taux,c,mu,eps,n):
-
-    Res = np.zeros_like(F)
-    Phase = np.zeros_like(F)
-    App_ImpZ= np.zeros_like(F,dtype='complex_')
-
-    for i in range(0,len(F)):
-
-        UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n)
-
-        App_ImpZ[i] = EH[0,1]/EH[1,1]
-
-        Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
-        Phase[i] = np.angle(App_ImpZ[i], deg = True)
-
-    return Res,Phase
-
-#Evaluate Up, Down components, E and H field, for a frequency range,
-#a discretized depth range and a time range (use to calculate envelope)
-def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample):
-
-    topc = top(H)
-
-    layer = np.zeros(len(zsample),dtype=np.int)-1
-
-    Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
-    Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
-    Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
-    Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
-    UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_')
-
-    for i in range(0,n+1,1):
-        layer = layer+(zsample>=topc[i])*1
-
-    for j in range(0,len(F)):
-
-        UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n)
-
-        for p in range(0,len(zsample)):
-
-            UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]])
-
-            for q in range(0,len(tsample)):
-
-                Exzt[p,q]  = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F)
-                Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F)
-                Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
-                Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
-
-    return  Exzt,Hyzt,Uz,Dz,UDaux,layer
-
-
-#Function to Plot Apparent Resistivity and Phase
-def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
-
-    Res, Phase = appres(F,H,sig,chg,taux,c,mu,eps,n)
-
-    fig,ax = plt.subplots(1,2,figsize=(16,10))
-
-    ax[0].scatter(Res,F,color='black')
-    ax[0].set_xscale('Log')
-    ax[0].set_yscale('Log')
-    ax[0].set_xlim([10.**(np.log10(Res.min())-1.),10.**(np.log10(Res.max())+1.)])
-    ax[0].set_ylim([F.min(),F.max()])
-    ax[0].set_xlabel('Apparent Resistivity (Ohm*m)',fontsize=16.,color="black")
-    ax[0].set_ylabel('Frequency (Hz)',fontsize=16.)
-    ax[0].grid(which='major')
-
-    ax0 = ax[0].twiny()
-
-    ax0.set_xlim([0.,90.])
-    ax0.set_ylim([F.min(),F.max()])
-    ax0.scatter(Phase,F,color='purple')
-    ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple")
-
-    zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.)
-
-    ax[0].tick_params(labelsize=16)
-    ax[1].tick_params(labelsize=16)
-    ax0.tick_params(labelsize=16)
-
-    if PlotEnvelope:
-
-        widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0)
-        fenvelope1n=np.ones(100)*fenvelope
-        ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black')
-
-        tc=np.arange(0.,1./fenvelope,0.01/(fenvelope))
-        Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc)
-
-        ax1=ax[1].twiny()
-
-        ax[1].tick_params(labelsize=16)
-        ax1.tick_params(labelsize=16)
-
-        ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16)
-
-        ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16)
-
-        ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))),
-                      np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))),
-                      color='blue', alpha=0.1)
-
-        ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))),
-                      np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))),
-                      color='red', alpha=0.1)
-
-        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc)
-        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
-        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
-    else:
-        print 'No envelop (if True, might be slow)'
-        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc)
-        ax[1].get_xaxis().set_ticks([])
-
-    plt.show()
-
-#Interactive MT for Notebook
-def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope):
-
-    frangn=frange(-5,5,100.)
-    sig3= np.array([0.,0.001,0.1, 0.001])
-    thick3 = np.array([120000.,50.,50.])
-    eps3=np.array([1.,1.,1.,1])
-    mu3=np.array([1.,1.,1.,1])
-    chg3=np.array([0.,0.1,0.,0.2])
-    chg3_0=np.array([0.,0.1,0.,0.])
-    taux3=np.array([0.,0.1,0.,0.1])
-    c3=np.array([1.,1.,1.,1.])
-
-    sig3[1]=sigl1
-    sig3[1]=10.**sig3[1]
-    sig3[2]=sigl2
-    sig3[2]=10.**sig3[2]
-    sig3[3]=sigl3
-    sig3[3]=10.**sig3[3]
-    mu3[1]=mul1
-    mu3[2]=mul2
-    mu3[3]=mul3
-    eps3[1]=epsl1
-    eps3[2]=epsl2
-    eps3[3]=epsl3
-    thick3[1]=h1
-    thick3[2]=h2
-
-    PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
-
-
-def run(n=3,plotIt=True):
-    # something to make a plot
-
-    F = frange(-5.,5.,20)
-    H = thick(50.,100.,n)
-    sign = sig(-5.,0.,n)
-    mun = mu(1.,2.,n)
-    epsn = eps(1.,9.,n)
-    chg = np.zeros_like(sign)
-    taux = np.zeros_like(sign)
-    c = np.zeros_like(sign)
-
-    Res, Phase = appres(F,H,sign,chg,taux,c,mun,epsn,n)
-
-    if plotIt:
-
-        PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True)
-
-    return Res, Phase
-
-if __name__ == '__main__':
-    run()
-
-
-
-
-

SimPEG/Examples/MT_3D_Foward.py

@@ -2,7 +2,7 @@
 
 # Import
 import SimPEG as simpeg
-from SimPEG import NSEM
+from SimPEG import MT
 import numpy as np
 try:
     from pymatsolver import MumpsSolver as Solver
@@ -12,7 +12,7 @@ except:
 def run(plotIt=True, nFreq=1):
     """
         MT: 3D: Forward
-        =======================
+        ===============
 
         Forward model 3D MT data.
 
@@ -37,25 +37,24 @@ def run(plotIt=True, nFreq=1):
     for loc in rx_loc:
         # NOTE: loc has to be a (1,3) np.ndarray otherwise errors accure
         for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(NSEM.Rx(simpeg.mkvc(loc,2).T,rxType))
+            rxList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
     # Source list
     srcList =[]
     for freq in np.logspace(3,-3,nFreq):
-        srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
+        srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
     # Survey MT
-    survey = NSEM.Survey(srcList)
+    survey = MT.Survey(srcList)
 
     ## Setup the problem object
-    problem = NSEM.Problem3D_ePrimSec(M, sigmaPrimary=sigBG)
+    problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG, Solver=Solver)
     problem.pair(survey)
-    problem.Solver = Solver
 
     # Calculate the data
     fields = problem.fields(sig)
     dataVec = survey.eval(fields)
 
     # Make the data
-    mtData = NSEM.Data(survey,dataVec)
+    mtData = MT.Data(survey, dataVec)
     # Add plots
     if plotIt:
         pass

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,8 +2,12 @@ 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
+        =======================
+
         Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below
         a topographic surface.
 
@@ -13,27 +17,25 @@ def run(plotIt=False, nx = 5, ny = 5):
     xtopo = np.linspace(mesh.gridN[:,0].min(), mesh.gridN[:,0].max())
     topo = 0.4*np.sin(xtopo*5) # define a topographic surface
 
-    Topo = np.hstack([Utils.mkvc(xtopo,2),Utils.mkvc(topo,2)]) #make it an array
+    Topo = np.hstack([Utils.mkvc(xtopo,2), Utils.mkvc(topo,2)]) #make it an array
 
-    indcc = surface2ind_topo(mesh, Topo,'CC')
+    indcc = surface2ind_topo(mesh, Topo, 'CC')
 
     if plotIt:
         from matplotlib.pylab import plt
         from scipy.interpolate import interp1d
-        fig, ax = plt.subplots(1,1,figsize=(6,6))
+        fig, ax = plt.subplots(1,1, figsize=(6,6))
         mesh.plotGrid(ax=ax, nodes=True, centers=True)
         ax.plot(xtopo,topo,'k',linewidth=1)
-        # ax.plot(mesh.vectorNx, interp1d(xtopo,topo)(mesh.vectorNx),'--k',linewidth=3)
         ax.plot(mesh.vectorCCx, interp1d(xtopo,topo)(mesh.vectorCCx),'--k',linewidth=3)
 
-
         aveN2CC = Utils.sdiag(mesh.aveN2CC.T.sum(1))*mesh.aveN2CC.T
         a = aveN2CC * indcc
         a[a > 0] = 1.
         a[a < 0.25] = np.nan
         a = a.reshape(mesh.vnN, order='F')
         masked_array = np.ma.array(a, mask=np.isnan(a))
-        ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap = plt.cm.gray,alpha=0.2)
+        ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap=plt.cm.gray, alpha=0.2)
         plt.show()
 
 

SimPEG/Examples/__init__.py

@@ -1,29 +1,30 @@
 # Run this file to add imports.
 
 ##### AUTOIMPORTS #####
-import EM_FDEM_1D_Inversion
-import Mesh_QuadTree_Creation
-import EM_TDEM_1D_Inversion
-import Mesh_QuadTree_FaceDiv
-import Mesh_Tensor_Creation
-import FLOW_Richards_1D_Celia1990
+import DC_Analytic_Dipole
 import DC_Forward_PseudoSection
-import Mesh_Operators_CahnHilliard
-import Mesh_Basic_Types
+import EM_FDEM_1D_Inversion
+import EM_FDEM_Analytic_MagDipoleWholespace
+import EM_Schenkel_Morrison_Casing
+import EM_TDEM_1D_Inversion
+import FLOW_Richards_1D_Celia1990
 import Inversion_IRLS
 import Inversion_Linear
-import EM_Schenkel_Morrison_Casing
-import MT_3D_Foward
+import Maps_ComboMaps
+import Maps_Mesh2Mesh
 import Mesh_Basic_ForwardDC
-import MT_1D_ForwardAndInversion
-import Utils_surface2ind_topo
-import MT_1D_analytic_nlayer_Earth
-import EM_FDEM_Analytic_MagDipoleWholespace
 import Mesh_Basic_PlotImage
-import DC_Analytic_Dipole
+import Mesh_Basic_Types
+import Mesh_Operators_CahnHilliard
+import Mesh_QuadTree_Creation
+import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
+import Mesh_Tensor_Creation
+import MT_1D_ForwardAndInversion
+import MT_3D_Foward
+import Utils_surface2ind_topo
 
-__examples__ = ["EM_FDEM_1D_Inversion", "Mesh_QuadTree_Creation", "EM_TDEM_1D_Inversion", "Mesh_QuadTree_FaceDiv", "Mesh_Tensor_Creation", "FLOW_Richards_1D_Celia1990", "DC_Forward_PseudoSection", "Mesh_Operators_CahnHilliard", "Mesh_Basic_Types", "Inversion_IRLS", "Inversion_Linear", "EM_Schenkel_Morrison_Casing", "MT_3D_Foward", "Mesh_Basic_ForwardDC", "MT_1D_ForwardAndInversion", "Utils_surface2ind_topo", "MT_1D_analytic_nlayer_Earth", "EM_FDEM_Analytic_MagDipoleWholespace", "Mesh_Basic_PlotImage", "DC_Analytic_Dipole", "Mesh_QuadTree_HangingNodes"]
+__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 #####
 
@@ -39,7 +40,7 @@ if __name__ == '__main__':
 
     # Create the examples dir in the docs folder.
     fName = os.path.realpath(__file__)
-    docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'examples'])
+    docExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['docs', 'content', 'examples'])
     shutil.rmtree(docExamplesDir)
     os.makedirs(docExamplesDir)
 
@@ -96,12 +97,12 @@ if __name__ == '__main__':
     from SimPEG import Examples
     Examples.%s.run()
 
-.. literalinclude:: ../../SimPEG/Examples/%s.py
+.. literalinclude:: ../../../SimPEG/Examples/%s.py
     :language: python
     :linenos:
 """%(name,doc,name,name)
 
-        rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'examples', name + '.rst']))
+        rst = os.path.sep.join((filePath.split(os.path.sep)[:-3] + ['docs', 'content', 'examples', name + '.rst']))
 
         print 'Creating: %s.rst'%name
         f = open(rst, 'w')

SimPEG/FLOW/Richards/Empirical.py

@@ -31,7 +31,7 @@ class NonLinearMap(object):
         """
             :param numpy.array u: fields
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
 
             The *transform* changes the model into the physical property.
@@ -44,7 +44,7 @@ class NonLinearMap(object):
         """
             :param numpy.array u: fields
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
 
             The *transform* changes the model into the physical property.

SimPEG/MT/BaseMT.py

@@ -0,0 +1,132 @@
+from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
+from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
+from SurveyMT import Survey, Data
+from FieldsMT import BaseMTFields
+
+
+class BaseMTProblem(BaseFDEMProblem):
+    """
+        Base class for all Natural source problems.
+    """
+
+    def __init__(self, mesh, **kwargs):
+        BaseFDEMProblem.__init__(self, mesh, **kwargs)
+        Utils.setKwargs(self, **kwargs)
+    # Set the default pairs of the problem
+    surveyPair = Survey
+    dataPair = Data
+    fieldsPair = BaseMTFields
+
+    # Set the solver
+    Solver = SimpegSolver
+    solverOpts = {}
+
+    verbose = False
+    # Notes:
+    # Use the forward and devs from BaseFDEMProblem
+    # Might need to add more stuff here.
+
+    ## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
+    def Jvec(self, m, v, f=None):
+        """
+        Function to calculate the data sensitivities dD/dm times a vector.
+
+            :param numpy.ndarray m (nC, 1) - conductive model
+            :param numpy.ndarray v (nC, 1) - random vector
+            :param MTfields object (optional) - MT fields object, if not given it is calculated
+            :rtype: MTdata object
+            :return: Data sensitivities wrt m
+        """
+
+        # Calculate the fields
+        if f is None:
+           f= self.fields(m)
+        # Set current model
+        self.curModel = m
+        # Initiate the Jv object
+        Jv = self.dataPair(self.survey)
+
+        # Loop all the frequenies
+        for freq in self.survey.freqs:
+            dA_du = self.getA(freq) #
+
+            dA_duI = self.Solver(dA_du, **self.solverOpts)
+
+            for src in self.survey.getSrcByFreq(freq):
+                # We need fDeriv_m = df/du*du/dm + df/dm
+                # Construct du/dm, it requires a solve
+                # NOTE: need to account for the 2 polarizations in the derivatives.
+                f_src = f[src,:]
+                # dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
+                dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns.
+                dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
+                if dRHS_dm is None:
+                    du_dm = dA_duI * ( -dA_dm )
+                else:
+                    du_dm = dA_duI * ( -dA_dm + dRHS_dm )
+                # Calculate the projection derivatives
+                for rx in src.rxList:
+                    # Get the projection derivative
+                    # v should be of size 2*nE (for 2 polarizations)
+                    PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
+                    Jv[src, rx] = PDeriv_u(mkvc(du_dm))
+            dA_duI.clean()
+        # Return the vectorized sensitivities
+        return mkvc(Jv)
+
+    def Jtvec(self, m, v, f=None):
+        """
+        Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
+
+            :param numpy.ndarray m (nC, 1) - conductive model
+            :param numpy.ndarray v (nD, 1) - vector
+            :param MTfields object u (optional) - MT fields object, if not given it is calculated
+            :rtype: MTdata object
+            :return: Data sensitivities wrt m
+        """
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        # Ensure v is a data object.
+        if not isinstance(v, self.dataPair):
+            v = self.dataPair(self.survey, v)
+
+        Jtv = np.zeros(m.size)
+
+        for freq in self.survey.freqs:
+            AT = self.getA(freq).T
+
+            ATinv = self.Solver(AT, **self.solverOpts)
+
+            for src in self.survey.getSrcByFreq(freq):
+                ftype = self._fieldType + 'Solution'
+                f_src = f[src, :]
+
+                for rx in src.rxList:
+                    # Get the adjoint evalDeriv
+                    # PTv needs to be nE,
+                    PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
+                    # Get the
+                    dA_duIT = ATinv * PTv
+                    dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
+                    dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
+                    # Make du_dmT
+                    if dRHS_dmT is None:
+                        du_dmT = -dA_dmT
+                    else:
+                        du_dmT = -dA_dmT + dRHS_dmT
+                    # Select the correct component
+                    # du_dmT needs to be of size nC,
+                    real_or_imag = rx.projComp
+                    if real_or_imag == 'real':
+                        Jtv +=  du_dmT.real
+                    elif real_or_imag == 'imag':
+                        Jtv +=  -du_dmT.real
+                    else:
+                        raise Exception('Must be real or imag')
+            # Clean the factorization, clear memory.
+            ATinv.clean()
+        return Jtv

SimPEG/MT/FieldsMT.py

@@ -4,21 +4,18 @@ import sys
 from numpy.lib import recfunctions as recFunc
 from SimPEG.EM.Utils import omega
 
-
 ##############
 ### Fields ###
 ##############
-class BaseNSEMFields(Problem.Fields):
-    """Field Storage for a NSEM survey."""
+class BaseMTFields(Problem.Fields):
+    """Field Storage for a MT survey."""
     knownFields = {}
     dtype = complex
 
-###########
-# 1D Fields
-###########
-class Fields1D_ePrimSec(BaseNSEMFields):
+
+class Fields1D_e(BaseMTFields):
     """
-    Fields storage for the 1D NSEM solution.
+    Fields storage for the 1D MT solution.
     """
     knownFields = {'e_1dSolution':'F'}
     aliasFields = {
@@ -31,119 +28,7 @@ class Fields1D_ePrimSec(BaseNSEMFields):
                   }
 
     def __init__(self,mesh,survey,**kwargs):
-        BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
-
-    def _ePrimary(self, eSolution, srcList):
-        ePrimary = np.zeros_like(eSolution)
-        for i, src in enumerate(srcList):
-            ep = src.ePrimary(self.survey.prob)
-            if ep is not None:
-                ePrimary[:,i] = ep[:,-1]
-        return ePrimary
-
-    def _eSecondary(self, eSolution, srcList):
-        return eSolution
-
-    def _e(self, eSolution, srcList):
-        return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
-
-    def _eDeriv_u(self, src, du_dm_v, adjoint = False):
-
-
-        return Utils.Identity()*du_dm_v
-
-    def _eDeriv_m(self, src, v, adjoint = False):
-        # assuming primary does not depend on the model
-        return Utils.Zero()
-
-    def _bPrimary(self, eSolution, srcList):
-        bPrimary = np.zeros([self.survey.mesh.nE,eSolution.shape[1]], dtype = complex)
-        for i, src in enumerate(srcList):
-            bp = src.bPrimary(self.survey.prob)
-            if bp is not None:
-                bPrimary[:,i] += bp[:,-1]
-        return bPrimary
-
-    def _bSecondary(self, eSolution, srcList):
-        C = self.mesh.nodalGrad
-        b = (C * eSolution)
-        for i, src in enumerate(srcList):
-            b[:,i] *= - 1./(1j*omega(src.freq))
-            # There is no magnetic source in the MT problem
-            # S_m, _ = src.eval(self.survey.prob)
-            # if S_m is not None:
-            #     b[:,i] += 1./(1j*omega(src.freq)) * S_m
-        return b
-
-    def _b(self, eSolution, srcList):
-        return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
-
-    def _bSecondaryDeriv_u(self, src, v, adjoint = False):
-        C = self.mesh.nodalGrad
-        if adjoint:
-            return - 1./(1j*omega(src.freq)) * (C.T * v)
-        return - 1./(1j*omega(src.freq)) * (C * v)
-
-    def _bSecondaryDeriv_m(self, src, v, adjoint = False):
-        # Doesn't depend on m
-        # _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
-        # S_eDeriv = S_eDeriv(v)
-        # if S_eDeriv is not None:
-        #     return 1./(1j * omega(src.freq)) * S_eDeriv
-        return None
-
-    def _bDeriv_u(self, src, v, adjoint=False):
-        # Primary does not depend on u
-        return self._bSecondaryDeriv_u(src, v, adjoint)
-
-    def _bDeriv_m(self, src, v, adjoint=False):
-        # Assuming the primary does not depend on the model
-        return self._bSecondaryDeriv_m(src, v, adjoint)
-
-    def _fDeriv_u(self, src, v, adjoint=False):
-        """
-        Derivative of the fields object wrt u.
-
-        :param NSEMsrc src: NSEM source
-        :param numpy.ndarray v: random vector of f_sol.size
-        This function stacks the fields derivatives appropriately
-
-        return a vector of size (nreEle+nrbEle)
-        """
-
-        de_du = v #Utils.spdiag(np.ones((self.nF,)))
-        db_du = self._bDeriv_u(src, v, adjoint)
-        # Return the stack
-        # This doesn't work...
-        return np.vstack((de_du,db_du))
-
-    def _fDeriv_m(self, src, v, adjoint=False):
-        """
-        Derivative of the fields object wrt m.
-
-        This function stacks the fields derivatives appropriately
-        """
-        return None
-
-
-class Fields1D_eTotal(BaseNSEMFields):
-    """
-    Fields storage for the 1D NSEM solution solved with for a total domain formulation.
-
-    Used in conjuction with Problem1D_eTotal.
-    """
-    knownFields = {'e_1dSolution':'F'}
-    aliasFields = {
-                    'e_1d' : ['e_1dSolution','F','_e'],
-                    'e_1dPrimary' : ['e_1dSolution','F','_ePrimary'],
-                    'e_1dSecondary' : ['e_1dSolution','F','_eSecondary'],
-                    'b_1d' : ['e_1dSolution','E','_b'],
-                    'b_1dPrimary' : ['e_1dSolution','E','_bPrimary'],
-                    'b_1dSecondary' : ['e_1dSolution','E','_bSecondary']
-                  }
-
-    def __init__(self,mesh,survey,**kwargs):
-        BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
+        BaseMTFields.__init__(self,mesh,survey,**kwargs)
 
     def _ePrimary(self, eSolution, srcList):
         ePrimary = np.zeros_like(eSolution)
@@ -214,7 +99,7 @@ class Fields1D_eTotal(BaseNSEMFields):
         """
         Derivative of the fields object wrt u.
 
-        :param NSEMsrc src: NSEM source
+        :param MTsrc src: MT source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 
@@ -235,18 +120,9 @@ class Fields1D_eTotal(BaseNSEMFields):
         """
         return None
 
-
-###########
-# 2D Fields
-###########
-
-
-###########
-# 3D Fields
-###########
-class Fields3D_ePrimSec(BaseNSEMFields):
+class Fields3D_e(BaseMTFields):
     """
-    Fields storage for the 3D NSEM solution. Labels polarizations by px and py.
+    Fields storage for the 3D MT solution. Labels polarizations by px and py.
 
         :param SimPEG object mesh: The solution mesh
         :param SimPEG object survey: A survey object
@@ -271,7 +147,7 @@ class Fields3D_ePrimSec(BaseNSEMFields):
                   }
 
     def __init__(self,mesh,survey,**kwargs):
-        BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
+        BaseMTFields.__init__(self,mesh,survey,**kwargs)
 
     def _e_pxPrimary(self, e_pxSolution, srcList):
         e_pxPrimary = np.zeros_like(e_pxSolution)
@@ -352,7 +228,7 @@ class Fields3D_ePrimSec(BaseNSEMFields):
         b = (C * e_pxSolution)
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
-            # There is no magnetic source in the NSEM problem
+            # There is no magnetic source in the MT problem
             # S_m, _ = src.eval(self.survey.prob)
             # if S_m is not None:
             #     b[:,i] += 1./(1j*omega(src.freq)) * S_m
@@ -363,7 +239,7 @@ class Fields3D_ePrimSec(BaseNSEMFields):
         b = (C * e_pySolution)
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
-            # There is no magnetic source in the NSEM problem
+            # There is no magnetic source in the MT problem
             # S_m, _ = src.eval(self.survey.prob)
             # if S_m is not None:
             #     b[:,i] += 1./(1j*omega(src.freq)) * S_m
@@ -426,7 +302,7 @@ class Fields3D_ePrimSec(BaseNSEMFields):
         """
         Derivative of the fields object wrt u.
 
-        :param NSEMsrc src: NSEM source
+        :param MTsrc src: MT source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 
@@ -443,7 +319,7 @@ class Fields3D_ePrimSec(BaseNSEMFields):
         """
         Derivative of the fields object wrt u.
 
-        :param NSEMsrc src: NSEM source
+        :param MTsrc src: MT source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 

SimPEG/MT/Problem1D/Probs.py

@@ -0,0 +1,291 @@
+from SimPEG.EM.Utils import omega
+from SimPEG import mkvc
+from scipy.constants import mu_0
+from SimPEG.MT.BaseMT import BaseMTProblem
+from SimPEG.MT.SurveyMT import Survey, Data
+from SimPEG.MT.FieldsMT import Fields1D_e
+from SimPEG.MT.Utils.MT1Danalytic import getEHfields
+import numpy as np
+import multiprocessing, sys, time
+
+
+class eForm_psField(BaseMTProblem):
+    """
+    A MT problem soving a e formulation and primary/secondary fields decomposion.
+
+    By eliminating the magnetic flux density using
+
+        .. math ::
+
+            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
+
+
+    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
+
+    .. math ::
+        \\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
+    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
+
+    The primary field is estimated from a background model (commonly half space ).
+
+
+    """
+    # From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
+    _fieldType = 'e_1d'
+    _eqLocs    = 'EF'
+    _sigmaPrimary = None
+
+
+    def __init__(self, mesh, **kwargs):
+        BaseMTProblem.__init__(self, mesh, **kwargs)
+        self.fieldsPair = Fields1D_e
+        # self._sigmaPrimary = sigmaPrimary
+    @property
+    def MeMui(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MeMui', None) is None:
+            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
+        return self._MeMui
+
+    @property
+    def MfSigma(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MfSigma', None) is None:
+            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
+        return self._MfSigma
+
+    @property
+    def sigmaPrimary(self):
+        """
+        A background model, use for the calculation of the primary fields.
+
+        """
+        return self._sigmaPrimary
+
+    @sigmaPrimary.setter
+    def sigmaPrimary(self, val):
+        # Note: TODO add logic for val, make sure it is the correct size.
+        self._sigmaPrimary = val
+
+    def getA(self, freq):
+        """
+            Function to get the A matrix.
+
+            :param float freq: Frequency
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+
+        # Note: need to use the code above since in the 1D problem I want
+        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
+        # Possible that _fieldType and _eqLocs can fix this
+        MeMui = self.MeMui
+        MfSigma = self.MfSigma
+        C = self.mesh.nodalGrad
+        # Make A
+        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
+        # Either return full or only the inner part of A
+        return A
+
+    def getADeriv_m(self, freq, u, v, adjoint=False):
+        """
+        The derivative of A wrt sigma
+        """
+
+        dsig_dm = self.curModel.sigmaDeriv
+        MeMui = self.MeMui
+        #
+        u_src = u['e_1dSolution']
+        dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
+        if adjoint:
+            return 1j * omega(freq) * (  dMfSigma_dm.T * v )
+        # Note: output has to be nN/nF, not nC/nE.
+        # v should be nC
+        return 1j * omega(freq) * ( dMfSigma_dm * v )
+
+    def getRHS(self, freq):
+        """
+            Function to return the right hand side for the system.
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
+            :return: RHS for 1 polarizations, primary fields
+        """
+
+        # Get sources for the frequncy(polarizations)
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_e = Src.S_e(self)
+        return -1j * omega(freq) * S_e
+
+    def getRHSDeriv_m(self, freq, v, adjoint=False):
+        """
+        The derivative of the RHS wrt sigma
+        """
+
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
+        return -1j * omega(freq) * S_eDeriv
+
+    def fields(self, m):
+        '''
+        Function to calculate all the fields for the model m.
+
+        :param np.ndarray (nC,) m: Conductivity model
+        '''
+        # Set the current model
+        self.curModel = m
+
+        F = Fields1D_e(self.mesh, self.survey)
+        for freq in self.survey.freqs:
+            if self.verbose:
+                startTime = time.time()
+                print 'Starting work for {:.3e}'.format(freq)
+                sys.stdout.flush()
+            A = self.getA(freq)
+            rhs  = self.getRHS(freq)
+            Ainv = self.Solver(A, **self.solverOpts)
+            e_s = Ainv * rhs
+
+            # Store the fields
+            Src = self.survey.getSrcByFreq(freq)[0]
+            # NOTE: only store the e_solution(secondary), all other components calculated in the fields object
+            F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
+
+            # Note curl e = -iwb so b = -curl e /iw
+            # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
+            # F[Src, 'b_1d'] = b[:,1]
+            if self.verbose:
+                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                sys.stdout.flush()
+        return F
+
+# Note this is not fully functional.
+# Missing:
+# Fields class corresponding to the fields
+# Update Jvec and Jtvec to include all the derivatives components
+# Other things ...
+class eForm_TotalField(BaseMTProblem):
+    """
+    A MT problem solving a e formulation and a Total bondary domain decompostion.
+
+    Solves the equation:
+
+    Math:
+
+
+    """
+
+    # From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
+    _fieldType = 'e'
+    _eqLocs    = 'EF'
+
+
+    def __init__(self, mesh, **kwargs):
+        BaseMTProblem.__init__(self, mesh, **kwargs)
+    @property
+    def MeMui(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MeMui', None) is None:
+            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
+        return self._MeMui
+
+    @property
+    def MfSigma(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MfSigma', None) is None:
+            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
+        return self._MfSigma
+
+    def getA(self, freq, full=False):
+        """
+            Function to get the A matrix.
+
+            :param float freq: Frequency
+            :param logic full: Return full A or the inner part
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+
+        MeMui = self.MeMui
+        MfSigma = self.MfSigma
+        # Note: need to use the code above since in the 1D problem I want
+        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
+        # Possible that _fieldType and _eqLocs can fix this
+        # MeMui = self.MfMui
+        # MfSigma = self.MfSigma
+        C = self.mesh.nodalGrad
+        # Make A
+        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
+        # Either return full or only the inner part of A
+        if full:
+            return A
+        else:
+            return A[1:-1,1:-1]
+
+    def getADeriv_m(self, freq, u, v, adjoint=False):
+        raise NotImplementedError('getADeriv is not implemented')
+
+    def getRHS(self, freq):
+        """
+            Function to return the right hand side for the system.
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
+            :return: RHS for both polarizations, primary fields
+        """
+        # Get sources for the frequency
+        # NOTE: Need to use the source information, doesn't really apply in 1D
+        src = self.survey.getSrcByFreq(freq)
+        # Get the full A
+        A = self.getA(freq,full=True)
+        # Define the outer part of the solution matrix
+        Aio = A[1:-1,[0,-1]]
+        Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
+        Etot = (Ed + Eu)
+        sourceAmp = 1.0
+        Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
+        ## Note: The analytic solution is derived with e^iwt
+        eBC = np.r_[Etot[0],Etot[-1]]
+        # The right hand side
+
+        return -Aio*eBC, eBC
+
+    def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
+        raise NotImplementedError('getRHSDeriv not implemented yet')
+        return None
+
+    def fields(self, m):
+        '''
+        Function to calculate all the fields for the model m.
+
+        :param np.ndarray (nC,) m: Conductivity model
+        :param np.ndarray (nC,) m_back: Background conductivity model
+        '''
+        self.curModel = m
+        # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
+
+        F = Fields1D_e(self.mesh, self.survey)
+        for freq in self.survey.freqs:
+            if self.verbose:
+                startTime = time.time()
+                print 'Starting work for {:.3e}'.format(freq)
+                sys.stdout.flush()
+            A = self.getA(freq)
+            rhs, e_o = self.getRHS(freq)
+            Ainv = self.Solver(A, **self.solverOpts)
+            e_i = Ainv * rhs
+            e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
+            # Store the fields
+            Src = self.survey.getSrcByFreq(freq)
+            # NOTE: only store e fields
+            F[Src, 'e_1dSolution'] = e[:,0]
+            if self.verbose:
+                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                sys.stdout.flush()
+        return F

SimPEG/MT/Problem1D/__init__.py

@@ -0,0 +1 @@
+from Probs import eForm_TotalField, eForm_psField

SimPEG/MT/Problem2D/__init__.py

@@ -0,0 +1 @@
+pass

SimPEG/MT/Problem3D/Probs.py

@@ -0,0 +1,138 @@
+from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver
+from SimPEG.EM.Utils import omega
+from scipy.constants import mu_0
+from SimPEG.MT.BaseMT import BaseMTProblem
+from SimPEG.MT.SurveyMT import Survey, Data
+from SimPEG.MT.FieldsMT import Fields3D_e
+import multiprocessing, sys, time
+
+
+
+class eForm_ps(BaseMTProblem):
+    """
+    A MT problem solving a e formulation and a primary/secondary fields decompostion.
+
+    By eliminating the magnetic flux density using
+
+        .. math ::
+
+            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
+
+
+    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
+
+    .. math ::
+        \\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
+    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
+
+    The primary field is estimated from a background model (commonly as a 1D model).
+
+    """
+
+    # From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
+    _fieldType = 'e'
+    _eqLocs    = 'FE'
+    fieldsPair = Fields3D_e
+    _sigmaPrimary = None
+
+    def __init__(self, mesh, **kwargs):
+        BaseMTProblem.__init__(self, mesh, **kwargs)
+
+    @property
+    def sigmaPrimary(self):
+        """
+        A background model, use for the calculation of the primary fields.
+
+        """
+        return self._sigmaPrimary
+    @sigmaPrimary.setter
+    def sigmaPrimary(self, val):
+        # Note: TODO add logic for val, make sure it is the correct size.
+        self._sigmaPrimary = val
+
+    def getA(self, freq):
+        """
+            Function to get the A system.
+
+            :param float freq: Frequency
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+        Mmui = self.MfMui
+        Msig = self.MeSigma
+        C = self.mesh.edgeCurl
+
+        return C.T*Mmui*C + 1j*omega(freq)*Msig
+
+    def getADeriv_m(self, freq, u, v, adjoint=False):
+        """
+        Calculate the derivative of A wrt m.
+
+        """
+
+        # This considers both polarizations and returns a nE,2 matrix for each polarization
+        if adjoint:
+            dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v
+        else:
+            # Need a nE,2 matrix to be returned
+            dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
+        return 1j * omega(freq) * dMe_dsigV
+
+
+    def getRHS(self, freq):
+        """
+        Function to return the right hand side for the system.
+
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
+            :return: RHS for both polarizations, primary fields
+        """
+
+        # Get sources for the frequncy(polarizations)
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_e = Src.S_e(self)
+        return -1j * omega(freq) * S_e
+
+    def getRHSDeriv_m(self, freq, v, adjoint=False):
+        """
+        The derivative of the RHS with respect to sigma
+        """
+
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
+        return -1j * omega(freq) * S_eDeriv
+
+    def fields(self, m):
+        '''
+        Function to calculate all the fields for the model m.
+
+        :param np.ndarray (nC,) m: Conductivity model
+        '''
+        # Set the current model
+        self.curModel = m
+
+        F = Fields3D_e(self.mesh, self.survey)
+        for freq in self.survey.freqs:
+            if self.verbose:
+                startTime = time.time()
+                print 'Starting work for {:.3e}'.format(freq)
+                sys.stdout.flush()
+            A = self.getA(freq)
+            rhs  = self.getRHS(freq)
+            # Solve the system
+            Ainv = self.Solver(A, **self.solverOpts)
+            e_s = Ainv * rhs
+
+            # Store the fields
+            Src = self.survey.getSrcByFreq(freq)[0]
+            # Store the fieldss
+            F[Src, 'e_pxSolution'] = e_s[:,0]
+            F[Src, 'e_pySolution'] = e_s[:,1]
+            # Note curl e = -iwb so b = -curl/iw
+
+            if self.verbose:
+                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                sys.stdout.flush()
+            Ainv.clean()
+        return F
+

SimPEG/MT/Problem3D/__init__.py

@@ -0,0 +1 @@
+from Probs import eForm_ps

SimPEG/MT/SrcMT.py

@@ -11,9 +11,9 @@ import sys
 ###   Sources ###
 #################
 
-class BaseNSEMSrc(FDEMBaseSrc):
+class BaseMTSrc(FDEMBaseSrc):
     '''
-    Sources for the NSEM problem.
+    Sources for the MT problem.
     Use the SimPEG BaseSrc, since the source fields share properties with the transmitters.
 
     :param float freq: The frequency of the source
@@ -29,28 +29,28 @@ class BaseNSEMSrc(FDEMBaseSrc):
         FDEMBaseSrc.__init__(self, rxList)
 
 # 1D sources
-class polxy_1DhomotD(BaseNSEMSrc):
+class polxy_1DhomotD(BaseMTSrc):
     """
-    NSEM source for both polarizations (x and y) for the total Domain.
+    MT source for both polarizations (x and y) for the total Domain.
 
     It calculates fields calculated based on conditions on the boundary of the domain.
     """
     def __init__(self, rxList, freq):
-        BaseNSEMSrc.__init__(self, rxList, freq)
+        BaseMTSrc.__init__(self, rxList, freq)
 
 
     # TODO: need to add the  primary fields calc and source terms into the problem.
 
 # Need to implement such that it works for all dims.
-class polxy_1Dprimary(BaseNSEMSrc):
+class polxy_1Dprimary(BaseMTSrc):
     """
-    NSEM source for both polarizations (x and y) given a 1D primary models.
+    MT source for both polarizations (x and y) given a 1D primary models.
     It assigns fields calculated from the 1D model as fields in the full space of the problem.
     """
     def __init__(self, rxList, freq):
         # assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
         self.sigma1d = None
-        BaseNSEMSrc.__init__(self, rxList, freq)
+        BaseMTSrc.__init__(self, rxList, freq)
         # Hidden property of the ePrimary
         self._ePrimary = None
 
@@ -128,15 +128,15 @@ class polxy_1Dprimary(BaseNSEMSrc):
             # v should be nC size
             return MsigmaDeriv * v
 
-class polxy_3Dprimary(BaseNSEMSrc):
+class polxy_3Dprimary(BaseMTSrc):
     """
-    NSEM source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
+    MT source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
     as fields in the full space of the problem.
     """
     def __init__(self, rxList, freq):
         # assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
         self.sigmaPrimary = None
-        BaseNSEMSrc.__init__(self, rxList, freq)
+        BaseMTSrc.__init__(self, rxList, freq)
         # Hidden property of the ePrimary
         self._ePrimary = None
 

SimPEG/MT/SurveyMT.py

@@ -4,7 +4,7 @@ from SimPEG.EM.Utils import omega
 from scipy.constants import mu_0
 from numpy.lib import recfunctions as recFunc
 from Utils import rec2ndarr
-import SrcNSEM
+import SrcMT
 import sys
 
 #################
@@ -63,9 +63,9 @@ class Rx(SimPEGsurvey.BaseRx):
         '''
         Project the fields to natural source data.
 
-            :param SrcNSEM src: The source of the fields to project
+            :param SrcMT src: The source of the fields to project
             :param SimPEG.Mesh mesh:
-            :param FieldsNSEM f: Natural source fields object to project
+            :param FieldsMT f: Natural source fields object to project
         '''
 
         ## NOTE: Assumes that e is on t
@@ -143,9 +143,9 @@ class Rx(SimPEGsurvey.BaseRx):
         """
         The derivative of the projection wrt u
 
-        :param NSEMsrc src: NSEM source
+        :param MTsrc src: MT source
         :param TensorMesh mesh: Mesh defining the topology of the problem
-        :param NSEMfields f: NSEM fields object of the source
+        :param MTfields f: MT fields object of the source
         :param numpy.ndarray v: Random vector of size
         """
 
@@ -390,12 +390,12 @@ class Rx(SimPEGsurvey.BaseRx):
 #################
 class Survey(SimPEGsurvey.BaseSurvey):
     """
-        Survey class for NSEM. Contains all the sources associated with the survey.
+        Survey class for MT. Contains all the sources associated with the survey.
 
         :param list srcList: List of sources associated with the survey
 
     """
-    srcPair = SrcNSEM.BaseNSEMSrc
+    srcPair = SrcMT.BaseMTSrc
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
@@ -443,7 +443,7 @@ class Survey(SimPEGsurvey.BaseSurvey):
 #################
 class Data(SimPEGsurvey.Data):
     '''
-    Data class for NSEMdata. Stores the data vector indexed by the survey.
+    Data class for MTdata. Stores the data vector indexed by the survey.
 
     :param SimPEG survey object survey:
     :param v vector of the data in order matching of the survey
@@ -461,7 +461,7 @@ class Data(SimPEGsurvey.Data):
 
     def toRecArray(self,returnType='RealImag'):
         '''
-        Function that returns a numpy.recarray for a SimpegNSEM impedance data object.
+        Function that returns a numpy.recarray for a SimpegMT impedance data object.
 
         :param str returnType: Switches between returning a rec array where the impedance is split to real and imaginary ('RealImag') or is a complex ('Complex')
 
@@ -483,7 +483,7 @@ class Data(SimPEGsurvey.Data):
                 locs = np.hstack((np.array([[0.0]]),locs))
             tArrRec = np.concatenate((src.freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],12))),axis=1).view(dtRI)
             # np.array([(src.freq,rx.locs[0,0],rx.locs[0,1],rx.locs[0,2],np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ) for rx in src.rxList],dtype=dtRI)
-            # Get the type and the value for the DataNSEM object as a list
+            # Get the type and the value for the DataMT object as a list
             typeList = [[rx.rxType.replace('z1d','zyx'),self[src,rx]] for rx in src.rxList]
             # Insert the values to the temp array
             for nr,(key,val) in enumerate(typeList):
@@ -517,17 +517,17 @@ class Data(SimPEGsurvey.Data):
     @classmethod
     def fromRecArray(cls, recArray, srcType='primary'):
         """
-        Class method that reads in a numpy record array to NSEMdata object.
+        Class method that reads in a numpy record array to MTdata object.
 
         Only imports the impedance data.
 
         """
         if srcType=='primary':
-            src = SrcNSEM.polxy_1Dprimary
+            src = SrcMT.polxy_1Dprimary
         elif srcType=='total':
-            src = SrcNSEM.polxy_1DhomotD
+            src = SrcMT.polxy_1DhomotD
         else:
-            raise NotImplementedError('{:s} is not a valid source type for NSEMdata')
+            raise NotImplementedError('{:s} is not a valid source type for MTdata')
 
         # Find all the frequencies in recArray
         uniFreq = np.unique(recArray['freq'])

SimPEG/MT/Utils/MT1Danalytic.py

@@ -3,7 +3,7 @@
 import numpy as np, SimPEG as simpeg
 from scipy.constants import mu_0, epsilon_0 as eps_0
 
-def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
+def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
     '''Analytic solution for MT 1D layered earth. Returns E and H fields.
 
     :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
@@ -12,7 +12,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
     :param numpy array, vector zd: location to calculate EH fields at
     :param bollean, scaleUD: scales the output to be 1 at the top, increases numeracal stability.
 
-    Assumes a halfspace with the same conductive as the deepest cell.
+    Assumes a halfspace with the same conductive as the last cell below.
 
     '''
     # Note add an error check for the mesh and sigma are the same size.
@@ -29,7 +29,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
 
     # Initiate the propagation matrix, in the order down up.
     UDp = np.zeros((2,m1d.nC+1),dtype=complex)
-    UDp[1,0] = scaleValue # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
+    UDp[1,0] = 1. # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
     # Loop over all the layers, starting at the bottom layer
     for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
         # Calculate
@@ -38,9 +38,9 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
         # Build the propagation matrix
 
         # Convert fields to down/up going components in layer below current layer
-        Pj1 = np.array([[1,1],[yp1,-yp1]],dtype=complex)
+        Pj1 = np.array([[1,1],[yp1,-yp1]])
         # Convert fields to down/up going components in current layer
-        Pjinv = 1./2*np.array([[1,zp],[1,-zp]],dtype=complex)
+        Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
         # Propagate down and up components through the current layer
         elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]])
 
@@ -48,14 +48,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
         UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
 
         if scaleUD:
-            # Scale the values such that 1 at the top
-            scaleVal = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
-            if np.any(np.isnan(scaleVal)):
-                # If there is a nan (thickness very great), rebuild the move up cell
-                scaleVal = np.zeros_like(UDp[:,lnr+1::-1],dtype=complex)
-                scaleVal[1,0] = scaleValue
-
-            UDp[:,lnr+1::-1] = scaleVal
+             UDp[:,lnr+1::-1] = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
 
     # Calculate the fields
     Ed = np.empty((zd.size,),dtype=complex)

SimPEG/MT/Utils/__init__.py

@@ -0,0 +1,4 @@
+from MT1Dsolutions import *  # Add the names of the functions
+from MT1Danalytic import *
+from dataUtils import *
+from ediFilesUtils import *

SimPEG/MT/Utils/dataUtils.py

@@ -5,25 +5,25 @@ import numpy.lib.recfunctions as recFunc
 from scipy.constants import mu_0
 from scipy import interpolate as sciint
 
-def getAppRes(NSEMdata):
+def getAppRes(MTdata):
     # Make impedance
     zList = []
-    for src in NSEMdata.survey.srcList:
+    for src in MTdata.survey.srcList:
         zc = [src.freq]
         for rx in src.rxList:
             if 'i' in rx.rxType:
                 m=1j
             else:
                 m = 1
-            zc.append(m*NSEMdata[src,rx])
+            zc.append(m*MTdata[src,rx])
         zList.append(zc)
     return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
 
-def rotateData(NSEMdata,rotAngle):
+def rotateData(MTdata, rotAngle):
     '''
     Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
     '''
-    recData = NSEMdata.toRecArray('Complex')
+    recData = MTdata.toRecArray('Complex')
     impData = rec2ndarr(recData[['zxx','zxy','zyx','zyy']],complex)
     # Make the rotation matrix
     # c,s,zxx,zxy,zyx,zyy = sympy.symbols('c,s,zxx,zxy,zyx,zyy')
@@ -40,27 +40,27 @@ def rotateData(NSEMdata,rotAngle):
     for nr,comp in enumerate(['zxx','zxy','zyx','zyy']):
         outRec[comp] = rotData[:,nr]
 
-    from SimPEG import NSEM
-    return NSEM.Data.fromRecArray(outRec)
+    from SimPEG import MT
+    return MT.Data.fromRecArray(outRec)
 
 
-def appResPhs(freq,z):
+def appResPhs(freq, z):
     app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
     app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
     return app_res, app_phs
 
-def skindepth(rho,freq):
+def skindepth(rho, freq):
     ''' Function to calculate the skindepth of EM waves'''
     return np.sqrt( (rho*((1/(freq * mu_0 * np.pi )))))
 
-def rec2ndarr(x,dt=float):
+def rec2ndarr(x, dt=float):
     return x.view((dt, len(x.dtype.names)))
 
-def makeAnalyticSolution(mesh,model,elev,freqs):
-    from SimPEG import NSEM
+def makeAnalyticSolution(mesh, model, elev, freqs):
+    from SimPEG import MT
     data1D = []
     for freq in freqs:
-        anaEd, anaEu, anaHd, anaHu = NSEM.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
+        anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
         anaE = anaEd+anaEu
         anaH = anaHd+anaHu
 
@@ -70,8 +70,8 @@ def makeAnalyticSolution(mesh,model,elev,freqs):
     dataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zyx',complex)])
     return dataRec
 
-def plotMT1DModelData(problem,models,symList=None):
-    from SimPEG import NSEM
+def plotMT1DModelData(problem, models, symList=None):
+    from SimPEG import MT
     # Setup the figure
     fontSize = 15
 
@@ -79,7 +79,7 @@ def plotMT1DModelData(problem,models,symList=None):
     axM = fig.add_axes([0.075,.1,.25,.875])
     axM.set_xlabel('Resistivity [Ohm*m]',fontsize=fontSize)
     axM.set_xlim(1e-1,1e5)
-    # axM.set_ylim(-10000,5000)
+    axM.set_ylim(-10000,5000)
     axM.set_ylabel('Depth [km]',fontsize=fontSize)
     axR = fig.add_axes([0.42,.575,.5,.4])
     axR.set_xscale('log')
@@ -132,94 +132,38 @@ def plotMT1DModelData(problem,models,symList=None):
     freq = simpeg.mkvc(data1D['freq'],2)
     res, phs = appResPhs(freq,allData)
 
-    if False:
-        stdCol = 'gray'
-        axRtw = axR.twinx()
-        axRtw.set_ylabel('Std of log10',color=stdCol)
-        [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
-        axPtw = axP.twinx()
-        axPtw.set_ylabel('Std ',color=stdCol)
-        [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
-        axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
-        axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
+    stdCol = 'gray'
+    axRtw = axR.twinx()
+    axRtw.set_ylabel('Std of log10',color=stdCol)
+    [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
+    axPtw = axP.twinx()
+    axPtw.set_ylabel('Std ',color=stdCol)
+    [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
+    axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
+    axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
 
     # Fix labels and ticks
 
-    # yMtick = [l/1000 for l in axM.get_yticks().tolist()]
-    # axM.set_yticklabels(yMtick)
+    yMtick = [l/1000 for l in axM.get_yticks().tolist()]
+    axM.set_yticklabels(yMtick)
     [ l.set_rotation(90) for l in axM.get_yticklabels()]
     [ l.set_rotation(90) for l in axR.get_yticklabels()]
-    # [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
-    # [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
+    [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
+    [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
     for ax in [axM,axR,axP]:
         ax.xaxis.set_tick_params(labelsize=fontSize)
         ax.yaxis.set_tick_params(labelsize=fontSize)
     return fig
 
-def plotImpAppRes(dataArrays,plotLoc,textStr=[]):
-    ''' Plots amplitude impedance and phase'''
-    # fig = plt.figure(1,(7, 7))
-    import plotDataTypes as pDt
-    # axes = ImageGrid(fig, (0.05,0.05,0.875,0.875),nrows_ncols = (2, 2),axes_pad = 0.25,add_all=True,share_all=True,label_mode = "L")
-    # Make the figure and axes
-    fig,axT=plt.subplots(2,2,sharex=True)
-    axes = axT.ravel()
-    fig.set_size_inches((13.5,7.0))
-    fig.suptitle('{:s}\nStation at: {:.1f}x ; {:.1f}y'.format(textStr,plotLoc[0],plotLoc[1]))
-    # Have to deal with axes
-    # Set log
-    for ax in axes.ravel():
-        ax.set_xscale('log')
-
-    axes[0].invert_xaxis()
-    axes[0].set_yscale('log')
-    axes[2].set_yscale('log')
-    # Set labels
-    axes[2].set_xlabel('Frequency [Hz]')
-    axes[3].set_xlabel('Frequency [Hz]')
-    axes[0].set_ylabel('Apperent resistivity [Ohm m]')
-    axes[1].set_ylabel('Apperent phase [degrees]')
-    axes[1].set_ylim(-180,180)
-    axes[2].set_ylabel('Impedance amplitude [V/A]')
-    axes[3].set_ylim(-180,180)
-    axes[3].set_ylabel('Impedance angle [degrees]')
-
-
-    # Plot the data
-    for nr,dataArray in enumerate(dataArrays):
-        if nr==1:
-            parSym = '*'
-        else:
-            parSym = 's'
-        # app res
-        pDt.plotIsoStaImpedance(axes[0],plotLoc,dataArray,'zxy',par='res',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[0],plotLoc,dataArray,'zyx',par='res',pSym=parSym)
-        # app phs
-        pDt.plotIsoStaImpedance(axes[1],plotLoc,dataArray,'zxy',par='phs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[1],plotLoc,dataArray,'zyx',par='phs',pSym=parSym)
-        # imp abs
-        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zxx',par='abs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zxy',par='abs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zyx',par='abs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zyy',par='abs',pSym=parSym)
-            # imp abs
-        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zxx',par='phs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zxy',par='phs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zyx',par='phs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zyy',par='phs',pSym=parSym)
-
-    return fig,axes
-
-
 def printTime():
     import time
     print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
 
-def convert3Dto1Dobject(NSEMdata,rxType3D='zyx'):
-    from SimPEG import NSEM
+def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
+    from SimPEG import MT
     # Find the unique locations
     # Need to find the locations
-    recDataTemp = NSEMdata.toRecArray()
+    recDataTemp = MTdata.toRecArray()
     # Check if survey.std has been assigned.
     ## NEED TO: write this...
     # Calculte and add the DET of the tensor to the recArray
@@ -241,24 +185,24 @@ def convert3Dto1Dobject(NSEMdata,rxType3D='zyx'):
         # Make the receiver list
         rx1DList = []
         for rxType in ['z1dr','z1di']:
-            rx1DList.append(NSEM.Rx(simpeg.mkvc(loc,2).T,rxType))
+            rx1DList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
         # Source list
         locrecData = recData[np.sqrt(np.sum( (rec2ndarr(recData[['x','y','z']]).data - loc )**2,axis=1)) < 1e-5]
         dat1DList = []
         src1DList = []
         for freq in locrecData['freq']:
-            src1DList.append(NSEM.SrcNSEM.src_polxy_1Dprimary(rx1DList,freq))
+            src1DList.append(MT.SrcMT.src_polxy_1Dprimary(rx1DList,freq))
             for comp  in ['r','i']:
                 dat1DList.append( corr * locrecData[rxType3D+comp][locrecData['freq']== freq].data )
 
         # Make the survey
-        sur1D = NSEM.Survey(src1DList)
+        sur1D = MT.Survey(src1DList)
 
         # Make the data
         dataVec = np.hstack(dat1DList)
-        dat1D = NSEM.Data(sur1D,dataVec)
+        dat1D = MT.Data(sur1D,dataVec)
         sur1D.dobs = dataVec
-        # Need to take NSEMdata.survey.std and split it as well.
+        # Need to take MTdata.survey.std and split it as well.
         std=0.05
         sur1D.std =  np.abs(sur1D.dobs*std) #+ 0.01*np.linalg.norm(sur1D.dobs)
         mtData1DList.append(dat1D)
@@ -266,29 +210,29 @@ def convert3Dto1Dobject(NSEMdata,rxType3D='zyx'):
     # Return the the list of data.
     return mtData1DList
 
-def resampleNSEMdataAtFreq(NSEMdata,freqs):
+def resampleMTdataAtFreq(MTdata,freqs):
     """
-    Function to resample NSEMdata at set of frequencies
+    Function to resample MTdata at set of frequencies
 
     """
-    from SimPEG import NSEM
+    from SimPEG import MT
     # Make a rec array
-    NSEMrec = NSEMdata.toRecArray().data
+    MTrec = MTdata.toRecArray().data
 
     # Find unique locations
-    uniLoc = np.unique(NSEMrec[['x','y','z']])
-    uniFreq = NSEMdata.survey.freqs
+    uniLoc = np.unique(MTrec[['x','y','z']])
+    uniFreq = MTdata.survey.freqs
     # Get the comps
-    dNames = NSEMrec.dtype
+    dNames = MTrec.dtype
 
     # Loop over all the locations and interpolate
     for loc in uniLoc:
         # Find the index of the station
-        ind = np.sqrt(np.sum((rec2ndarr(NSEMrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
+        ind = np.sqrt(np.sum((rec2ndarr(MTrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
         # Make a temporary recArray and interpolate all the components
         tArrRec = np.concatenate((simpeg.mkvc(freqs,2),np.ones((len(freqs),1))*rec2ndarr(loc),np.nan*np.ones((len(freqs),12))),axis=1).view(dNames)
         for comp in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            int1d = sciint.interp1d(NSEMrec[ind]['freq'],NSEMrec[ind][comp],bounds_error=False)
+            int1d = sciint.interp1d(MTrec[ind]['freq'],MTrec[ind][comp],bounds_error=False)
             tArrRec[comp] = simpeg.mkvc(int1d(freqs),2)
 
         # Join together
@@ -297,5 +241,5 @@ def resampleNSEMdataAtFreq(NSEMdata,freqs):
         except NameError as e:
             outRecArr = tArrRec
 
-    # Make the NSEMdata and return
-    return NSEM.Data.fromRecArray(outRecArr)
+    # Make the MTdata and return
+    return MT.Data.fromRecArray(outRecArr)

SimPEG/MT/Utils/ediFilesUtils.py

@@ -2,7 +2,7 @@
 from SimPEG import mkvc
 from scipy.constants import mu_0
 from numpy.lib import recfunctions as recFunc
-from SimPEG.NSEM.Utils.dataUtils import rec2ndarr
+from SimPEG.MT.Utils.dataUtils import rec2ndarr
 
 # Import modules
 import numpy as np

SimPEG/MT/Utils/sourceUtils.py

@@ -12,7 +12,7 @@ def homo1DModelSource(mesh,freq,sigma_1d):
 
     '''
     # import
-    from SimPEG.NSEM.Utils import get1DEfields
+    from SimPEG.MT.Utils import get1DEfields
     # Get a 1d solution for a halfspace background
     if mesh.dim == 1:
         mesh1d = mesh
@@ -77,7 +77,7 @@ def analytic1DModelSource(mesh,freq,sigma_1d):
 
     '''
     # import
-    from SimPEG.NSEM.Utils import getEHfields
+    from SimPEG.MT.Utils import getEHfields
     # Get a 1d solution for a halfspace background
     if mesh.dim == 1:
         mesh1d = mesh

SimPEG/MT/Utils/srcUtils.py

@@ -0,0 +1,46 @@
+import SimPEG as simpeg, numpy as np
+
+def homo1DModelSource(mesh,freq,m_back):
+    '''
+        Function that calculates and return background fields for a 3D mesh and model.
+        The calculuations use 1D field solution for a vertical slice throught model (south-western most column),
+        which is assigned at the fields everywhere for the respective polarizations.2
+
+        :param Simpeg mesh object mesh: Holds information on the discretization
+        :param float freq: The frequency to solve at
+        :param np.array m_back: Background model of conductivity to base the calculations on.
+        :rtype: numpy.ndarray (mesh.nE,2)
+        :return: eBG_bp, E fields for the background model at both polarizations.
+
+    '''
+
+    # import
+    from SimPEG.MT.Utils import get1DEfields
+    # Get a 1d solution for a halfspace background
+    mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
+    # Note: Everything is using e^iwt
+    e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq)
+    # Setup x (east) polarization (_x)
+    ex_px = np.zeros(mesh.vnEx,dtype=complex)
+    ey_px = np.zeros((mesh.nEy,1),dtype=complex)
+    ez_px = np.zeros((mesh.nEz,1),dtype=complex)
+    # Assign the source to ex_x
+    for i in np.arange(mesh.vnEx[0]):
+        for j in np.arange(mesh.vnEx[1]):
+            ex_px[i,j,:] = -e0_1d
+    eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
+    # Setup y (north) polarization (_py)
+    ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
+    ey_py = np.zeros(mesh.vnEy, dtype='complex128')
+    ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
+    # Assign the source to ey_py
+
+    for i in np.arange(mesh.vnEy[0]):
+        for j in np.arange(mesh.vnEy[1]):
+            ey_py[i,j,:] = e0_1d
+    # ey_py[1:-1,1:-1,1:-1] = 0
+    eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
+
+    # Return the electric fields
+    eBG_bp = np.hstack((eBG_px,eBG_py))
+    return eBG_bp

SimPEG/MT/__init__.py

@@ -0,0 +1,5 @@
+import Utils
+from SurveyMT import Rx, Survey, Data
+from FieldsMT import Fields1D_e, Fields3D_e
+import Problem1D, Problem2D, Problem3D
+import SrcMT

SimPEG/Maps.py

@@ -41,8 +41,8 @@ class IdentityMap(object):
             If this is a meshless mapping (i.e. nP is defined independently)
             the shape will be the the shape (nP,nP).
 
-            :rtype: (int,int)
-            :return: shape of the operator as a tuple
+            :rtype: tuple
+            :return: shape of the operator as a tuple (int,int)
         """
         if self._nP is not None:
             return (self.nP, self.nP)
@@ -86,7 +86,7 @@ class IdentityMap(object):
             The derivative of the transformation.
 
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
 
         """
@@ -216,7 +216,7 @@ class ExpMap(IdentityMap):
     def deriv(self, m):
         """
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
 
             The *transform* changes the model into the physical property.
@@ -366,7 +366,7 @@ class SurjectVertical1D(IdentityMap):
     def deriv(self, m):
         """
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
         """
         repNum = self.mesh.vnC[:self.mesh.dim-1].prod()
@@ -427,7 +427,7 @@ class Surject2Dto3D(IdentityMap):
     def deriv(self, m):
         """
             :param numpy.array m: model
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: derivative of transformed model
         """
         inds = self * np.arange(self.nP)
@@ -502,7 +502,9 @@ 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]

SimPEG/Mesh/BaseMesh.py

@@ -7,8 +7,8 @@ class BaseMesh(object):
     BaseMesh does all the counting you don't want to do.
     BaseMesh should be inherited by meshes with a regular structure.
 
-    :param numpy.array,list n: number of cells in each direction (dim, )
-    :param numpy.array,list x0: Origin of the mesh (dim, )
+    :param numpy.array n: (or list) number of cells in each direction (dim, )
+    :param numpy.array x0: (or list) Origin of the mesh (dim, )
 
     """
 
@@ -34,8 +34,8 @@ class BaseMesh(object):
         """
         Origin of the mesh
 
-        :rtype: numpy.array (dim, )
-        :return: x0
+        :rtype: numpy.array
+        :return: x0, (dim, )
         """
         return self._x0
 
@@ -116,8 +116,8 @@ class BaseMesh(object):
         """
         Total number of edges in each direction
 
-        :rtype: numpy.array (dim, )
-        :return: [nEx, nEy, nEz]
+        :rtype: numpy.array
+        :return: [nEx, nEy, nEz], (dim, )
 
         .. plot::
             :include-source:
@@ -173,8 +173,8 @@ class BaseMesh(object):
         """
         Total number of faces in each direction
 
-        :rtype: numpy.array (dim, )
-        :return: [nFx, nFy, nFz]
+        :rtype: numpy.array
+        :return: [nFx, nFy, nFz], (dim, )
 
         .. plot::
             :include-source:
@@ -200,8 +200,8 @@ class BaseMesh(object):
         """
         Face Normals
 
-        :rtype: numpy.array (sum(nF), dim)
-        :return: normals
+        :rtype: numpy.array
+        :return: normals, (sum(nF), dim)
         """
         if self.dim == 2:
             nX = np.c_[np.ones(self.nFx), np.zeros(self.nFx)]
@@ -218,8 +218,8 @@ class BaseMesh(object):
         """
         Edge Tangents
 
-        :rtype: numpy.array (sum(nE), dim)
-        :return: normals
+        :rtype: numpy.array
+        :return: normals, (sum(nE), dim)
         """
         if self.dim == 2:
             tX = np.c_[np.ones(self.nEx), np.zeros(self.nEx)]
@@ -236,8 +236,9 @@ class BaseMesh(object):
         Given a vector, fV, in cartesian coordinates, this will project it onto the mesh using the normals
 
         :param numpy.array fV: face vector with shape (nF, dim)
-        :rtype: numpy.array with shape (nF, )
-        :return: projected face vector
+        :rtype: numpy.array
+        :return: projected face vector, (nF, )
+
         """
         assert isinstance(fV, np.ndarray), 'fV must be an ndarray'
         assert len(fV.shape) == 2 and fV.shape[0] == self.nF and fV.shape[1] == self.dim, 'fV must be an ndarray of shape (nF x dim)'
@@ -248,8 +249,9 @@ class BaseMesh(object):
         Given a vector, eV, in cartesian coordinates, this will project it onto the mesh using the tangents
 
         :param numpy.array eV: edge vector with shape (nE, dim)
-        :rtype: numpy.array with shape (nE, )
-        :return: projected edge vector
+        :rtype: numpy.array
+        :return: projected edge vector, (nE, )
+
         """
         assert isinstance(eV, np.ndarray), 'eV must be an ndarray'
         assert len(eV.shape) == 2 and eV.shape[0] == self.nE and eV.shape[1] == self.dim, 'eV must be an ndarray of shape (nE x dim)'
@@ -295,7 +297,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Total number of cells in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: [nCx, nCy, nCz]
         """
         return np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None])
@@ -335,7 +337,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Total number of nodes in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: [nNx, nNy, nNz]
         """
         return np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None])
@@ -345,7 +347,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of x-edges in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnEx
         """
         return np.array([x for x in [self.nCx, self.nNy, self.nNz] if not x is None])
@@ -355,7 +357,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of y-edges in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnEy or None if dim < 2
         """
         return None if self.dim < 2 else np.array([x for x in [self.nNx, self.nCy, self.nNz] if not x is None])
@@ -365,7 +367,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of z-edges in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnEz or None if dim < 3
         """
         return None if self.dim < 3 else np.array([x for x in [self.nNx, self.nNy, self.nCz] if not x is None])
@@ -375,7 +377,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of x-faces in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnFx
         """
         return np.array([x for x in [self.nNx, self.nCy, self.nCz] if not x is None])
@@ -385,7 +387,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of y-faces in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnFy or None if dim < 2
         """
         return None if self.dim < 2 else np.array([x for x in [self.nCx, self.nNy, self.nCz] if not x is None])
@@ -395,7 +397,7 @@ class BaseRectangularMesh(BaseMesh):
         """
         Number of z-faces in each direction
 
-        :rtype: numpy.array (dim, )
+        :rtype: numpy.array
         :return: vnFz or None if dim < 3
         """
         return None if self.dim < 3 else np.array([x for x in [self.nCx, self.nCy, self.nNz] if not x is None])

SimPEG/Mesh/CylMesh.py

@@ -68,8 +68,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         """
         Number of x-faces in each direction
 
-        :rtype: numpy.array (dim, )
-        :return: vnFx
+        :rtype: numpy.array 
+        :return: vnFx, (dim, )
         """
         return self.vnC
 
@@ -78,8 +78,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         """
         Number of y-edges in each direction
 
-        :rtype: numpy.array (dim, )
-        :return: vnEy or None if dim < 2
+        :rtype: numpy.array 
+        :return: vnEy or None if dim < 2, (dim, )
         """
         nNx = self.nNx if self.isSymmetric else self.nNx - 1
         return np.r_[nNx, self.nCy, self.nNz]
@@ -89,8 +89,8 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         """
         Number of z-edges in each direction
 
-        :rtype: numpy.array (dim, )
-        :return: vnEz or None if nCy > 1
+        :rtype: numpy.array 
+        :return: vnEz or None if nCy > 1, (dim, )
         """
         if self.isSymmetric:
             return np.r_[self.nNx, self.nNy, self.nCz]

SimPEG/Mesh/InnerProducts.py

@@ -16,7 +16,7 @@ class InnerProducts(object):
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
             :param bool doFast: do a faster implementation if available.
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: M, the inner product matrix (nF, nF)
         """
         return self._getInnerProduct('F', prop=prop, invProp=invProp, invMat=invMat, doFast=doFast)
@@ -27,7 +27,7 @@ class InnerProducts(object):
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
             :param bool doFast: do a faster implementation if available.
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: M, the inner product matrix (nE, nE)
         """
         return self._getInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat, doFast=doFast)
@@ -39,7 +39,7 @@ class InnerProducts(object):
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
             :param bool doFast: do a faster implementation if available.
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: M, the inner product matrix (nE, nE)
         """
         assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"
@@ -115,13 +115,12 @@ class InnerProducts(object):
             :param bool doFast: do a faster implementation if available.
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
-            :rtype: function
             :return: dMdmu(u), the derivative of the inner product matrix (u)
 
             Given u, dMdmu returns (nF, nC*nA)
 
-            :param np.ndarray u: vector that multiplies dMdmu
-            :rtype: scipy.csr_matrix
+            :param numpy.ndarray u: vector that multiplies dMdmu
+            :rtype: scipy.sparse.csr_matrix
             :return: dMdmu, the derivative of the inner product matrix for a certain u
         """
         return self._getInnerProductDeriv(prop, 'F', doFast=doFast, invProp=invProp, invMat=invMat)
@@ -133,7 +132,7 @@ class InnerProducts(object):
             :param bool doFast: do a faster implementation if available.
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
         """
         return self._getInnerProductDeriv(prop, 'E', doFast=doFast, invProp=invProp, invMat=invMat)
@@ -145,7 +144,7 @@ class InnerProducts(object):
             :param bool doFast: do a faster implementation if available.
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (nE, nC*nA)
         """
         fast = None
@@ -169,7 +168,7 @@ class InnerProducts(object):
             :param numpy.array v: vector to multiply (required in the general implementation)
             :param list P: list of projection matrices
             :param str projType: 'F' for faces 'E' for edges
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: dMdm, the derivative of the inner product matrix (n, nC*nA)
         """
         assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"

SimPEG/Mesh/MeshIO.py

@@ -6,13 +6,11 @@ class TensorMeshIO(object):
     @classmethod
     def readUBC(TensorMesh, fileName):
         """
-            Read UBC GIF 3DTensor mesh and generate 3D Tensor mesh in simpegTD
+            Read UBC GIF 3D tensor mesh and generate 3D TensorMesh in SimPEG.
 
-            Input:
-            :param fileName, path to the UBC GIF mesh file
-
-            Output:
-            :param SimPEG TensorMesh object
+            :param string fileName: path to the UBC GIF mesh file
+            :rtype: TensorMesh
+            :return: The tensor mesh for the fileName.
         """
 
         # Interal function to read cell size lines for the UBC mesh files.
@@ -48,11 +46,9 @@ class TensorMeshIO(object):
             Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model
 
             Input:
-            :param vtrFileName, path to the vtr model file to write to
-
-            Output:
-            :return SimPEG TensorMesh object
-            :return SimPEG model dictionary
+            :param string fileName: path to the vtr model file to read
+            :rtype: tuple
+            :return: (TensorMesh, modelDictionary)
 
         """
         # Import
@@ -102,9 +98,8 @@ class TensorMeshIO(object):
         Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.
 
         Input:
-        :param str, path to the output vtk file
-        :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
-        :param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells
+        :param string fileName: path to the output vtk file
+        :param dict models: dictionary of numpy.array - Name('s) and array('s). Match number of cells
 
         """
         # Import
@@ -140,6 +135,7 @@ class TensorMeshIO(object):
                 vtkObj.GetCellData().AddArray(vtkDoubleArr)
             # Set the active scalar
             vtkObj.GetCellData().SetActiveScalars(models.keys()[0])
+        # vtkObj.Update()
 
         # Check the extension of the fileName
         ext = os.path.splitext(fileName)[1]
@@ -156,61 +152,14 @@ class TensorMeshIO(object):
         vtrWriteFilter.SetFileName(fileName)
         vtrWriteFilter.Update()
 
-    def _toVTRObj(mesh,models=None):
-        """
-        Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.
-
-        Input:
-        :param str, path to the output vtk file
-        :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
-        :param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells
-
-        """
-        # Import
-        from vtk import vtkRectilinearGrid as rectGrid, VTK_VERSION
-        from vtk.util.numpy_support import numpy_to_vtk
-
-        # Deal with dimensionalities
-        if mesh.dim >= 1:
-            vX = mesh.vectorNx
-            xD = mesh.nNx
-            yD,zD = 1,1
-            vY, vZ = np.array([0,0])
-        if mesh.dim >= 2:
-            vY = mesh.vectorNy
-            yD = mesh.nNy
-        if mesh.dim == 3:
-            vZ = mesh.vectorNz
-            zD = mesh.nNz
-        # Use rectilinear VTK grid.
-        # Assign the spatial information.
-        vtkObj = rectGrid()
-        vtkObj.SetDimensions(xD,yD,zD)
-        vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
-        vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
-        vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
-
-        # Assign the model('s) to the object
-        if models is not None:
-            for item in models.iteritems():
-                # Convert numpy array
-                vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
-                vtkDoubleArr.SetName(item[0])
-                vtkObj.GetCellData().AddArray(vtkDoubleArr)
-            # Set the active scalar
-            vtkObj.GetCellData().SetActiveScalars(models.keys()[0])
-        return vtkObj
 
     def readModelUBC(mesh, fileName):
         """
             Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg
 
-            Input:
-            :param fileName, path to the UBC GIF mesh file to read
-            :param mesh, TensorMesh object, mesh that coresponds to the model
-
-            Output:
-            :return numpy array, model with TensorMesh ordered
+            :param string fileName: path to the UBC GIF mesh file to read
+            :rtype: numpy.ndarray
+            :return: model with TensorMesh ordered
         """
         f = open(fileName, 'r')
         model = np.array(map(float, f.readlines()))
@@ -226,8 +175,7 @@ class TensorMeshIO(object):
             Writes a model associated with a SimPEG TensorMesh
             to a UBC-GIF format model file.
 
-            :param str fileName: File to write to
-            :param simpeg.Mesh.TensorMesh mesh: The mesh
+            :param string fileName: File to write to
             :param numpy.ndarray model: The model
         """
 
@@ -244,8 +192,8 @@ class TensorMeshIO(object):
         """
             Writes a SimPEG TensorMesh to a UBC-GIF format mesh file.
 
-            :param str fileName: File to write to
-            :param simpeg.Mesh.TensorMesh mesh: The mesh
+            :param string fileName: File to write to
+            :param dict models: A dictionary of the models
 
         """
         assert mesh.dim == 3
@@ -274,9 +222,8 @@ class TreeMeshIO(object):
         """
             Write UBC ocTree mesh and model files from a simpeg ocTree mesh and model.
 
-            :param str fileName: File to write to
-            :param simpeg.Mesh.TreeMesh mesh: The mesh
-            :param dictionary models: The models in a dictionary, where the keys is the name of the of the model file
+            :param string fileName: File to write to
+            :param dict models: The models in a dictionary, where the keys is the name of the of the model file
         """
 
         # Calculate information to write in the file.
@@ -329,10 +276,9 @@ class TreeMeshIO(object):
 
             Input:
             :param str meshFile: path to the UBC GIF OcTree mesh file to read
+            :rtype: SimPEG.Mesh.TreeMesh
+            :return: The octree mesh
 
-            Output:
-            :return SimPEG.Mesh.TreeMesh mesh: The octree mesh
-            :return list of ndarray's: models as a list of numpy array's
         """
 
         ## Read the file lines
@@ -378,11 +324,9 @@ class TreeMeshIO(object):
         """
             Read UBC OcTree model and get vector
 
-            Input:
-            :param fileName, path to the UBC GIF model file to read
-
-            Output:
-            :return numpy array, OcTree model
+            :param string fileName: path to the UBC GIF model file to read
+            :rtype: numpy.ndarray
+            :return: OcTree model
         """
 
         if type(fileName) is list:

SimPEG/Mesh/TensorMesh.py

@@ -198,8 +198,8 @@ class BaseTensorMesh(BaseMesh):
         Determines if a set of points are inside a mesh.
 
         :param numpy.ndarray pts: Location of points to test
-        :rtype numpy.ndarray
-        :return inside, numpy array of booleans
+        :rtype numpy.ndarray:
+        :return: inside, numpy array of booleans
         """
         pts = Utils.asArray_N_x_Dim(pts, self.dim)
 
@@ -221,7 +221,7 @@ class BaseTensorMesh(BaseMesh):
 
         :param numpy.ndarray loc: Location of points to interpolate to
         :param str locType: What to interpolate (see below)
-        :rtype: scipy.sparse.csr.csr_matrix
+        :rtype: scipy.sparse.csr_matrix
         :return: M, the interpolation matrix
 
         locType can be::
@@ -289,7 +289,7 @@ class BaseTensorMesh(BaseMesh):
             :param bool returnP: returns the projection matrices
             :param bool invProp: inverts the material property
             :param bool invMat: inverts the matrix
-            :rtype: scipy.csr_matrix
+            :rtype: scipy.sparse.csr_matrix
             :return: M, the inner product matrix (nF, nF)
         """
         assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges"

SimPEG/Mesh/TreeMesh.py

@@ -1875,7 +1875,7 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
 
         :param numpy.ndarray locs: Location of points to interpolate to
         :param str locType: What to interpolate (see below)
-        :rtype: scipy.sparse.csr.csr_matrix
+        :rtype: scipy.sparse.csr_matrix
         :return: M, the interpolation matrix
 
         locType can be::

SimPEG/NSEM/ProblemNSEM.py

@@ -1,560 +0,0 @@
-from SimPEG.EM.Utils.EMUtils import omega, mu_0
-from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
-from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
-from SurveyNSEM import Survey, Data
-from FieldsNSEM import BaseNSEMFields, Fields1D_ePrimSec, Fields3D_ePrimSec
-from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields
-import time, sys
-
-class BaseNSEMProblem(BaseFDEMProblem):
-    """
-        Base class for all Natural source problems.
-    """
-
-    def __init__(self, mesh, **kwargs):
-        BaseFDEMProblem.__init__(self, mesh, **kwargs)
-        Utils.setKwargs(self, **kwargs)
-    # Set the default pairs of the problem
-    surveyPair = Survey
-    dataPair = Data
-    fieldsPair = BaseNSEMFields
-
-    # Set the solver
-    Solver = SimpegSolver
-    solverOpts = {}
-
-    verbose = False
-    # Notes:
-    # Use the forward and devs from BaseFDEMProblem
-    # Might need to add more stuff here.
-
-    ## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
-    def Jvec(self, m, v, f=None):
-        """
-        Function to calculate the data sensitivities dD/dm times a vector.
-
-            :param numpy.ndarray m (nC, 1) - conductive model
-            :param numpy.ndarray v (nC, 1) - random vector
-            :param NSEMfields object (optional) - NSEM fields object, if not given it is calculated
-            :rtype: NSEMdata object
-            :return: Data sensitivities wrt m
-        """
-
-        # Calculate the fields
-        if f is None:
-           f= self.fields(m)
-        # Set current model
-        self.curModel = m
-        # Initiate the Jv object
-        Jv = self.dataPair(self.survey)
-
-        # Loop all the frequenies
-        for freq in self.survey.freqs:
-            dA_du = self.getA(freq) #
-
-            dA_duI = self.Solver(dA_du, **self.solverOpts)
-
-            for src in self.survey.getSrcByFreq(freq):
-                # We need fDeriv_m = df/du*du/dm + df/dm
-                # Construct du/dm, it requires a solve
-                # NOTE: need to account for the 2 polarizations in the derivatives.
-                u_src = f[src,:] # u should be a vector by definition. Need to fix this...
-                # dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
-                dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
-                dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
-                if dRHS_dm is None:
-                    du_dm = dA_duI * ( -dA_dm )
-                else:
-                    du_dm = dA_duI * ( -dA_dm + dRHS_dm )
-                # Calculate the projection derivatives
-                for rx in src.rxList:
-                    # Get the projection derivative
-                    # v should be of size 2*nE (for 2 polarizations)
-                    PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
-                    Jv[src, rx] = PDeriv_u(mkvc(du_dm))
-            dA_duI.clean()
-        # Return the vectorized sensitivities
-        return mkvc(Jv)
-
-    def Jtvec(self, m, v, f=None):
-        """
-        Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
-
-            :param numpy.ndarray m (nC, 1) - conductive model
-            :param numpy.ndarray v (nD, 1) - vector
-            :param NSEMfields object f (optional) - NSEM fields object, if not given it is calculated
-            :rtype: NSEMdata object
-            :return: Data sensitivities wrt m
-        """
-
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-
-        # Ensure v is a data object.
-        if not isinstance(v, self.dataPair):
-            v = self.dataPair(self.survey, v)
-
-        Jtv = np.zeros(m.size)
-
-        for freq in self.survey.freqs:
-            AT = self.getA(freq).T
-
-            ATinv = self.Solver(AT, **self.solverOpts)
-
-            for src in self.survey.getSrcByFreq(freq):
-                ftype = self._solutionType
-                f_src = f[src, :] # Need to fix this...
-
-                for rx in src.rxList:
-                    # Get the adjoint evalDeriv
-                    # PTv needs to be nE,
-                    PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
-                    # Get the
-                    dA_duIT = ATinv * PTv
-                    dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
-                    dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
-                    # Make du_dmT
-                    if dRHS_dmT is None:
-                        du_dmT = -dA_dmT
-                    else:
-                        du_dmT = -dA_dmT + dRHS_dmT
-                    # Select the correct component
-                    # du_dmT needs to be of size nC,
-                    real_or_imag = rx.projComp
-                    if real_or_imag == 'real':
-                        Jtv +=  du_dmT.real
-                    elif real_or_imag == 'imag':
-                        Jtv +=  -du_dmT.real
-                    else:
-                        raise Exception('Must be real or imag')
-            # Clean the factorization, clear memory.
-            ATinv.clean()
-        return Jtv
-
-###################################
-## 1D problems
-###################################
-
-class Problem1D_ePrimSec(BaseNSEMProblem):
-    """
-    A NSEM problem soving a e formulation and primary/secondary fields decomposion.
-
-    By eliminating the magnetic flux density using
-
-        .. math ::
-
-            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
-
-
-    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
-
-    .. math ::
-        \\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
-    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
-
-    The primary field is estimated from a background model (commonly half space ).
-
-
-    """
-
-    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
-    _solutionType = 'e_1dSolution'
-    _formulation  = 'EF'
-    fieldsPair = Fields1D_ePrimSec
-
-    # Initiate properties
-    _sigmaPrimary = None
-
-
-    def __init__(self, mesh, **kwargs):
-        BaseNSEMProblem.__init__(self, mesh, **kwargs)
-        # self._sigmaPrimary = sigmaPrimary
-    @property
-    def MeMui(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MeMui', None) is None:
-            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
-        return self._MeMui
-
-    @property
-    def MfSigma(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MfSigma', None) is None:
-            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
-        return self._MfSigma
-
-    @property
-    def sigmaPrimary(self):
-        """
-        A background model, use for the calculation of the primary fields.
-
-        """
-        return self._sigmaPrimary
-
-    @sigmaPrimary.setter
-    def sigmaPrimary(self, val):
-        # Note: TODO add logic for val, make sure it is the correct size.
-        self._sigmaPrimary = val
-
-    def getA(self, freq):
-        """
-            Function to get the A matrix.
-
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-
-        # Note: need to use the code above since in the 1D problem I want
-        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
-        # Possible that _fieldType and _eqLocs can fix this
-        MeMui = self.MeMui
-        MfSigma = self.MfSigma
-        C = self.mesh.nodalGrad
-        # Make A
-        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
-        # Either return full or only the inner part of A
-        return A
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        """
-        The derivative of A wrt sigma
-        """
-
-        dsig_dm = self.curModel.sigmaDeriv
-        MeMui = self.MeMui
-        #
-        u_src = u['e_1dSolution']
-        dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
-        if adjoint:
-            return 1j * omega(freq) * (  dMfSigma_dm.T * v )
-        # Note: output has to be nN/nF, not nC/nE.
-        # v should be nC
-        return 1j * omega(freq) * ( dMfSigma_dm * v )
-
-    def getRHS(self, freq):
-        """
-            Function to return the right hand side for the system.
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
-            :return: RHS for 1 polarizations, primary fields
-        """
-
-        # Get sources for the frequncy(polarizations)
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_e = Src.S_e(self)
-        return -1j * omega(freq) * S_e
-
-    def getRHSDeriv_m(self, freq, v, adjoint=False):
-        """
-        The derivative of the RHS wrt sigma
-        """
-
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
-        return -1j * omega(freq) * S_eDeriv
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        '''
-        # Set the current model
-        self.curModel = m
-        # Make the fields object
-        F = self.fieldsPair(self.mesh, self.survey)
-        # Loop over the frequencies
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs  = self.getRHS(freq)
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_s = Ainv * rhs
-
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)[0]
-            # NOTE: only store the e_solution(secondary), all other components calculated in the fields object
-            F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
-
-            # Note curl e = -iwb so b = -curl e /iw
-            # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
-            # F[Src, 'b_1d'] = b[:,1]
-            if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
-                sys.stdout.flush()
-        return F
-
-# Note this is not fully functional.
-# Missing:
-# Fields class corresponding to the fields
-# Update Jvec and Jtvec to include all the derivatives components
-# Other things ...
-class Problem1D_eTotal(BaseNSEMProblem):
-    """
-    A NSEM problem solving a e formulation and a Total bondary domain decompostion.
-
-    Solves the equation:
-
-    Math:
-        Have to do this...
-        Not implement correctly.......
-    """
-
-    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
-    _solutionType = 'e_1dSolution'
-    _formulation  = 'EF'
-    # fieldsPair = Fields1D_eTotal
-
-    def __init__(self, mesh, **kwargs):
-        BaseNSEMProblem.__init__(self, mesh, **kwargs)
-    @property
-    def MeMui(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MeMui', None) is None:
-            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
-        return self._MeMui
-
-    @property
-    def MfSigma(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MfSigma', None) is None:
-            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
-        return self._MfSigma
-
-    def getA(self, freq, full=False):
-        """
-            Function to get the A matrix.
-
-            :param float freq: Frequency
-            :param logic full: Return full A or the inner part
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-
-        MeMui = self.MeMui
-        MfSigma = self.MfSigma
-        # Note: need to use the code above since in the 1D problem I want
-        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
-        # Possible that _fieldType and _eqLocs can fix this
-        # MeMui = self.MfMui
-        # MfSigma = self.MfSigma
-        C = self.mesh.nodalGrad
-        # Make A
-        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
-        # Either return full or only the inner part of A
-        if full:
-            return A
-        else:
-            return A[1:-1,1:-1]
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        raise NotImplementedError('getADeriv is not implemented')
-
-    def getRHS(self, freq):
-        """
-            Function to return the right hand side for the system.
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
-            :return: RHS for both polarizations, primary fields
-        """
-        # Get sources for the frequency
-        # NOTE: Need to use the source information, doesn't really apply in 1D
-        src = self.survey.getSrcByFreq(freq)
-        # Get the full A
-        A = self.getA(freq,full=True)
-        # Define the outer part of the solution matrix
-        Aio = A[1:-1,[0,-1]]
-        Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
-        Etot = (Ed + Eu)
-        sourceAmp = 1.0
-        Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
-        ## Note: The analytic solution is derived with e^iwt
-        eBC = np.r_[Etot[0],Etot[-1]]
-        # The right hand side
-
-        return -Aio*eBC, eBC
-
-    def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
-        raise NotImplementedError('getRHSDeriv not implemented yet')
-        return None
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        :param np.ndarray (nC,) m_back: Background conductivity model
-        '''
-
-
-        self.curModel = m
-        # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
-
-        F = Fields1D_eTotal(self.mesh, self.survey)
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs, e_o = self.getRHS(freq)
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_i = Ainv * rhs
-            e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)
-            # NOTE: only store e fields
-            F[Src, 'e_1dSolution'] = e[:,0]
-            if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
-                sys.stdout.flush()
-        return F
-
-
-###################################
-## 3D problems
-###################################
-class Problem3D_ePrimSec(BaseNSEMProblem):
-    """
-    A NSEM problem solving a e formulation and a primary/secondary fields decompostion.
-
-    By eliminating the magnetic flux density using
-
-        .. math ::
-
-            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
-
-
-    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
-
-    .. math ::
-        \\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
-    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
-
-    The primary field is estimated from a background model (commonly as a 1D model).
-
-    """
-
-    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
-    _solutionType = [ 'e_pxSolution', 'e_pySolution']  # Forces order on the object
-    _formulation  = 'EB'
-    fieldsPair = Fields3D_ePrimSec
-
-    # Initiate properties
-    _sigmaPrimary = None
-
-    def __init__(self, mesh, **kwargs):
-        BaseNSEMProblem.__init__(self, mesh, **kwargs)
-
-    @property
-    def sigmaPrimary(self):
-        """
-        A background model, use for the calculation of the primary fields.
-
-        """
-        return self._sigmaPrimary
-    @sigmaPrimary.setter
-    def sigmaPrimary(self, val):
-        # Note: TODO add logic for val, make sure it is the correct size.
-        self._sigmaPrimary = val
-
-    def getA(self, freq):
-        """
-            Function to get the A system.
-
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-        Mmui = self.MfMui
-        Msig = self.MeSigma
-        C = self.mesh.edgeCurl
-
-        return C.T*Mmui*C + 1j*omega(freq)*Msig
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        """
-        Calculate the derivative of A wrt m.
-
-        """
-        # Fix u to be a matrix nE,2
-        # This considers both polarizations and returns a nE,2 matrix for each polarization
-        if adjoint:
-            dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v
-        else:
-            # Need a nE,2 matrix to be returned
-            dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
-        return 1j * omega(freq) * dMe_dsigV
-
-
-    def getRHS(self, freq):
-        """
-        Function to return the right hand side for the system.
-
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
-            :return: RHS for both polarizations, primary fields
-        """
-
-        # Get sources for the frequncy(polarizations)
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_e = Src.S_e(self)
-        return -1j * omega(freq) * S_e
-
-    def getRHSDeriv_m(self, freq, v, adjoint=False):
-        """
-        The derivative of the RHS with respect to sigma
-        """
-
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
-        return -1j * omega(freq) * S_eDeriv
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        '''
-        # Set the current model
-        self.curModel = m
-
-        F = self.fieldsPair(self.mesh, self.survey)
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs  = self.getRHS(freq)
-            # Solve the system
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_s = Ainv * rhs
-
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)[0]
-            # Store the fields
-            # Use self._solutionType
-            F[Src, 'e_pxSolution'] = e_s[:,0]
-            F[Src, 'e_pySolution'] = e_s[:,1]
-            # Note curl e = -iwb so b = -curl/iw
-
-            if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
-                sys.stdout.flush()
-            Ainv.clean()
-        return F

SimPEG/NSEM/Utils/__init__.py

@@ -1,5 +0,0 @@
-from MT1Dsolutions import get1DEfields  # Add the names of the functions
-from MT1Danalytic import getEHfields, getImpedance
-from dataUtils import *
-from ediFilesUtils import *
-from testUtils import *

SimPEG/NSEM/Utils/testUtils.py

@@ -1,198 +0,0 @@
-import unittest
-import sys
-from scipy.constants import mu_0
-import SimPEG as simpeg
-
-from SimPEG.Utils import meshTensor
-import numpy as np
-
-np.random.seed(1100)
-# Define the tolerances
-TOLr = 5e-2
-TOLp = 5e-2
-
-
-def getAppResPhs(NSEMdata):
-    # Make impedance
-    from SimPEG.NSEM.Utils import appResPhs
-    zList = []
-    for src in NSEMdata.survey.srcList:
-        zc = [src.freq]
-        for rx in src.rxList:
-            if 'i' in rx.rxType:
-                m=1j
-            else:
-                m = 1
-            zc.append(m*NSEMdata[src,rx])
-        zList.append(zc)
-    return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
-
-
-def setup1DSurvey(sigmaHalf,tD=True,structure=False):
-    from SimPEG import NSEM
-    # Frequency
-    nFreq = 33
-    freqs = np.logspace(3,-3,nFreq)
-    # Make the mesh
-    ct = 5
-    air = meshTensor([(ct,25,1.3)])
-    # coreT0 = meshTensor([(ct,15,1.2)])
-    # coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
-    core = np.concatenate( (  np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
-    bot = meshTensor([(core[0],20,-1.3)])
-    x0 = -np.array([np.sum(np.concatenate((core,bot)))])
-    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
-    # Make the model
-    sigma = np.zeros(m1d.nC) + sigmaHalf
-    sigma[m1d.gridCC > 0 ] = 1e-8
-    sigmaBack = sigma.copy()
-    # Add structure
-    if structure:
-        shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600)
-        deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000)
-        sigma[shallow] = 1
-        sigma[deep] = 0.1
-
-    rxList = []
-    for rxType in ['z1dr','z1di']:
-        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
-    # Source list
-    srcList =[]
-    if tD:
-        for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq))
-    else:
-        for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
-
-    survey = NSEM.Survey(srcList)
-    return survey, sigma, m1d
-
-
-def setupSimpegNSEM_ePrimSec(inputSetup,comp='Imp',singleFreq=False,expMap=True):
-    from SimPEG import NSEM
-
-    M,freqs,sig,sigBG,rx_loc = inputSetup
-    # Make a receiver list
-    rxList = []
-    if comp == 'All':
-        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(NSEM.Rx(rx_loc,rxType))
-    elif comp == 'Imp':
-        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
-            rxList.append(NSEM.Rx(rx_loc,rxType))
-    elif comp == 'Tip':
-        for rxType in ['tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(NSEM.Rx(rx_loc,rxType))
-    else:
-        rxList.append(NSEM.Rx(rx_loc,comp))
-    # Source list
-    srcList =[]
-
-    if singleFreq:
-        srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,singleFreq))
-    else:
-        for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
-    # Survey NSEM
-    survey = NSEM.Survey(srcList)
-
-    ## Setup the problem object
-    sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
-    if expMap:
-        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= np.log(sigma1d) )
-        problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
-        problem.curModel = np.log(sig)
-    else:
-        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= sigma1d)
-        problem.curModel = sig
-    problem.pair(survey)
-    problem.verbose = False
-    try:
-        from pymatsolver import MumpsSolver
-        problem.Solver = MumpsSolver
-    except:
-        pass
-
-    return (survey, problem)
-
-def getInputs():
-    """
-    Function that returns Mesh, freqs, rx_loc, elev.
-    """
-    # Make a mesh
-    # M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360])
-    # M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model
-    M = simpeg.Mesh.TensorMesh([[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,8,-1.5),(200.,8),(200,8,1.5)]], x0=['C','C','C'])# Setup the model
-    # Set the frequencies
-    freqs = np.logspace(1,-3,5)
-    elev = 0
-
-    ## Setup the the survey object
-    # Receiver locations
-    rx_x, rx_y = np.meshgrid(np.arange(-350,350,200),np.arange(-350,350,200))
-    rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
-
-    return M, freqs, rx_loc, elev
-
-def random(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Backround
-    sigBG = np.ones(M.nC)*conds
-    # Add randomness to the model (10% of the value).
-    sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 )
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def halfSpace(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    # conds = [1e-2]
-    groundInd = ccM[:,2] < elev
-    sig = np.zeros(M.nC) + 1e-8
-    sig[groundInd] = conds
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def blockInhalfSpace(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    # conds = [1e-2]
-    groundInd = ccM[:,2] < elev
-    sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds)
-    sig[~groundInd] = 1e-8
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds[1]
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def twoLayer(conds):
-    ''' Returns a 2 layer model based on the conductivity values given'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    groundInd = ccM[:,2] < elev
-    botInd = ccM[:,2] < -3000
-    sig = np.zeros(M.nC) + 1e-8
-    sig[groundInd] = conds[1]
-    sig[botInd] = conds[0]
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds[1]
-
-
-    return (M, freqs, sig, sigBG, rx_loc)
-

SimPEG/NSEM/__init__.py

@@ -1,5 +0,0 @@
-import Utils
-from SurveyNSEM import Rx, Survey, Data
-from FieldsNSEM import Fields1D_ePrimSec, Fields3D_ePrimSec
-from ProblemNSEM import Problem1D_ePrimSec, Problem3D_ePrimSec
-import SrcNSEM

SimPEG/Optimization.py

@@ -131,7 +131,7 @@ class Minimize(object):
 
         Minimizes the function (evalFunction) starting at the location x0.
 
-        :param def evalFunction: function handle that evaluates: f, g, H = F(x)
+        :param callable evalFunction: function handle that evaluates: f, g, H = F(x)
         :param numpy.ndarray x0: starting location
         :rtype: numpy.ndarray
         :return: x, the last iterate of the optimization algorithm
@@ -372,8 +372,8 @@ class Minimize(object):
             Else, a modifySearchDirectionBreak call is preformed.
 
             :param numpy.ndarray p: searchDirection
-            :rtype: numpy.ndarray,bool
-            :return: (xt, passLS)
+            :rtype: tuple
+            :return: (xt, passLS) numpy.ndarray, bool
         """
         # Projected Armijo linesearch
         self._LS_t = 1
@@ -408,8 +408,8 @@ class Minimize(object):
             evalFunction returns a False indicating the break was not caught.
 
             :param numpy.ndarray p: searchDirection
-            :rtype: numpy.ndarray,bool
-            :return: (xt, breakCaught)
+            :rtype: tuple
+            :return: (xt, breakCaught) numpy.ndarray, bool
         """
         self.printDone(inLS=True)
         print 'The linesearch got broken. Boo.'

SimPEG/PropMaps.py

@@ -187,7 +187,7 @@ class _PropMapMetaClass(type):
             attrs[attr + 'Model'] = prop._getModelProperty()
             attrs[attr + 'Deriv'] = prop._getModelDerivProperty()
 
-        return type(name.replace('PropMap', 'PropModel'), (PropModel, ), attrs)
+        return type('PropModel', (PropModel, ), attrs)
 
 
 class PropMap(object):

SimPEG/Regularization.py

@@ -10,7 +10,7 @@ class RegularizationMesh(object):
     are not necessarily true differential operators, but are constructed from
     a SimPEG Mesh.
 
-    :param Mesh mesh: problem mesh
+    :param BaseMesh mesh: problem mesh
     :param numpy.array indActive: bool array, size nC, that is True where we have active cells. Used to reduce the operators so we regularize only on active cells
     """
 
@@ -383,8 +383,8 @@ class BaseRegularization(object):
 
         :param numpy.array m: geophysical model
         :param numpy.array v: vector to multiply
-        :rtype: scipy.sparse.csr_matrix or numpy.ndarray
-        :return: WtW or WtW*v
+        :rtype: scipy.sparse.csr_matrix
+        :return: WtW, or if v is supplied WtW*v (numpy.ndarray)
 
         The regularization is:
 
@@ -650,8 +650,8 @@ class Tikhonov(Simple):
     Note if the key word argument `mrefInSmooth` is False, then mref is not
     included in the smoothness contribution.
 
-    :param Mesh mesh: SimPEG mesh
-    :param Maps mapping: regularization mapping, takes the model from model space to the thing you want to regularize
+    :param BaseMesh mesh: SimPEG mesh
+    :param IdentityMap mapping: regularization mapping, takes the model from model space to the thing you want to regularize
     :param numpy.ndarray indActive: active cell indices for reducing the size of differential operators in the definition of a regularization mesh
     :param bool mrefInSmooth: (default = False) put mref in the smoothness component?
     :param float alpha_s: (default 1e-6) smallness weight
@@ -671,7 +671,7 @@ class Tikhonov(Simple):
     alpha_yy     = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'],    "Weight for the second derivative in the y direction")
     alpha_zz     = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'],    "Weight for the second derivative in the z direction")
 
-    def __init__(self, mesh, mapping=None, indActive = None, **kwargs):
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
         BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
 
     @property

SimPEG/Survey.py

@@ -311,7 +311,6 @@ class BaseSurvey(object):
         if f is None: f = self.prob.fields(m)
         return Utils.mkvc(self.eval(f))
 
-
     @Utils.count
     def eval(self, f):
         """eval(f)
@@ -322,7 +321,7 @@ class BaseSurvey(object):
 
                 d_\\text{pred} = \mathbf{P} f(m)
         """
-        raise NotImplemented('eval is not yet implemented.')
+        raise NotImplementedError('eval is not yet implemented.')
 
     @Utils.count
     def evalDeriv(self, f):
@@ -334,7 +333,7 @@ class BaseSurvey(object):
 
                 \\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P}
         """
-        raise NotImplemented('eval is not yet implemented.')
+        raise NotImplementedError('eval is not yet implemented.')
 
     @Utils.count
     def residual(self, m, f=None):

SimPEG/Tests.py

@@ -237,7 +237,7 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
         Compares error decay of 0th and 1st order Taylor approximation at point
         x0 for a randomized search direction.
 
-        :param lambda fctn: function handle
+        :param callable fctn: function handle
         :param numpy.array x0: point at which to check derivative
         :param int num: number of times to reduce step length, h
         :param bool plotIt: if you would like to plot

SimPEG/Utils/ModelBuilder.py

@@ -7,11 +7,11 @@ def addBlock(gridCC, modelCC, p0, p1, blockProp):
     """
         Add a block to an exsisting cell centered model, modelCC
 
-        :param numpy.array, gridCC: mesh.gridCC is the cell centered grid
-        :param numpy.array, modelCC: cell centered model
-        :param numpy.array, p0: bottom, southwest corner of block
-        :param numpy.array, p1: top, northeast corner of block
-        :blockProp float, blockProp: property to assign to the model 
+        :param numpy.array gridCC: mesh.gridCC is the cell centered grid
+        :param numpy.array modelCC: cell centered model
+        :param numpy.array p0: bottom, southwest corner of block
+        :param numpy.array p1: top, northeast corner of block
+        :blockProp float blockProp: property to assign to the model
 
         :return numpy.array, modelBlock: model with block
     """
@@ -147,7 +147,7 @@ def getIndicesSphere(center,radius,ccMesh):
 
     if dimMesh == 1:
        # Define the reference points
-        
+
         ind  = np.abs(center[0] - ccMesh[:,0]) < radius
 
     elif dimMesh == 2:
@@ -222,14 +222,14 @@ def layeredModel(ccMesh, layerTops, layerValues):
 
         :param numpy.array ccMesh: cell-centered mesh
         :param numpy.array layerTops: z-locations of the tops of each layer
-        :param numpy.array layerValue: values of the property to assign for each layer (starting at the top) 
+        :param numpy.array layerValue: values of the property to assign for each layer (starting at the top)
         :rtype: numpy.array
-        :return: M, layered model on the mesh 
+        :return: M, layered model on the mesh
     """
 
     descending = np.linalg.norm(sorted(layerTops, reverse=True) - layerTops) < 1e-20
 
-    # TODO: put an error check to make sure that there is an ordering... needs to work with inf elts 
+    # TODO: put an error check to make sure that there is an ordering... needs to work with inf elts
     # assert ascending or descending, "Layers must be listed in either ascending or descending order"
 
     # start from bottom up
@@ -253,10 +253,10 @@ def layeredModel(ccMesh, layerTops, layerValues):
     model = np.zeros(ccMesh.shape[0])
 
     for i, top in enumerate(layerTops):
-        zind = z <= top 
+        zind = z <= top
         model[zind] = layerValues[i]
 
-    return model 
+    return model
 
 
 
@@ -265,9 +265,9 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=None):
         Create a random model by convolving a kernel with a
         uniformly distributed model.
 
-        :param int,tuple shape: shape of the model.
+        :param tuple shape: shape of the model.
         :param int seed: pick which model to produce, prints the seed if you don't choose.
-        :param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernel that is used.
+        :param numpy.ndarray anisotropy: this is the (3 x n) blurring kernel that is used.
         :param int its: number of smoothing iterations
         :param list bounds: bounds on the model, len(list) == 2
         :rtype: numpy.ndarray

SimPEG/Utils/SolverUtils.py

@@ -13,7 +13,7 @@ def _checkAccuracy(A, b, X, accuracyTol):
         warnings.warn(msg, RuntimeWarning)
 
 
-def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
+def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6, name=None):
     """
     Wraps a direct Solver.
 
@@ -72,11 +72,11 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
         if factorize and hasattr(self.solver, 'clean'):
             return self.solver.clean()
 
-    return type(fun.__name__+'_Wrapped', (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
+    return type(name if name is not None else fun.__name__, (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
 
 
 
-def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5):
+def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5, name=None):
     """
     Wraps an iterative Solver.
 
@@ -128,13 +128,13 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5):
     def clean(self):
         pass
 
-    return type(fun.__name__+'_Wrapped', (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
+    return type(name if name is not None else fun.__name__, (object,), {"__init__": __init__, "clean": clean, "__mul__": __mul__})
 
 
 from scipy.sparse import linalg
-Solver   = SolverWrapD(linalg.spsolve, factorize=False)
-SolverLU = SolverWrapD(linalg.splu, factorize=True)
-SolverCG = SolverWrapI(linalg.cg)
+Solver   = SolverWrapD(linalg.spsolve, factorize=False, name="Solver")
+SolverLU = SolverWrapD(linalg.splu, factorize=True, name="SolverLU")
+SolverCG = SolverWrapI(linalg.cg, name="SolverCG")
 
 
 class SolverDiag(object):

SimPEG/Utils/interputils.py

@@ -25,7 +25,7 @@ def interpmat(locs, x, y=None, z=None):
         :param numpy.ndarray x: Tensor vector of 1st dimension of grid.
         :param numpy.ndarray y: Tensor vector of 2nd dimension of grid. None by default.
         :param numpy.ndarray z: Tensor vector of 3rd dimension of grid. None by default.
-        :rtype: scipy.sparse.csr.csr_matrix
+        :rtype: scipy.sparse.csr_matrix
         :return: Interpolation matrix
 
         .. plot::

SimPEG/Utils/matutils.py

@@ -27,7 +27,7 @@ def mkvc(x, numDims=1):
 
     if isinstance(x, Zero):
         return x
-    
+
     assert isinstance(x, np.ndarray), "Vector must be a numpy array"
 
     if numDims == 1:
@@ -355,9 +355,9 @@ def diagEst(matFun, n, k=None, approach='Probing'):
         2. Ones : random +/- 1 entries
         3. Random : random vectors
 
-    :param lambda (numpy.array) matFun: matrix to estimate the diagonal of
-    :param int64 n: size of the vector that should be used to compute matFun(v)
-    :param int64 k: number of vectors to be used to estimate the diagonal
+    :param callable matFun: takes a (numpy.array) and multiplies it by a matrix to estimate the diagonal
+    :param int n: size of the vector that should be used to compute matFun(v)
+    :param int k: number of vectors to be used to estimate the diagonal
     :param str approach: approach to be used for getting vectors
     :rtype: numpy.array
     :return: est_diag(A)
@@ -422,9 +422,9 @@ class Zero(object):
     def __ge__(self, v):return 0 >= v
     def __gt__(self, v):return 0 > v
 
-    @property 
+    @property
     def transpose(self): return Zero()
-    
+
     @property
     def T(self): return Zero()
 

SimPEG/Utils/meshutils.py

@@ -83,7 +83,7 @@ def closestPoints(mesh, pts, gridLoc='CC'):
     """
         Move a list of points to the closest points on a grid.
 
-        :param simpeg.Mesh.BaseMesh mesh: The mesh
+        :param BaseMesh mesh: The mesh
         :param numpy.ndarray pts: Points to move
         :param string gridLoc: ['CC', 'N', 'Fx', 'Fy', 'Fz', 'Ex', 'Ex', 'Ey', 'Ez']
         :rtype: numpy.ndarray
@@ -104,16 +104,20 @@ def closestPoints(mesh, pts, gridLoc='CC'):
 
 def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
     """
-        Extracts Core Mesh from Global mesh
-        xyzlim: 2D array [ndim x 2]
-        mesh: SimPEG mesh
-        This function ouputs:
-            - actind: corresponding boolean index from global to core
-            - meshcore: core SimPEG mesh
-        Warning: 1D and 2D has not been tested
+    Extracts Core Mesh from Global mesh
+
+    :param numpy.ndarray xyzlim: 2D array [ndim x 2]
+    :param BaseMesh mesh: The mesh
+
+    This function ouputs::
+
+        - actind: corresponding boolean index from global to core
+        - meshcore: core SimPEG mesh
+
+    Warning: 1D and 2D has not been tested
     """
     from SimPEG import Mesh
-    if mesh.dim ==1:
+    if mesh.dim == 1:
         xyzlim = xyzlim.flatten()
         xmin, xmax = xyzlim[0], xyzlim[1]
 
@@ -125,11 +129,11 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
 
         x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
 
-        meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
+        meshCore = Mesh.TensorMesh([hx, hy], x0=x0)
 
         actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax)
 
-    elif mesh.dim ==2:
+    elif mesh.dim == 2:
         xmin, xmax = xyzlim[0,0], xyzlim[0,1]
         ymin, ymax = xyzlim[1,0], xyzlim[1,1]
 
@@ -144,12 +148,12 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
 
         x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
 
-        meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
+        meshCore = Mesh.TensorMesh([hx, hy], x0=x0)
 
         actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
                & (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
 
-    elif mesh.dim==3:
+    elif mesh.dim == 3:
         xmin, xmax = xyzlim[0,0], xyzlim[0,1]
         ymin, ymax = xyzlim[1,0], xyzlim[1,1]
         zmin, zmax = xyzlim[2,0], xyzlim[2,1]
@@ -168,7 +172,7 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
 
         x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5, zc[0]-hz[0]*0.5]
 
-        meshCore = Mesh.TensorMesh([hx, hy, hz] ,x0=x0)
+        meshCore = Mesh.TensorMesh([hx, hy, hz], x0=x0)
 
         actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
                & (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \

SimPEG/__init__.py

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

docs/Makefile

@@ -2,7 +2,7 @@
 #
 
 # You can set these variables from the command line.
-SPHINXOPTS    =
+SPHINXOPTS    = -n -w warnings.txt
 SPHINXBUILD   = sphinx-build
 PAPER         =
 BUILDDIR      = _build

docs/_templates/layout.html

@@ -0,0 +1,22 @@
+{# Import the theme's layout. #}
+{% extends "!layout.html" %}
+
+{% block extrahead %}
+{{ super() }}
+
+<meta name="description" content="Simulation and Parameter Estimation in Geophysics">
+<meta name="author" content="SimPEG Developers">
+<meta name="keywords" content="python, geophysics, inversion, electromagnetics, magnetotellurics, magnetics, gravity, DC, flow inverse problems, open source, finite volume">
+
+
+<script>
+  (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
+  (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
+  m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
+  })(window,document,'script','https://www.google-analytics.com/analytics.js','ga');
+
+  ga('create', 'UA-45185336-1', 'auto');
+  ga('send', 'pageview');
+
+</script>
+{% endblock %}

docs/api_FiniteVolume.rst

@@ -1,19 +0,0 @@
-.. _api_FiniteVolume:
-
-Finite Volume
-*************
-
-Any numerical implementation requires the discretization of continuous functions into discrete approximations. These approximations are typically organized in a mesh, which defines boundaries, locations, and connectivity. Of specific interest to geophysical simulations, we require that averaging, interpolation and differential operators be defined for any mesh. In SimPEG, we have implemented a staggered mimetic finite volume approach (`Hyman and Shashkov, 1999 <http://math.lanl.gov/~mac/papers/numerics/HS99B.pdf>`_). This approach requires the definitions of variables at either cell-centers, nodes, faces, or edges as seen in the figure below.
-
-.. image:: images/finitevolrealestate.png
-   :width: 400 px
-   :alt: FiniteVolume
-   :align: center
-
-
-.. toctree::
-   :maxdepth: 2
-
-   api_Mesh
-   api_DiffOps
-   api_InnerProducts

docs/api_MeshCode.rst

@@ -1,36 +0,0 @@
-.. _api_MeshCode:
-
-Tensor Mesh
-===========
-
-.. automodule:: SimPEG.Mesh.TensorMesh
-    :show-inheritance:
-    :members:
-    :undoc-members:
-
-
-Cylindrical Mesh
-================
-
-.. automodule:: SimPEG.Mesh.CylMesh
-    :show-inheritance:
-    :members:
-    :undoc-members:
-
-
-Tree Mesh
-=========
-
-.. autoclass:: SimPEG.Mesh.TreeMesh.TreeMesh
-    :show-inheritance:
-    :members:
-    :undoc-members:
-
-
-Curvilinear Mesh
-================
-
-.. automodule:: SimPEG.Mesh.CurvilinearMesh
-    :show-inheritance:
-    :members:
-    :undoc-members:

docs/app.yaml

@@ -0,0 +1,95 @@
+# application: simpegdocs
+# version: 1
+runtime: python27
+api_version: 1
+threadsafe: yes
+
+handlers:
+
+# favicon
+- url: /images/logo-block\.ico
+  static_files: /images/logo-block.ico
+  upload: /images/logo-block\.ico
+
+# all css
+- url: /(.*\.css)
+  mime_type: text/css
+  static_files: _build/html/\1
+  upload: _build/html/(.*\.css)
+
+# webfonts
+- url: /(.*\.(eot|svg|ttf|woff|woff2|otf))
+  static_files: _build/html/\1
+  upload: _build/html/(.*\.(eot|svg|ttf|woff|woff2|otf))
+
+# javascript
+- url: /(.*\.js)
+  mime_type: text/javascript
+  static_files: _build/html/\1
+  upload: _build/html/(.*\.js)
+
+# plain text source
+- url: /(.*\.txt)
+  mime_type: text/plain
+  static_files: _build/html/\1
+  upload: _build/html/(.*\.txt)
+
+# images
+- url: /_images/(.*\.(gif|png|jpg|ico))
+  static_files: _build/html/_images/\1
+  upload: _build/html/_images/(.*\.(gif|png|jpg|ico))
+
+# redirect en/latest traffic
+- url: /en/latest/(.*\.html)
+  script: simpegdocs.app
+
+# raw html
+- url: /(.*\.html)
+  mime_type: text/html
+  static_files: _build/html/\1
+  upload: _build/html/(.*\.html)
+
+# serve index files
+- url: /(.+)/
+  static_files: _build/html/\1/index.html
+  upload: _build/html/(.+)/index.html
+
+- url: /(.+)
+  static_files: _build/html/\1/index.html
+  upload: _build/html/(.+)/index.html
+
+- url: /
+  static_files: _build/html/index.html
+  upload: _build/html/index.html
+
+- url: .*
+  script: simpegdocs.app
+
+# Recommended file skipping declaration from the GAE tutorials
+skip_files:
+  - ^(.*/)?app\.yaml
+  - ^(.*/)?app\.yml
+  - ^(.*/)?#.*#
+  - ^(.*/)?.*~
+  - ^(.*/)?.*\.py[co]
+  - ^(.*/)?.*/RCS/.*
+  - ^(.*/)?\..*
+  - ^(.*/)?tests$
+  - ^(.*/)?test$
+  - ^test/(.*/)?
+  - ^COPYING.LESSER
+  - ^README\..*
+  - \.gitignore
+  - ^\.git/.*
+  - \.*\.lint$
+  - ^(.*/)?.*\.doctree$
+
+libraries:
+- name: webapp2
+  version: "2.5.2"
+- name: PIL
+  version: "1.1.7"
+- name: numpy
+  version: "latest"
+- name: jinja2
+  version: "latest"

docs/conf.py

@@ -28,7 +28,7 @@ sys.path.append('../')
 
 # Add any Sphinx extension module names here, as strings. They can be extensions
 # coming with Sphinx (named 'sphinx.ext.*') or your custom ones.
-extensions = ['sphinx.ext.todo', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode', 'sphinx.ext.autodoc', 'matplotlib.sphinxext.plot_directive']
+extensions = ['sphinx.ext.todo', 'sphinx.ext.mathjax', 'sphinx.ext.viewcode', 'sphinx.ext.autodoc', 'sphinx.ext.intersphinx', 'matplotlib.sphinxext.plot_directive']
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']
@@ -44,16 +44,16 @@ master_doc = 'index'
 
 # General information about the project.
 project = u'SimPEG'
-copyright = u'2013, SimPEG Developers'
+copyright = u'2013 - 2016, SimPEG Developers'
 
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
 # built documents.
 #
 # The short X.Y version.
-version = '0.1.10'
+version = '0.1.12'
 # The full version, including alpha/beta/rc tags.
-release = '0.1.10'
+release = '0.1.12'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
@@ -124,12 +124,12 @@ except Exception, e:
 # The name of an image file (within the static path) to use as favicon of the
 # docs.  This file should be a Windows icon file (.ico) being 16x16 or 32x32
 # pixels large.
-#html_favicon = None
+html_favicon = './images/logo-block.ico'
 
 # Add any paths that contain custom static files (such as style sheets) here,
 # relative to this directory. They are copied after the builtin static files,
 # so a file named "default.css" will overwrite the builtin "default.css".
-html_static_path = ['_static']
+html_static_path = []
 
 # If not '', a 'Last updated on:' timestamp is inserted at every page bottom,
 # using the given strftime format.
@@ -229,6 +229,12 @@ man_pages = [
 # If true, show URL addresses after external links.
 #man_show_urls = False
 
+# Intersphinx
+intersphinx_mapping = {'python': ('http://docs.python.org/2', None),
+                       'numpy': ('http://docs.scipy.org/doc/numpy/', None),
+                       'scipy': ('http://docs.scipy.org/doc/scipy/reference/', None),
+                       'matplotlib': ('http://matplotlib.sourceforge.net/', None)}
+
 
 # -- Options for Texinfo output ------------------------------------------------
 
@@ -251,3 +257,36 @@ texinfo_documents = [
 #texinfo_show_urls = 'footnote'
 
 autodoc_member_order = 'bysource'
+
+def supress_nonlocal_image_warn():
+    import sphinx.environment
+    sphinx.environment.BuildEnvironment.warn_node = _supress_nonlocal_image_warn
+
+def _supress_nonlocal_image_warn(self, msg, node):
+    from docutils.utils import get_source_line
+
+    if not msg.startswith('nonlocal image URI found:'):
+        self._warnfunc(msg, '%s:%s' % get_source_line(node))
+
+supress_nonlocal_image_warn()
+
+
+nitpick_ignore = [
+    ('py:class', 'IdentityMap'),
+    ('py:class', 'BaseSurvey'),
+    ('py:class', 'BaseSrc'),
+    ('py:class', 'BaseRx'),
+    ('py:class', 'Survey'),
+    ('py:class', 'FieldsFDEM'),
+    ('py:class', 'Fields3D_e'),
+    ('py:class', 'Fields3D_b'),
+    ('py:class', 'Fields3D_j'),
+    ('py:class', 'Fields3D_h'),
+    ('py:class', 'SurveyTDEM'),
+    ('py:class', 'SrcTDEM'),
+    ('py:class', 'EMPropMap'),
+    ('py:class', 'Data'),
+    ('py:class', 'SurveyDC'),
+    ('py:class', 'BaseMTFields'),
+    ('py:class', 'SolverLU'),
+]

docs/content/api_core/api_Examples.rst

@@ -7,7 +7,7 @@ Examples
     :maxdepth: 1
     :glob:
 
-    examples/*
+    ../examples/*
 
 
 External Notebooks

docs/content/api_core/api_FiniteVolume.rst

@@ -0,0 +1,27 @@
+.. _api_FiniteVolume:
+
+Finite Volume
+*************
+
+Any numerical implementation requires the discretization of continuous
+functions into discrete approximations. These approximations are typically
+organized in a mesh, which defines boundaries, locations, and connectivity. Of
+specific interest to geophysical simulations, we require that averaging,
+interpolation and differential operators be defined for any mesh. In SimPEG,
+we have implemented a staggered mimetic finite volume approach (`Hyman and
+Shashkov, 1999 <http://math.lanl.gov/~mac/papers/numerics/HS99B.pdf>`_). This
+approach requires the definitions of variables at either cell-centers, nodes,
+faces, or edges as seen in the figure below.
+
+.. image:: ../../images/finitevolrealestate.png
+   :width: 400 px
+   :alt: FiniteVolume
+   :align: center
+
+
+.. toctree::
+   :maxdepth: 2
+
+   api_Mesh
+   api_DiffOps
+   api_InnerProducts

docs/content/api_core/api_ForwardProblem.rst

@@ -52,13 +52,15 @@ We can take the derivative of the PDE:
 
     \nabla_m c(m, u) \partial m + \nabla_u c(m, u) \partial u = 0
 
-If the forward problem is invertible, then we can rearrange for \\(\\frac{\\partial u}{\\partial m}\\):
+If the forward problem is invertible, then we can rearrange for
+\\(\\frac{\\partial u}{\\partial m}\\):
 
 .. math::
 
     J = - P \left( \nabla_u c(m, u) \right)^{-1} \nabla_m c(m, u)
 
-This can often be computed given a vector (i.e. \\(J(v)\\)) rather than stored, as \\(J\\) is a large dense matrix.
+This can often be computed given a vector (i.e. \\(J(v)\\)) rather than
+stored, as \\(J\\) is a large dense matrix.
 
 
 
@@ -67,13 +69,45 @@ The API
 
 Problem
 -------
-.. automodule:: SimPEG.Problem
+
+.. autoclass:: SimPEG.Problem.BaseProblem
+    :members:
+    :undoc-members:
+
+.. autoclass:: SimPEG.Problem.BaseTimeProblem
+    :members:
+    :undoc-members:
+
+Fields
+------
+
+.. autoclass:: SimPEG.Fields.Fields
+    :members:
+    :undoc-members:
+
+.. autoclass:: SimPEG.Fields.TimeFields
     :members:
     :undoc-members:
 
 Survey
 ------
-.. automodule:: SimPEG.Survey
+
+.. autoclass:: SimPEG.Survey.BaseSurvey
     :members:
     :undoc-members:
 
+.. autoclass:: SimPEG.Survey.BaseSrc
+    :members:
+    :undoc-members:
+
+.. autoclass:: SimPEG.Survey.BaseRx
+    :members:
+    :undoc-members:
+
+.. autoclass:: SimPEG.Survey.BaseTimeRx
+    :members:
+    :undoc-members:
+
+.. autoclass:: SimPEG.Survey.Data
+    :members:
+    :undoc-members:

docs/content/api_core/api_InnerProducts.rst

@@ -4,7 +4,10 @@
 Inner Products
 **************
 
-By using the weak formulation of many of the PDEs in geophysical applications, we can rapidly develop discretizations. Much of this work, however, needs a good understanding of how to approximate inner products on our discretized meshes. We will define the inner product as:
+By using the weak formulation of many of the PDEs in geophysical applications,
+we can rapidly develop discretizations. Much of this work, however, needs a
+good understanding of how to approximate inner products on our discretized
+meshes. We will define the inner product as:
 
 .. math::
 
@@ -14,12 +17,15 @@ where a and b are either scalars or vectors.
 
 .. note::
 
-    The InnerProducts class is a base class providing inner product matrices for meshes and cannot run on its own.
+    The InnerProducts class is a base class providing inner product matrices
+    for meshes and cannot run on its own.
 
 
 Example problem for DC resistivity
 ----------------------------------
-We will start with the formulation of the Direct Current (DC) resistivity problem in geophysics.
+
+We will start with the formulation of the Direct Current (DC) resistivity
+problem in geophysics.
 
 
 .. math::
@@ -28,12 +34,13 @@ We will start with the formulation of the Direct Current (DC) resistivity proble
 
         \nabla\cdot \vec{j} = q
 
-In the following discretization, \\\( \\sigma \\\) and \\\( \\phi \\\)
-will be discretized on the cell-centers and the flux, \\\(\\vec{j}\\\),
+In the following discretization, :math:`\sigma` and :math:`\phi`
+will be discretized on the cell-centers and the flux, :math:`\vec{j}`,
 will be on the faces. We will use the weak formulation to discretize
 the DC resistivity equation.
 
-We can define in weak form by integrating with a general face function \\\(\\vec{f}\\\):
+We can define in weak form by integrating with a general face function
+:math:`\vec{f}`:
 
 .. math::
 
@@ -61,9 +68,16 @@ We can then discretize for every cell:
 
 .. note::
 
-    We have discretized the dot product above, but remember that we do not really have a single vector \\\(\\mathbf{J}\\\), but approximations of \\\(\\vec{j}\\\) on each face of our cell. In 2D that means 2 approximations of \\\(\\mathbf{J}_x\\\) and 2 approximations of \\\(\\mathbf{J}_y\\\). In 3D we also have 2 approximations of \\\(\\mathbf{J}_z\\\).
+    We have discretized the dot product above, but remember that we do not
+    really have a single vector :math:`\mathbf{J}`, but approximations of
+    :math:`\vec{j}` on each face of our cell. In 2D that means 2
+    approximations of :math:`\mathbf{J}_x` and 2 approximations of
+    :math:`\mathbf{J}_y`. In 3D we also have 2 approximations of
+    :math:`\mathbf{J}_z`.
 
-Regardless of how we choose to approximate this dot product, we can represent this in vector form (again this is for every cell), and will generalize for the case of anisotropic (tensor) sigma.
+Regardless of how we choose to approximate this dot product, we can represent
+this in vector form (again this is for every cell), and will generalize for
+the case of anisotropic (tensor) sigma.
 
 .. math::
 
@@ -71,14 +85,17 @@ Regardless of how we choose to approximate this dot product, we can represent th
     -\phi^{\top} v_{\text{cell}} \mathbf{D}_{\text{cell}} \mathbf{F})
     + \text{BC}
 
-We multiply by  square-root of volume on each side of the tensor conductivity to keep symmetry in the system. Here \\\(\\mathbf{J}_c\\\) is the Cartesian \\\(\\mathbf{J}\\\) (on the faces that we choose to use in our approximation) and must be calculated differently depending on the mesh:
+We multiply by  square-root of volume on each side of the tensor conductivity
+to keep symmetry in the system. Here :math:`\mathbf{J}_c` is the Cartesian
+:math:`\mathbf{J}` (on the faces that we choose to use in our approximation)
+and must be calculated differently depending on the mesh:
 
 .. math::
     \mathbf{J}_c = \mathbf{Q}_{(i)}\mathbf{J}_\text{TENSOR} \\
     \mathbf{J}_c = \mathbf{N}_{(i)}^{-1}\mathbf{Q}_{(i)}\mathbf{J}_\text{Curv}
 
-Here the \\\(i\\\) index refers to where we choose to approximate this integral, as discussed in the note above.
-We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix \\\( \\mathbf{Q}_{(i)} \\\) to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like:
+Here the :math:`i` index refers to where we choose to approximate this integral, as discussed in the note above.
+We will approximate this integral by taking the fluxes clustered around every node of the cell, there are 8 combinations in 3D, and 4 in 2D. We will use a projection matrix :math:`\mathbf{Q}_{(i)}` to pick the appropriate fluxes. So, now that we have 8 approximations of this integral, we will just take the average. For the TensorMesh, this looks like:
 
 .. math::
 
@@ -107,10 +124,12 @@ By defining the faceInnerProduct (8 combinations of fluxes in 3D, 4 in 2D, 2 in
         \sum_{i=1}^{2^d}
         \mathbf{P}_{(i)}^{\top} \Sigma^{-1} \mathbf{P}_{(i)}
 
-Where \\\(d\\\) is the dimension of the mesh.
-The \\\( \\mathbf{M}^f \\\) is returned when given the input of \\\( \\Sigma^{-1} \\\).
+Where :math:`d` is the dimension of the mesh.
+The :math:`\mathbf{M}^f` is returned when given the input of :math:`\Sigma^{-1}`.
 
-Here each \\( \\mathbf{P} \\in \\mathbb{R}^{(d*nC, nF)} \\\) is a combination of the projection, volume, and any normalization to Cartesian coordinates (where the dot product is well defined):
+Here each :math:`\mathbf{P} ~ \in ~ \mathbb{R}^{(d*nC, nF)}` is a combination
+of the projection, volume, and any normalization to Cartesian coordinates
+(where the dot product is well defined):
 
 .. math::
 
@@ -129,7 +148,10 @@ If ``returnP=True`` is requested in any of these methods the projection matrices
     # In 1D
     P = [P0, P1]
 
-The derivation for ``edgeInnerProducts`` is exactly the same, however, when we approximate the integral using the fields around each node, the projection matrices look a bit different because we have 12 edges in 3D instead of just 6 faces. The interface to the code is exactly the same.
+The derivation for ``edgeInnerProducts`` is exactly the same, however, when we
+approximate the integral using the fields around each node, the projection
+matrices look a bit different because we have 12 edges in 3D instead of just 6
+faces. The interface to the code is exactly the same.
 
 
 Defining Tensor Properties
@@ -137,7 +159,8 @@ Defining Tensor Properties
 
 **For 3D:**
 
-Depending on the number of columns (either 1, 3, or 6) of mu, the material property is interpreted as follows:
+Depending on the number of columns (either 1, 3, or 6) of mu, the material
+property is interpreted as follows:
 
 .. math::
 
@@ -188,13 +211,16 @@ Which is nice and easy to invert if necessary, however, in the fully anisotropic
 Taking Derivatives
 ------------------
 
-We will take the derivative of the fully anisotropic tensor for a 3D mesh, the other cases are easier and will not be discussed here. Let us start with one part of the sum which makes up \\\(\\mathbf{M}^f_\\Sigma\\\) and take the derivative when this is multiplied by some vector \\\(\\mathbf{v}\\\):
+We will take the derivative of the fully anisotropic tensor for a 3D mesh, the
+other cases are easier and will not be discussed here. Let us start with one
+part of the sum which makes up :math:`\mathbf{M}^f_\Sigma` and take the
+derivative when this is multiplied by some vector :math:`\mathbf{v}`:
 
 .. math::
 
     \mathbf{P}^\top \boldsymbol{\Sigma} \mathbf{Pv}
 
-Here we will let \\\( \\mathbf{Pv} = \\mathbf{y} \\\) and \\\(\\mathbf{y}\\\) will have the form:
+Here we will let :math:`\mathbf{Pv} = \mathbf{y}` and :math:`\mathbf{y}` will have the form:
 
 .. math::
 
@@ -233,7 +259,9 @@ Here we will let \\\( \\mathbf{Pv} = \\mathbf{y} \\\) and \\\(\\mathbf{y}\\\) wi
         \end{matrix}
     \right]
 
-Now it is easy to take the derivative with respect to any one of the parameters, for example, \\\(\\frac{\\partial}{\\partial\\boldsymbol{\\sigma}_1}\\\)
+Now it is easy to take the derivative with respect to any one of the
+parameters, for example,
+:math:`\frac{\partial}{\partial\boldsymbol{\sigma}_1}`
 
 .. math::
     \frac{\partial}{\partial \boldsymbol{\sigma}_1}\left(\mathbf{P}^\top\Sigma\mathbf{y}\right)
@@ -247,7 +275,8 @@ Now it is easy to take the derivative with respect to any one of the parameters,
         \end{matrix}
     \right]
 
-Whereas \\\(\\frac{\\partial}{\\partial\\boldsymbol{\\sigma}_4}\\\), for example, is:
+Whereas :math:`\frac{\partial}{\partial\boldsymbol{\sigma}_4}`, for
+example, is:
 
 .. math::
     \frac{\partial}{\partial \boldsymbol{\sigma}_4}\left(\mathbf{P}^\top\Sigma\mathbf{y}\right)
@@ -261,11 +290,12 @@ Whereas \\\(\\frac{\\partial}{\\partial\\boldsymbol{\\sigma}_4}\\\), for example
         \end{matrix}
     \right]
 
-These are computed for each of the 8 projections, horizontally concatenated, and returned.
+These are computed for each of the 8 projections, horizontally concatenated,
+and returned.
 
 The API
 -------
 
-.. automodule:: SimPEG.Mesh.InnerProducts
+.. autoclass:: SimPEG.Mesh.InnerProducts.InnerProducts
     :members:
     :undoc-members:

docs/content/api_core/api_Inversion.rst

@@ -3,7 +3,7 @@
 InvProblem
 **********
 
-.. automodule:: SimPEG.InvProblem
+.. autoclass:: SimPEG.InvProblem.BaseInvProblem
     :show-inheritance:
     :members:
     :undoc-members:
@@ -12,7 +12,7 @@ InvProblem
 Inversion
 *********
 
-.. automodule:: SimPEG.Inversion
+.. autoclass:: SimPEG.Inversion.BaseInversion
     :show-inheritance:
     :members:
     :undoc-members:

docs/content/api_core/api_Maps.rst

@@ -27,7 +27,8 @@ back to conductivity. This is a relatively trivial example (we are just taking
 the exponential!) but by defining maps we can start to combine and manipulate
 exactly what we think about as our model, \\\(m\\\). In code, this looks like
 
-::
+.. code-block:: python
+    :linenos:
 
     M = Mesh.TensorMesh([100]) # Create a mesh
     expMap = Maps.ExpMap(M)    # Create a mapping
@@ -46,14 +47,15 @@ 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.Vertical1DMap`),
+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.Vertical1DMap(M)
+    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!
@@ -61,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.Vertical1DMap(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).
@@ -122,6 +106,8 @@ When these are used in the inverse problem, this is extremely important!!
 The API
 =======
 
+The :code:`IdentityMap` is the base class for all mappings, and it does absolutely nothing.
+
 .. autoclass:: SimPEG.Maps.IdentityMap
     :members:
     :undoc-members:
@@ -130,7 +116,6 @@ The API
 Common Maps
 ===========
 
-
 Exponential Map
 ---------------
 
@@ -148,7 +133,7 @@ lives (i.e. it varies logarithmically).
 Vertical 1D Map
 ---------------
 
-.. autoclass:: SimPEG.Maps.Vertical1DMap
+.. autoclass:: SimPEG.Maps.SurjectVertical1D
     :members:
     :undoc-members:
 
@@ -164,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
@@ -196,8 +160,8 @@ Mesh to Mesh Map
     :undoc-members:
 
 
-Some Extras
-===========
+Under the Hood
+==============
 
 Combo Map
 ---------

docs/content/api_core/api_Mesh.rst

@@ -188,6 +188,6 @@ other types of meshes in this SimPEG framework.
 The API
 =======
 
-.. automodule:: SimPEG.Mesh.BaseMesh
+.. autoclass:: SimPEG.Mesh.BaseMesh.BaseMesh
     :members:
     :undoc-members:

docs/content/api_core/api_MeshCode.rst

@@ -0,0 +1,68 @@
+.. _api_MeshCode:
+
+Tensor Mesh
+===========
+
+.. autoclass:: SimPEG.Mesh.TensorMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+Cylindrical Mesh
+================
+
+.. autoclass:: SimPEG.Mesh.CylMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+Tree Mesh
+=========
+
+.. autoclass:: SimPEG.Mesh.TreeMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+Curvilinear Mesh
+================
+
+.. autoclass:: SimPEG.Mesh.CurvilinearMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Base Rectangular Mesh
+=====================
+
+.. autoclass:: SimPEG.Mesh.BaseMesh.BaseRectangularMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+Base Tensor Mesh
+================
+
+.. autoclass:: SimPEG.Mesh.TensorMesh.BaseTensorMesh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Mesh IO
+=======
+
+.. automodule:: SimPEG.Mesh.MeshIO
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+
+Mesh Viewing
+============
+
+.. automodule:: SimPEG.Mesh.View
+    :members:
+    :undoc-members:
+    :show-inheritance:

docs/content/api_core/api_PropMaps.rst

@@ -0,0 +1,29 @@
+SimPEG PropMaps
+***************
+
+The API
+=======
+
+Property
+--------
+
+.. autoclass:: SimPEG.PropMaps.Property
+    :members:
+    :undoc-members:
+
+
+PropMap
+-------
+
+.. autoclass:: SimPEG.PropMaps.PropMap
+    :members:
+    :undoc-members:
+
+
+PropModel
+---------
+
+.. autoclass:: SimPEG.PropMaps.PropModel
+    :members:
+    :undoc-members:
+

docs/content/api_core/api_Regularization.rst

@@ -91,10 +91,21 @@ The API
     :members:
     :undoc-members:
 
+.. autoclass:: SimPEG.Regularization.Simple
+    :show-inheritance:
+    :members:
 
 .. autoclass:: SimPEG.Regularization.Tikhonov
     :show-inheritance:
     :members:
 
+.. autoclass:: SimPEG.Regularization.Sparse
+    :show-inheritance:
+    :members:
+
+.. autoclass:: SimPEG.Regularization.RegularizationMesh
+    :show-inheritance:
+    :members:
+
 
 

docs/content/api_core/api_Solver.rst

@@ -46,6 +46,8 @@ The API
 =======
 
 .. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapD
+    :noindex:
 
 .. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapI
+    :noindex:
 

docs/content/api_core/api_Utilities.rst

@@ -6,5 +6,6 @@ Utilities
 
    api_Solver
    api_Maps
+   api_PropMaps
    api_Utils
    api_Tests

docs/content/api_core/api_Utils.rst

@@ -21,7 +21,7 @@ Solver Utilities
     :undoc-members:
 
 Curv Utilities
-=============
+==============
 
 .. automodule:: SimPEG.Utils.curvutils
     :members:
@@ -51,7 +51,9 @@ Interpolation Utilities
 Counter Utilities
 =================
 
-::
+.. code-block:: python
+    :linenos:
+
     class MyClass(object):
         def __init__(self, url):
             self.counter = Counter()
@@ -69,7 +71,9 @@ Counter Utilities
     for i in range(300): c.MySecondMethod()
     c.counter.summary()
 
-::
+
+.. code-block:: text
+    :linenos:
 
     Counters:
       MyClass.MyMethod                        :      100
@@ -77,6 +81,8 @@ Counter Utilities
     Times:                                        mean      sum
       MyClass.MySecondMethod                  : 1.70e-06, 5.10e-04,  300x
 
+
+
 The API
 -------
 

docs/content/api_core/api_bigPicture.rst

@@ -35,7 +35,7 @@ The Big Picture
 Defining a well-posed inverse problem and solving it is a complex task that requires many components that must interact. It is helpful
 to view this task as a workflow in which various elements are explicitly identified and integrated. The figure below outlines the inversion components that consists of inputs, implementation, and evaluation. The inputs are composed of the geophysical data, the equations which are a mathematical description of the governing physics, and prior knowledge or assumptions about the setting. The implementation consists of two broad categories: the forward simulation and the inversion. The **forward simulation** is the means by which we solve the governing equations given a model and the **inversion components** evaluate and update this model. We are considering a gradient based approach, which updates the model through an optimization routine. The output of this implementation is a model, which, prior to interpretation, must be evaluated. This requires considering, and often re-assessing, the choices and assumptions made in both the input and implementation stages.
 
-.. image:: InversionWorkflow-PreSimPEG.png
+.. image:: ../../images/InversionWorkflow-PreSimPEG.png
    :width: 400 px
    :alt: Components
    :align: center
@@ -46,24 +46,24 @@ A Comprehensive Framework
 
 There are an overwhelming amount of choices to be made as one works through the forward modeling and inversion process (see figure above). As a result, software implementations of this workflow often become complex and highly interdependent, making it difficult to interact with and to ask other scientists to pick up and change. Our approach to handling this complexity is to propose a framework, (see below), that compartmentalizes the implementation of inversions into various units. We present it in this specific modular style, as each unit contains a targeted subset of choices crucial to the inversion process.
 
-.. image:: InversionWorkflow.png
+.. image:: ../../images/InversionWorkflow.png
    :width: 400 px
    :alt: Framework
    :align: center
 
-The process of obtaining an acceptable model from an inversion generally requires the geophysicist to perform several iterations of the inversion workflow, rethinking and redesigning each piece of the framework to ensure it is appropriate in the current context. Inversions are experimental and empirical by nature and our software package is designed to facilitate this iterative process. To accomplish this, we have divided the inversion methodology into eight major components (See figure above). The (:class:`SimPEG.Mesh.BaseMesh`) class handles the discretization of the earth and also provides numerical operators. The forward simulation is split into two classes, the (:class:`SimPEG.Survey.BaseSurvey`) and the (:class:`SimPEG.Problem.BaseProblem`). The (:class:`SimPEG.Survey.BaseSurvey`) class handles the geometry of a geophysical problem as well as sources. The (:class:`SimPEG.Problem.BaseProblem`) class handles the simulation of the physics for the geophysical problem of interest. Although created independently, these two classes must be paired to form all of the components necessary for a geophysical forward simulation and calculation of the sensitivity. The (:class:`SimPEG.Problem.BaseProblem`) creates geophysical fields given a source from the (:class:`SimPEG.Survey.BaseSurvey`). The (:class:`SimPEG.Survey.BaseSurvey`) interpolates these fields to the receiver locations and converts them to the appropriate data type, for example, by selecting only the measured components of the field. Each of these operations may have associated derivatives with respect to the model and the computed field; these are included in the calculation of the sensitivity. For the inversion, a (:class:`SimPEG.DataMisfit.BaseDataMisfit`) is chosen to capture the goodness of fit of the predicted data and a (:class:`SimPEG.Regularization.BaseRegularization`) is chosen to handle the non-uniqueness. These inversion elements and an Optimization routine are combined into an inverse problem class (:class:`SimPEG.InvProblem.BaseInvProblem`). (:class:`SimPEG.InvProblem.BaseInvProblem`) is the mathematical statement that will be numerically solved by running an Inversion. The (:class:`SimPEG.Inversion.BaseInversion`) class handles organization and dispatch of directives between all of the various pieces of the framework.
+The process of obtaining an acceptable model from an inversion generally requires the geophysicist to perform several iterations of the inversion workflow, rethinking and redesigning each piece of the framework to ensure it is appropriate in the current context. Inversions are experimental and empirical by nature and our software package is designed to facilitate this iterative process. To accomplish this, we have divided the inversion methodology into eight major components (See figure above). The :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class handles the discretization of the earth and also provides numerical operators. The forward simulation is split into two classes, the :class:`SimPEG.Survey.BaseSurvey` and the :class:`SimPEG.Problem.BaseProblem`. The :class:`SimPEG.Survey.BaseSurvey` class handles the geometry of a geophysical problem as well as sources. The :class:`SimPEG.Problem.BaseProblem` class handles the simulation of the physics for the geophysical problem of interest. Although created independently, these two classes must be paired to form all of the components necessary for a geophysical forward simulation and calculation of the sensitivity. The :class:`SimPEG.Problem.BaseProblem` creates geophysical fields given a source from the :class:`SimPEG.Survey.BaseSurvey`. The :class:`SimPEG.Survey.BaseSurvey` interpolates these fields to the receiver locations and converts them to the appropriate data type, for example, by selecting only the measured components of the field. Each of these operations may have associated derivatives with respect to the model and the computed field; these are included in the calculation of the sensitivity. For the inversion, a :class:`SimPEG.DataMisfit.BaseDataMisfit` is chosen to capture the goodness of fit of the predicted data and a :class:`SimPEG.Regularization.BaseRegularization` is chosen to handle the non-uniqueness. These inversion elements and an Optimization routine are combined into an inverse problem class :class:`SimPEG.InvProblem.BaseInvProblem`. :class:`SimPEG.InvProblem.BaseInvProblem` is the mathematical statement that will be numerically solved by running an Inversion. The :class:`SimPEG.Inversion.BaseInversion` class handles organization and dispatch of directives between all of the various pieces of the framework.
 
-The arrows in the figure above indicate what each class takes as a primary argument. For example, both the (:class:`SimPEG.Problem.BaseProblem`) and (:class:`SimPEG.Regularization.BaseRegularization`) classes take a (:class:`SimPEG.Mesh.BaseMesh`) class as an argument. The diagram does not show class inheritance, as each of the base classes outlined have many subtypes that can be interchanged. The (:class:`SimPEG.Mesh.BaseMesh`) class, for example, could be a regular Cartesian mesh (:class:`SimPEG.Mesh.TensorMesh`) or a cylindrical coordinate mesh (:class:`SimPEG.Mesh.CylMesh`), which have many properties in common. These common features, such as both meshes being created from tensor products, can be exploited through inheritance of base classes, and differences can be expressed through subtype polymorphism. Please look at the documentation here for more in-depth information.
+The arrows in the figure above indicate what each class takes as a primary argument. For example, both the :class:`SimPEG.Problem.BaseProblem` and :class:`SimPEG.Regularization.BaseRegularization` classes take a :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class as an argument. The diagram does not show class inheritance, as each of the base classes outlined have many subtypes that can be interchanged. The :class:`SimPEG.Mesh.BaseMesh.BaseMesh` class, for example, could be a regular Cartesian mesh :class:`SimPEG.Mesh.TensorMesh` or a cylindrical coordinate mesh :class:`SimPEG.Mesh.CylMesh`, which have many properties in common. These common features, such as both meshes being created from tensor products, can be exploited through inheritance of base classes, and differences can be expressed through subtype polymorphism. Please look at the documentation here for more in-depth information.
 
 
-.. include:: ../CITATION.rst
+.. include:: ../../../CITATION.rst
 
 Authors
 -------
 
-.. include:: ../AUTHORS.rst
+.. include:: ../../../AUTHORS.rst
 
 License
 -------
 
-.. include:: ../LICENSE
+.. include:: ../../../LICENSE

docs/content/dc/index.rst

@@ -1,5 +1,3 @@
-.. _api_DC:
-
 .. math::
 
     \renewcommand{\div}{\nabla\cdot\,}
@@ -38,8 +36,16 @@
     \renewcommand {\u}  { {\vec u} }
     \newcommand{\I}{\vec{I}}
 
+
+Direct Current Resistivity
+**************************
+
+`SimPEG.DCIP` uses SimPEG as the framework for the forward and inverse
+direct current (DC) resistivity and induced polarization (IP) geophysical problems.
+
+
 DC resistivity survey
-*********************
+=====================
 
 Electrical resistivity of subsurface materials is measured by causing an electrical current to flow in the earth between one pair of electrodes while the voltage across a second pair of electrodes is measured. The result is an "apparent" resistivity which is a value representing the weighted average resistivity over a volume of the earth. Variations in this measurement are caused by variations in the soil, rock, and pore fluid electrical resistivity. Surveys require contact with the ground, so they can be labour intensive. Results are sometimes interpreted directly, but more commonly, 1D, 2D or 3D models are estimated using inversion procedures (`GPG <http://www.eos.ubc.ca/courses/eosc350/content/>`_).
 
@@ -55,7 +61,7 @@ As direct current (DC) implies, in DC resistivity survey, we assume steady-state
 
     \curl \e = 0
 
-Then by taking \\(\\curl\\) for the first equation, we have
+Then by taking \\(\\div\\) of the first equation, we have
 
 .. math::
 
@@ -137,13 +143,14 @@ Comparing to the analytic function:
 
 .. plot::
 
-    import simpegDC as DC
-    DC.Examples.Verification.run(plotIt=True)
+    from SimPEG import Examples
+    Examples.DC_Analytic_Dipole.run(plotIt=True)
 
-API
-===
 
-.. automodule:: simpegDC.BaseDC
+API for DC codes
+================
+
+.. automodule:: SimPEG.DCIP.BaseDC
     :show-inheritance:
     :members:
     :undoc-members:

docs/content/em/api_FDEM.rst

@@ -9,17 +9,28 @@
 Frequency Domain Electromagnetics
 *********************************
 
-Electromagnetic (EM) geophysical methods are used in a variety of applications from resource exploration, including for hydrocarbons and minerals, to environmental applications, such as groundwater monitoring. The primary physical property of interest in EM is electrical conductivity, which describes the ease with which electric current flows through a material.
+Electromagnetic (EM) geophysical methods are used in a variety of applications
+from resource exploration, including for hydrocarbons and minerals, to
+environmental applications, such as groundwater monitoring. The primary
+physical property of interest in EM is electrical conductivity, which
+describes the ease with which electric current flows through a material.
 
 
 Background
 ==========
 
-Electromagnetic phenomena are governed by Maxwell's equations. They describe the behavior of EM fields and fluxes. Electromagnetic theory for geophysical applications by Ward and Hohmann (1988) is a highly recommended resource on this topic.
+Electromagnetic phenomena are governed by Maxwell's equations. They describe
+the behavior of EM fields and fluxes. Electromagnetic theory for geophysical
+applications by Ward and Hohmann (1988) is a highly recommended resource on
+this topic.
 
 Fourier Transform Convention
 ----------------------------
-In order to examine Maxwell's equations in the frequency domain, we must first define our choice of harmonic time-dependence by choosing a Fourier transform convention. We use the :math:`e^{i \omega t}` convention, so we define our Fourier Transform pair as
+
+In order to examine Maxwell's equations in the frequency domain, we must first
+define our choice of harmonic time-dependence by choosing a Fourier transform
+convention. We use the :math:`e^{i \omega t}` convention, so we define our
+Fourier Transform pair as
 
 .. math ::
     F(\omega) = \int_{-\infty}^{\infty} f(t) e^{- i \omega t} dt \\
@@ -31,6 +42,7 @@ where :math:`\omega` is angular frequency, :math:`t` is time, :math:`F(\omega)`
 
 Maxwell's Equations
 ===================
+
 In the frequency domain, Maxwell's equations are given by
 
 .. math ::
@@ -104,19 +116,20 @@ The H-J formulation is in terms of the current density and the magnetic field:
 
 Discretizing
 ------------
+
 For both formulations, we use a finite volume discretization
 and discretize fields on cell edges, fluxes on cell faces and
 physical properties in cell centers. This is particularly
 important when using symmetry to reduce the dimensionality of a problem
 (for instance on a 2D CylMesh, there are :math:`r`, :math:`z` faces and :math:`\theta` edges)
 
-.. figure:: ../images/finitevolrealestate.png
+.. figure:: ../../images/finitevolrealestate.png
     :align: center
     :scale: 60 %
 
 For the two formulations, the discretization of the physical properties, fields and fluxes are summarized below.
 
-.. figure:: ../images/ebjhdiscretizations.png
+.. figure:: ../../images/ebjhdiscretizations.png
     :align: center
     :scale: 60 %
 
@@ -150,7 +163,7 @@ API
 FDEM Problem
 ------------
 
-.. automodule:: SimPEG.EM.FDEM.FDEM
+.. automodule:: SimPEG.EM.FDEM.ProblemFDEM
     :show-inheritance:
     :members:
     :undoc-members:
@@ -169,6 +182,11 @@ FDEM Survey
     :members:
     :undoc-members:
 
+.. automodule:: SimPEG.EM.FDEM.RxFDEM
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
 FDEM Fields
 -----------
 

docs/content/em/api_TDEM.rst

@@ -359,7 +359,7 @@ TDEM - B formulation
 Field Storage
 =============
 
-.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.FieldsTDEM
+.. autoclass:: SimPEG.EM.TDEM.BaseTDEM.FieldsTDEM
     :show-inheritance:
     :members:
     :undoc-members:

docs/content/em/api_basic.rst

@@ -0,0 +1,33 @@
+Overview of Electromagnetics in SimPEG
+**************************************
+
+
+The API
+=======
+
+Physical Properties
+-------------------
+
+.. autoclass:: SimPEG.EM.Base.EMPropMap
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
+Problem
+-------
+
+.. autoclass:: SimPEG.EM.Base.BaseEMProblem
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
+
+Survey
+------
+
+.. autoclass:: SimPEG.EM.Base.BaseEMSurvey
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
+

docs/content/em/index.rst

@@ -3,22 +3,23 @@ Electromagnetics
 ================
 
 `SimPEG.EM` uses SimPEG as the framework for the forward and inverse
-electromagnetics geophysical problems. 
+electromagnetics geophysical problems.
 
 To solve for predicted data, we follow the framework shown below. The model is
 what we invert for. This is mapped to a physical property on the simulation
 mesh. A source which is used to excite the system is specified. Having a model
 and a source, we can solve Maxwell's equations for fields. We sample these
-fields with recievers to give us predicted data. 
+fields with recievers to give us predicted data.
 
 
-.. image:: ../images/simpegEM_noMath.png
+.. image:: ../../images/simpegEM_noMath.png
     :scale: 50%
 
 
 .. toctree::
    :maxdepth: 2
 
+   api_basic
    api_FDEM
    api_TDEM
    api_Utils