cleanup in mapping tests

Parametric maps renamed to ParametricXXXMap. Some Pep8-ing
light pep8 formatting in tests
2026-09-13 13:03:14 +08:00 · 2016-07-20 09:24:32 -07:00 · 2016-07-18 10:31:25 -07:00 · 2016-07-16 16:26:36 -07:00 · 2016-07-16 16:26:14 -07:00 · 2016-07-16 05:42:06 -07:00
127 changed files with 2963 additions and 1489 deletions
@@ -1,4 +1,4 @@
 [bumpversion]
-current_version = 0.1.10
+current_version = 0.1.12
 files = setup.py SimPEG/__init__.py docs/conf.py

@@ -39,3 +39,5 @@ nosetests.xml
 *.sublime-workspace
 docs/_build/
 Makefile
+docs/warnings.txt
+.DS_Store
@@ -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 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,12 +53,28 @@ 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

+  - 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:
    - rowanc1@gmail.com
@@ -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

 ======
@@ -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::
@@ -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::
@@ -1,12 +1,16 @@
-from SimPEG import np
+from SimPEG import np, Utils
 import BaseDC as DC
 import BaseDC as IP
+import warnings

 def getActiveindfromTopo(mesh, topo):
 # def genActiveindfromTopo(mesh, topo):
    """
        Get active indices from topography
    """
+    warnings.warn(
+            "`getActiveindfromTopo` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead",
+            FutureWarning)
    from scipy.interpolate import NearestNDInterpolator
    if mesh.dim==3:
        nCxy = mesh.nCx*mesh.nCy
@@ -28,6 +32,9 @@ def gettopoCC(mesh, airind):
    """
        Get topography from active indices of mesh.
    """
+    warnings.warn(
+            "`gettopoCC` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead",
+            FutureWarning)
    mesh2D = Mesh.TensorMesh([mesh.hx, mesh.hy], mesh.x0[:2])
    zc = mesh.gridCC[:,2]
    AIRIND = airind.reshape((mesh.vnC[0]*mesh.vnC[1],mesh.vnC[2]), order='F')
@@ -118,34 +125,27 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"):

 def readUBC_DC2DModel(fileName):
    """
-        Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg
+    Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg

-        Input:
-        :param fileName, path to the UBC GIF 2D model file
-
-        Output:
-        :param SimPEG TensorMesh 2D object
-        :return
-
-        Created on Thu Nov 12 13:14:10 2015
-
-        @author: dominiquef
+    :param string fileName: path to the UBC GIF 2D model file
+    :rtype: TensorMesh
+    :return: SimPEG TensorMesh 2D object

    """
    from SimPEG import np, mkvc

    # Open fileand skip header... assume that we know the mesh already
-    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
+    obsfile = np.genfromtxt(fileName, delimiter=' \n', dtype=np.str, comments='!')

-    dim = np.array(obsfile[0].split(),dtype=float)
+    dim = np.array(obsfile[0].split(), dtype=float)

-    temp = np.array(obsfile[1].split(),dtype=float)
+    temp = np.array(obsfile[1].split(), dtype=float)

    if len(temp) > 1:
        model = np.zeros(dim)

        for ii in range(len(obsfile)-1):
-            mm = np.array(obsfile[ii+1].split(),dtype=float)
+            mm = np.array(obsfile[ii+1].split(), dtype=float)
            model[:,ii] = mm

        model = model[:,::-1]
@@ -153,10 +153,10 @@ def readUBC_DC2DModel(fileName):
    else:

        if len(obsfile[1:])==1:
-            mm = np.array(obsfile[1:].split(),dtype=float)
+            mm = np.array(obsfile[1:].split(), dtype=float)

        else:
-            mm = np.array(obsfile[1:],dtype=float)
+            mm = np.array(obsfile[1:], dtype=float)

        # Permute the second dimension to flip the order
        model = mm.reshape(dim[1],dim[0])
@@ -169,23 +169,19 @@ def readUBC_DC2DModel(fileName):

    return model

-def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None):
+
+def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt', clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None):
    """
-        Read list of 2D tx-rx location and plot a speudo-section of apparent
-        resistivity.
+    Read list of 2D tx-rx location and plot a speudo-section of apparent
+    resistivity.

-        Assumes flat topo for now...
+    Assumes flat topo for now...

-        Input:
-        :param d2D, z0
-        :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole)
-        :switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential)
-        Output:
-        :figure scatter plot overlayed on image
-
-        Edited Feb 17th, 2016
-
-        @author: dominiquef
+    :param SurveyDC DCsurvey:
+    :param string surveyType: Either 'pole-dipole'  | 'dipole-dipole'
+    :param string unitType: Either 'appResistivity' | 'appConductivity'  | 'volt'
+    :rtype: matplotlib.plt
+    :return: figure scatter plot overlayed on image

    """
    from SimPEG import np
@@ -218,39 +214,39 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl
        Cmid = (Tx[0][0] + Tx[1][0])/2
        Pmid = (Rx[0][:,0] + Rx[1][:,0])/2

-        # Change output for dtype
-        if dtype == 'volt':
+        # Change output for unitType
+        if unitType == 'volt':

            rho = np.hstack([rho,data])

        else:

            # Compute pant leg of apparent rho
-            if stype == 'pdp':
+            if surveyType == 'pole-dipole':

                leg =  data * 2*np.pi  * MA * ( MA + MN ) / MN

-            elif stype == 'dpdp':
+            elif surveyType == 'dipole-dipole':

                leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )

            else:
-                print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
+                print """unitType must be 'pole-dipole' | 'dipole-dipole' """
                break


-            if dtype == 'appc':
+            if unitType == 'appConductivity':

                leg = np.log10(abs(1./leg))
                rho = np.hstack([rho,leg])

-            elif dtype == 'appr':
+            elif unitType == 'appResistivity':

                leg = np.log10(abs(leg))
                rho = np.hstack([rho,leg])

            else:
-                print """dtype must be 'appr' | 'appc' | 'volt' """
+                print """unitType must be 'appResistivity' | 'appConductivity' | 'volt' """
                break

        midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
@@ -259,7 +255,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl
    # Grid points
    grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
    grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear')
-    
+
    # Scale the color scheme
    if clim == None:
        vmin, vmax = rho.min(), rho.max()
@@ -268,36 +264,37 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl

    # Plot data
    grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho)
+
+
+
    ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, vmin = vmin, vmax = vmax)
    plt.gca().tick_params(axis='both', which='major', labelsize=8)
-    
+
    if contour is not None:
        plt.contour(grid_x,grid_z,grid_rho,levels = contour,colors = 'r', vmin = vmin, vmax = vmax)
-        
+
    # Add scatter points
    axs.scatter(midx,midz,s=10,c=rho.T, vmin = vmin, vmax = vmax)
-    
+
    if colorbar:
-        
-        if dtype == 'volt':
+
+        if unitType == 'volt':
            cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal")

-        else:               
+        else:
            cbar = plt.colorbar(ph, ax = axs, format="$10^{%.1f}$",fraction=0.04,orientation="horizontal")
-        
-        cmin,cmax = cbar.get_clim()
-        ticks = np.linspace(cmin,cmax,3)
-        cbar.set_ticks(ticks)
-        cbar.ax.tick_params(labelsize=10)
-        
-        if cblabel:
-            if dtype == 'appc':
-                cbar.set_label("App.Cond",size=12)
-            elif dtype == 'appr':
-                cbar.set_label("App.Res.",size=12)
-            elif dtype == 'volt':
-                cbar.set_label("Potential (V)",size=12)

+            cmin,cmax = cbar.get_clim()
+            ticks = np.linspace(cmin,cmax,3)
+            cbar.set_ticks(ticks)
+            cbar.ax.tick_params(labelsize=10)
+
+            if unitType == 'appConductivity':
+                cbar.set_label("App.Cond",size=12)
+            elif unitType == 'appResistivity':
+                cbar.set_label("App.Res.",size=12)
+            elif unitType == 'volt':
+                cbar.set_label("Potential (V)",size=12)


    if not axlabel:
@@ -310,27 +307,24 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl

    return ph

-def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
+def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
    """
-        Load in endpoints and survey specifications to generate Tx, Rx location
-        stations.
+    Load in endpoints and survey specifications to generate Tx, Rx location
+    stations.

-        Assumes flat topo for now...
+    Assumes flat topo for now...

-        Input:
-        :param endl -> input endpoints [x1, y1, z1, x2, y2, z2]
-        :object mesh -> SimPEG mesh object
-        :switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient'
-        : param a, n -> pole seperation, number of rx dipoles per tx
+    :param numpy.array endl: input endpoints [[x1, y1] , [x2, y2]]
+    :param Mesh mesh: SimPEG mesh object
+    :param string surveyType: 'dipole-dipole' | 'pole-dipole' | 'gradient'
+    :param float AM_sep: transmitter (A) - receiver (M) seperation
+    :param float b: receiver dipole seperation
+    :param float nrx: pole seperation, number of rx dipoles per tx

-        Output:
-        :param Tx, Rx -> List objects for each tx location
-            Lines: P1x, P1y, P1z, P2x, P2y, P2z
+    :rtype: DC.Survey, Src, Rx
+    :returns: DC survey, Source

-        Created on Wed December 9th, 2015
-
-        @author: dominiquef
-        !! Require clean up to deal with DCsurvey
+    !! Require clean up to deal with DCsurvey
    """

    from SimPEG import np
@@ -346,17 +340,17 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
    dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
    dl_y = ( endl[1,1] - endl[0,1] ) / dl_len

-    nstn = np.floor( dl_len / a )
+    nstn = np.floor( dl_len / AM_sep )

    # Compute discrete pole location along line
-    stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
-    stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
+    stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*AM_sep
+    stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*AM_sep

    # Create line of P1 locations
    M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]

    # Create line of P2 locations
-    N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
+    N = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]

    ## Build list of Tx-Rx locations depending on survey type
    # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn]
@@ -366,14 +360,14 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
    SrcList = []


-    if stype != 'gradient':
+    if surveyType != 'gradient':

        for ii in range(0, int(nstn)-1):


-            if stype == 'dpdp':
+            if surveyType == 'dipole-dipole':
                tx = np.c_[M[ii,:],N[ii,:]]
-            elif stype == 'pdp':
+            elif surveyType == 'pole-dipole':
                tx = np.c_[M[ii,:],M[ii,:]]

            # Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]])
@@ -382,33 +376,33 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
            AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])

            # Number of receivers to fit
-            nstn = np.min([np.floor( (AB - b) / a ) , n])
+            nstn = np.min([np.floor( (AB - MN_sep) / AM_sep ) , nrx])

            # Check if there is enough space, else break the loop
            if nstn <= 0:
                continue

            # Compute discrete pole location along line
-            stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
-            stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
+            stn_x = N[ii,0] + dl_x*MN_sep + np.array(range(int(nstn)))*dl_x*AM_sep
+            stn_y = N[ii,1] + dl_y*MN_sep + np.array(range(int(nstn)))*dl_y*AM_sep

            # Create receiver poles
            # Create line of P1 locations
            P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]

            # Create line of P2 locations
-            P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
+            P2 = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]

            Rx.append(np.c_[P1,P2])
            rxClass = DC.RxDipole(P1, P2)
            Tx.append(tx)
-            if stype == 'dpdp':
+            if surveyType == 'dipole-dipole':
                srcClass = DC.SrcDipole([rxClass], M[ii,:],N[ii,:])
-            elif stype == 'pdp':
+            elif surveyType == 'pole-dipole':
                srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:])
            SrcList.append(srcClass)

-    elif stype == 'gradient':
+    elif surveyType == 'gradient':

        # Gradient survey only requires Tx at end of line and creates a square
        # grid of receivers at in the middle at a pre-set minimum distance
@@ -416,23 +410,23 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
        Tx.append(np.c_[M[0,:],N[-1,:]])

        # Get the edge limit of survey area
-        min_x = endl[0,0] + dl_x * b
-        min_y = endl[0,1] + dl_y * b
+        min_x = endl[0,0] + dl_x * MN_sep
+        min_y = endl[0,1] + dl_y * MN_sep

-        max_x = endl[1,0] - dl_x * b
-        max_y = endl[1,1] - dl_y * b
+        max_x = endl[1,0] - dl_x * MN_sep
+        max_y = endl[1,1] - dl_y * MN_sep

        box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
        box_w = box_l/2.

-        nstn = np.floor( box_l / a )
+        nstn = np.floor( box_l / AM_sep )

        # Compute discrete pole location along line
-        stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
-        stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
+        stn_x = min_x + np.array(range(int(nstn)))*dl_x*AM_sep
+        stn_y = min_y + np.array(range(int(nstn)))*dl_y*AM_sep

        # Define number of cross lines
-        nlin = int(np.floor( box_w / a ))
+        nlin = int(np.floor( box_w / AM_sep ))
        lind = range(-nlin,nlin+1)

        ngrad = nstn * len(lind)
@@ -441,12 +435,12 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
        for ii in range( len(lind) ):

            # Move line in perpendicular direction by dipole spacing
-            lxx = stn_x - lind[ii]*a*dl_y
-            lyy = stn_y + lind[ii]*a*dl_x
+            lxx = stn_x - lind[ii]*AM_sep*dl_y
+            lyy = stn_y + lind[ii]*AM_sep*dl_x


            M = np.c_[ lxx, lyy , np.ones(nstn).T*mesh.vectorNz[-1]]
-            N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
+            N = np.c_[ lxx+AM_sep*dl_x, lyy+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]

            rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N]

@@ -455,44 +449,38 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
        srcClass = DC.SrcDipole([rxClass], M[0,:], N[-1,:])
        SrcList.append(srcClass)
    else:
-        print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
+        print """surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """

    survey = DC.SurveyDC(SrcList)
    return survey, Tx, Rx

-def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
+
+def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
    """
        Write UBC GIF DCIP 2D or 3D observation file

-        Input:
-        :string     fileName -> including path where the file is written out
-        :DCsurvey   DC survey class object
-        :string     dtype ->  either '2D' | '3D'
-        :string     stype ->  either 'SURFACE' | 'GENERAL'
-
-        Output:
-        :param UBC2D-Data file
-        :return
-
-        Last edit: February 16th, 2016
-
-        @author: dominiquef
-
+        :param string fileName: including path where the file is written out
+        :param Survey DCsurvey: DC survey class object
+        :param string dim:  either '2D' | '3D'
+        :param string surveyType:  either 'SURFACE' | 'GENERAL'
+        :rtype: file
+        :return: UBC2D-Data file
    """
+
    from SimPEG import mkvc

-    assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'"
-    assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
-    
+    assert (dim=='2D') | (dim=='3D'), "Data must be either '2D' | '3D'"
+    assert (surveyType=='SURFACE') | (surveyType=='GENERAL') | (surveyType=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
+
    fid = open(fileName,'w')
-    
-    
+    fid.write('! ' + surveyType + ' FORMAT\n')
+
    if iptype!=0:
        fid.write('IPTYPE=%i\n'%iptype)
-        
+
    else:
        fid.write('! ' + stype + ' FORMAT\n')
-        
+
    count = 0

    for ii in range(DCsurvey.nSrc):
@@ -506,33 +494,33 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
        M = rx[0]
        N = rx[1]

-        # Adapt source-receiver location for dtype and stype
-        if dtype=='2D':
+        # Adapt source-receiver location for dim and surveyType
+        if dim=='2D':

-            if stype == 'SIMPLE':
+            if surveyType == 'SIMPLE':

                #fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
                A = np.repeat(tx[0,0],M.shape[0],axis=0)
                B = np.repeat(tx[0,1],M.shape[0],axis=0)
                M = M[:,0]
                N = N[:,0]
-                
+
                np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')


            else:

-                if stype == 'SURFACE':
+                if surveyType == 'SURFACE':

                    fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
                    M = M[:,0]
                    N = N[:,0]

-                if stype == 'GENERAL':
+                if surveyType == 'GENERAL':

                    # Flip sign for z-elevation to depth
                    tx[2::2,:] = -tx[2::2,:]
-                    
+
                    fid.writelines("%e " % ii for ii in mkvc(tx[::2,:]))
                    M = M[:,0::2]
                    N = N[:,0::2]
@@ -540,31 +528,31 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
                    # Flip sign for z-elevation to depth
                    M[:,1::2] = -M[:,1::2]
                    N[:,1::2] = -N[:,1::2]
-                    
+
                fid.write('%i\n'% nD)
                np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')

-        if dtype=='3D':
+        if dim=='3D':

-            if stype == 'SURFACE':
+            if surveyType == 'SURFACE':

                fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:]))
                M = M[:,0:2]
                N = N[:,0:2]

-            if stype == 'GENERAL':
+            if surveyType == 'GENERAL':

                fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))

            fid.write('%i\n'% nD)
            np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
            fid.write('\n')
-            
+
        count += nD

    fid.close()

-def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
+def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'):
    """
        Read DC survey and projects the coordinate system
        according to the flag = 'Xloc' | 'Yloc' | 'local' (default)
@@ -573,15 +561,9 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):

        The Z value is preserved, but Y coordinates zeroed.

-        Input:
-        :param survey3D
-
-        Output:
-        :figure survey2D
-
-        Edited April 6th, 2016
-
-        @author: dominiquef
+        :param DC.Survey survey3D: 3D simpeg DC survey
+        :rtype: DC.Survey
+        :return: survey2D

    """
    from SimPEG import np
@@ -666,39 +648,34 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
    DCsurvey2D.std = np.asarray(DCsurvey.std)

    return DCsurvey2D
-    
-def readUBC_DC3Dobs(fileName, dtype = 'DC'):
+
+def readUBC_DC3Dobs(fileName, rtype = 'DC'):
    """
        Read UBC GIF IP 3D observation file and generate survey

-        Input:
-        :param fileName, path to the UBC GIF 3D obs file
-
-        Output:
-        :param IPsurvey
-        :return
-
-        @author: dominiquef
+        :param string fileName:, path to the UBC GIF 3D obs file
+        :rtype: Survey
+        :return: DCIPsurvey

    """
    zflag = True # Flag for z value provided
-    
+
    # Load file
-    if dtype == 'IP':
+    if rtype == 'IP':
        obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='IPTYPE')
-        
-    elif dtype == 'DC':
+
+    elif rtype == 'DC':
        obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
-        
+
    else:
-        print "dtype must be 'DC'(default) | 'IP'"
-        
+        print "rtype must be 'DC'(default) | 'IP'"
+
    # Pre-allocate
    srcLists = []
    Rx = []
    d = []
    wd = []
-    
+

    # Countdown for number of obs/tx
    count = 0
@@ -717,7 +694,7 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'):
            # Check if z value is provided, if False -> nan
            if len(temp)==5:
                tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan]
-                
+
                zflag = False # Pass on the flag to the receiver loc

            else:
@@ -729,12 +706,12 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'):
        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string

        # Filter out negative IP
-#        if temp[-2] < 0:       
+#        if temp[-2] < 0:
 #            count = count -1
 #            print "Negative!"
-#            
+#
 #        else:
-            
+
        # If the Z-location is provided, otherwise put nan
        if zflag:

@@ -772,17 +749,9 @@ def readUBC_DC2Dobs(fileName):
        ------- NEEDS TO BE UPDATED ------
        Read UBC GIF 2D observation file and generate arrays for tx-rx location

-        Input:
-        :param fileName, path to the UBC GIF 2D model file
-
-        Output:
-        :param rx, tx
-        :return
-
-        Created on Thu Nov 12 13:14:10 2015
-
-        @author: dominiquef
-
+        :param string fileName: path to the UBC GIF 2D model file
+        :rtype: (DC.Src, DC.Rx, ??, ??)
+        :return: source_locs, rx_locs, ??, ??
    """

    from SimPEG import np
@@ -822,11 +791,9 @@ def readUBC_DC2Dpre(fileName):
        Read UBC GIF DCIP 2D observation file and generate arrays for tx-rx location

        Input:
-        :param fileName, path to the UBC GIF 3D obs file
-
-        Output:
-        DCsurvey
-        :return
+        :param string fileName: path to the UBC GIF 3D obs file
+        :rtype: DC.Survey
+        :return: DCsurvey

        Created on Mon March 9th, 2016 << Doug's 70th Birthday !! >>

@@ -888,12 +855,9 @@ def readUBC_DC2DMesh(fileName):
    """
        Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg

-        Input:
-        :param fileName, path to the UBC GIF mesh file
-
-        Output:
-        :param SimPEG TensorMesh 2D object
-        :return
+        :param string fileName: path to the UBC GIF mesh file
+        :rtype: Mesh.TensorMesh
+        :return: SimPEG TensorMesh 2D object

        Created on Thu Nov 12 13:14:10 2015

@@ -959,12 +923,9 @@ def xy_2_lineID(DCsurvey):
        they were collected. May need to generalize for random
        point locations, but will be more expensive

-        Input:
-        :param DCdict Vectors of station location
-
-        Output:
-        :param LineID Vector of integers
-        :return
+        :param numpy.array DCdict: Vectors of station location
+        :rtype: numpy.array
+        :return: LineID Vector of integers

        Created on Thu Feb 11, 2015

@@ -144,6 +144,7 @@ class BetaSchedule(InversionDirective):
            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
            self.invProb.beta /= self.coolingFactor

+
 class TargetMisfit(InversionDirective):

    chifact = 1.
@@ -166,7 +167,7 @@ class TargetMisfit(InversionDirective):



-class _SaveEveryIteration(InversionDirective):
+class SaveEveryIteration(InversionDirective):
    @property
    def name(self):
        if getattr(self, '_name', None) is None:
@@ -187,7 +188,7 @@ class _SaveEveryIteration(InversionDirective):
        self._fileName = value


-class SaveModelEveryIteration(_SaveEveryIteration):
+class SaveModelEveryIteration(SaveEveryIteration):
    """SaveModelEveryIteration"""

    def initialize(self):
@@ -197,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):
@@ -211,7 +212,7 @@ 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):
+class SaveOutputDictEveryIteration(SaveEveryIteration):
    """SaveOutputDictEveryIteration"""

    def initialize(self):
@@ -242,12 +243,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
        # Save the file as a npz
        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)

-
-# class UpdateReferenceModel(Parameter):
-
-#     mref0 = None
-
-#     def nextIter(self):
 #         mref = getattr(self, 'm_prev', None)
 #         if mref is None:
 #             if self.debug: print 'UpdateReferenceModel is using mref0'
@@ -258,56 +253,151 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 class Update_IRLS(InversionDirective):

    eps_min = None
+    eps = None
+    norms = [2.,2.,2.,2.]
    factor = None
    gamma = None
    phi_m_last = None
    phi_d_last = None
+    f_old = None
+    f_min_change = 1e-2
+    beta_tol = 5e-2
+    prctile = 95

+    # Solving parameter for IRLS (mode:2)
+    IRLSiter   = 0
+    minGNiter = 5
+    maxIRLSiter = 10
+    iterStart = 0
+
+    # Beta schedule
+    coolingFactor = 2.
+    coolingRate = 1
+
+    mode = 1
+
+    @property
+    def target(self):
+        if getattr(self, '_target', None) is None:
+            self._target = self.survey.nD*0.5
+        return self._target
+    @target.setter
+    def target(self, val):
+        self._target = val

    def initialize(self):

-        # Scale the regularization for changes in norm
-        if getattr(self, 'phi_m_last', None) is not None:
-
-            self.reg.curModel = self.invProb.curModel
-            self.reg.gamma = 1.
-            phim_new = self.reg.eval(self.invProb.curModel)
-            self.gamma = self.phi_m_last / phim_new
-
-        self.reg.curModel = self.invProb.curModel
-        self.reg.gamma = self.gamma
-
-        if getattr(self, 'phi_d_last', None) is None:
-            self.phi_d_last = self.invProb.phi_d
+        if self.mode == 1:
+            self.reg.norms = [2., 2., 2., 2.]

    def endIter(self):
-        # Cool the threshold parameter if required
-        if getattr(self, 'factor', None) is not None:
-            eps = self.reg.eps / self.factor

-            if getattr(self, 'eps_min', None) is not None:
-                self.reg.eps = np.max([self.eps_min,eps])
+        # After reaching target misfit with l2-norm, switch to IRLS (mode:2)
+        if self.invProb.phi_d < self.target and self.mode == 1:
+            print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
+
+            self.mode = 2
+
+            # 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
+            self.iterStart = self.opt.iter
+            self.phi_d_last = self.invProb.phi_d
+            self.phi_m_last = self.invProb.phi_m_last
+
+            self.reg.l2model = self.invProb.curModel
+            self.reg.curModel = self.invProb.curModel
+
+            if getattr(self, 'f_old', None) is None:
+                self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False)
+
+        # Beta Schedule
+        if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
+            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
+            self.invProb.beta /= self.coolingFactor
+
+
+        # Only update after GN iterations
+        if (self.opt.iter-self.iterStart) % self.minGNiter == 0 and self.mode==2:
+
+            self.IRLSiter += 1
+
+            phim_new = self.reg.eval(self.invProb.curModel)
+            self.f_change = np.abs(self.f_old - phim_new) / self.f_old
+
+            print "Regularization decrease: %6.3e" % (self.f_change)
+
+            # Check for maximum number of IRLS cycles
+            if self.IRLSiter == self.maxIRLSiter:
+                print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
+                self.opt.stopNextIteration = True
+                return
+
+            # Check if the function has changed enough
+            if self.f_change < self.f_min_change and self.IRLSiter > 1:
+                print "Minimum decrease in regularization. End of IRLS"
+                self.opt.stopNextIteration = True
+                return
            else:
-                self.reg.eps = eps
+                self.f_old = phim_new

-        # Get phi_m at the end of current iteration
-        self.phi_m_last = self.invProb.phi_m_last
+#            # Cool the threshold parameter if required
+#            if getattr(self, 'factor', None) is not None:
+#                eps = self.reg.eps / self.factor
+#
+#                if getattr(self, 'eps_min', None) is not None:
+#                    self.reg.eps = np.max([self.eps_min,eps])
+#                else:
+#                    self.reg.eps = eps

-         # Update the model used for the IRLS weights
-        self.reg.curModel = self.invProb.curModel
+            # Get phi_m at the end of current iteration
+            self.phi_m_last = self.invProb.phi_m_last

-        # Temporarely set gamma to 1. to get raw phi_m
-        self.reg.gamma = 1.
+            # Reset the regularization matrices so that it is
+            # recalculated for current model
+            self.reg._Wsmall = None
+            self.reg._Wx = None
+            self.reg._Wy = None
+            self.reg._Wz = None

-        # Compute new model objective function value
-        phim_new = self.reg.eval(self.invProb.curModel)
+             # Update the model used for the IRLS weights
+            self.reg.curModel = self.invProb.curModel

-        # Update gamma to scale the regularization between IRLS iterations
-        self.reg.gamma = self.phi_m_last / phim_new
+            # Temporarely set gamma to 1. to get raw phi_m
+            self.reg.gamma = 1.

-        # Set the weighting matrix to None so that it is recomputed next time
-        # it is called in the inversion
-        self.reg._W = None
+            # Compute new model objective function value
+            phim_new = self.reg.eval(self.invProb.curModel)
+
+            # Update gamma to scale the regularization between IRLS iterations
+            self.reg.gamma = self.phi_m_last / phim_new
+
+            # Reset the regularization matrices again for new gamma
+            self.reg._Wsmall = None
+            self.reg._Wx = None
+            self.reg._Wy = None
+            self.reg._Wz = None
+
+            # Check if misfit is within the tolerance, otherwise scale beta
+            val = self.invProb.phi_d / (self.survey.nD*0.5)
+
+            if np.abs(1.-val) > self.beta_tol:
+                self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d

 class Update_lin_PreCond(InversionDirective):
    """
@@ -360,19 +450,3 @@ class Update_Wj(InversionDirective):
            JtJdiag = JtJdiag / max(JtJdiag)

            self.reg.wght = JtJdiag
-
-class Scale_Beta(InversionDirective):
-    """
-        Instead of a linear cooling schedule, beta is allowed to change based
-        on the ratio between the target misfit and the current data misfit. The
-        update is done only if the misfit is outside some threshold bounds.
-    """
-    tol = 0.05
-
-    def endIter(self):
-
-        # Check if misfit is within the tolerance, otherwise adjust beta
-        val = self.invProb.phi_d / (self.survey.nD*0.5)
-
-        if np.abs(1.-val) > self.tol:
-            self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
@@ -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
+
+
+
+
+
@@ -2,3 +2,4 @@ from TDEM import hzAnalyticDipoleT
 from FDEM import hzAnalyticDipoleF
 from FDEMcasing import *
 from DC import DCAnalyticHalf, DCAnalyticSphere
+from FDEMDipolarfields import *
@@ -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 = {}
@@ -169,9 +169,7 @@ class BaseEMProblem(Problem.BaseProblem):

        dMeSigmaI_dI = -self.MeSigmaI**2
        dMe_dsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u)
-        dsig_dm = self.curModel.sigmaDeriv
-        return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm))
-        # return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u)
+        return dMeSigmaI_dI * ( dMe_dsig * self.curModel.sigmaDeriv )

    @property
    def MfRho(self):
@@ -187,8 +185,7 @@ class BaseEMProblem(Problem.BaseProblem):
        """
            Derivative of :code:`MfRho` with respect to the model.
        """
-        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv)
-        # self.curModel.rhoDeriv
+        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv

    @property
    def MfRhoI(self):
@@ -208,9 +205,7 @@ class BaseEMProblem(Problem.BaseProblem):

        dMfRhoI_dI = -self.MfRhoI**2
        dMf_drho = self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u)
-        return dMfRhoI_dI * ( dMf_drho * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) )
-
-        # return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
+        return dMfRhoI_dI * ( dMf_drho * self.curModel.rhoDeriv )

 class BaseEMSurvey(Survey.BaseSurvey):

@@ -222,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
@@ -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
@@ -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
@@ -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
@@ -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
        """
@@ -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)

@@ -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
@@ -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()
@@ -2,7 +2,7 @@ from SimPEG import Mesh, Utils, np, sp
 import SimPEG.DCIP as DC
 import time

-def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', plotIt=True):
+def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', unitType='appConductivity', plotIt=True):
    """
        DC Forward Simulation
        =====================
@@ -15,14 +15,14 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
        loc     = Location of spheres [[x1,y1,z1],[x2,y2,z2]]
        radi    = Radius of spheres [r1,r2]
        param   = Conductivity of background and two spheres [m0,m1,m2]
-        stype   = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
-        dtype   = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
+        surveyType = survey type 'pole-dipole' or 'dipole-dipole'
+        unitType = Data type "appResistivity" | "appConductivity"  | "volt"
        Created by @fourndo

    """

-    assert stype in ['pdp', 'dpdp'], "Source type (stype) must be pdp or dpdp (pole dipole or dipole dipole)"
-    assert dtype in ['appr', 'appc', 'volt'], "Data type (dtype) must be appr (app res) or appc (app cond) or volt (potential)"
+    assert surveyType in ['pole-dipole', 'dipole-dipole'], "Source type (surveyType) must be pdp or dpdp (pole dipole or dipole dipole)"
+    assert unitType in ['appResistivity', 'appConductivity', 'volt'], "Unit type (unitType) must be appResistivity or appConductivity or volt (potential)"

    if loc is None:
        loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]]
@@ -73,8 +73,8 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
    locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]]

    # We will handle the geometry of the survey for you and create all the combination of tx-rx along line
-    # [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2])
-    survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2])
+    # [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])
+    survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])

    # Define some global geometry
    dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )
@@ -118,8 +118,8 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
        rxloc_N = np.asarray(Rx[ii][:,3:])


-        # For usual cases "dpdp" or "gradient"
-        if stype == 'pdp':
+        # For usual cases 'dipole-dipole' or "gradient"
+        if surveyType == 'pole-dipole':
            # Create an "inifinity" pole
            tx =  np.squeeze(Tx[ii][:,0:1])
            tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
@@ -157,12 +157,12 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
        fig = plt.figure(figsize=(7,7))
        ax = plt.subplot(2,1,1, aspect='equal')
        # Plot the location of the spheres for reference
-        circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
-        circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
+        circle1=plt.Circle((loc[0,0], loc[2,0]), radi[0], color='w', fill=False, lw=3)
+        circle2=plt.Circle((loc[0,1], loc[2,1]), radi[1], color='k', fill=False, lw=3)
        ax.add_artist(circle1)
        ax.add_artist(circle2)

-        dat = mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y',
+        dat = mesh.plotSlice(np.log10(model), ax = ax, normal = 'Y',
                             ind = indy,grid=True, clim = np.log10([sig.min(),sig.max()]))

        ax.set_title('3-D model')
@@ -188,15 +188,13 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
        ax2 = plt.subplot(2,1,2, aspect='equal')

        # Plot the location of the spheres for reference
-        circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
-        circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
+        circle1=plt.Circle((loc[0,0], loc[2,0]), radi[0], color='w', fill=False, lw=3)
+        circle2=plt.Circle((loc[0,1], loc[2,1]), radi[1], color='k', fill=False, lw=3)
        ax2.add_artist(circle1)
        ax2.add_artist(circle2)

        # Add the speudo section
-        dat = DC.plot_pseudoSection(survey2D,ax2,stype=stype, dtype = dtype)
-
-        # plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
+        dat = DC.plot_pseudoSection(survey2D, ax2, surveyType=surveyType, unitType=unitType)        # plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
        # plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y')
        # plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k')
        ax2.set_title('Apparent Conductivity data')
@@ -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 ---------------------------
@@ -1,7 +1,7 @@
 from SimPEG import *


-def run(N=200, plotIt=True):
+def run(N=100, plotIt=True):
    """
        Inversion: Linear Problem
        =========================
@@ -18,6 +18,8 @@ def run(N=200, plotIt=True):
    mesh = Mesh.TensorMesh([N])

    m0 = np.ones(mesh.nC) * 1e-4
+    mref = np.zeros(mesh.nC)
+
    nk = 10
    jk = np.linspace(1.,nk,nk)
    p = -2.
@@ -40,67 +42,35 @@ def run(N=200, 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.wght = wr
-
    dmis = DataMisfit.l2_DataMisfit(survey)
    dmis.Wd = 1./wd

-    opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
-    invProb.curModel = m0
-
-    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
-    target = Directives.TargetMisfit()
-
    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

-#==============================================================================
-#     fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-#     dmdx = reg.mesh.cellDiffxStencil * mrec
-#     plt.plot(np.sort(dmdx))
-#==============================================================================
-
-    #reg.recModel = mrec
-    reg.wght = np.ones(mesh.nC)
    reg.mref = np.zeros(mesh.nC)
-    reg.eps_p = 5e-2
-    reg.eps_q = 1e-2
-    reg.norms   = [0., 0., 2., 2.]
-    reg.wght = wr
+    

-    opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.)
-    beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
-    #betaest = Directives.BetaEstimate_ByEig()
-    target = Directives.TargetMisfit()
-    IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
+    opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
+    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
+    update_Jacobi = Directives.Update_lin_PreCond()
+    
+    # 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=[beta,IRLS])
-
-    m0 = mrec
+    inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])

    # Run inversion
    mrec = inv.run(m0)
@@ -117,7 +87,7 @@ def run(N=200, plotIt=True):
        axes[0].set_title('Columns of matrix G')

        axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
-        axes[1].plot(mesh.vectorCCx, ml2, 'r-')
+        axes[1].plot(mesh.vectorCCx, reg.l2model, 'r-')
        #axes[1].legend(('True Model', 'Recovered Model'))
        axes[1].set_ylim(-1.0,1.25)

@@ -7,7 +7,7 @@ import matplotlib.pyplot as plt
 def run(plotIt=True):
    """
        MT: 1D: Inversion
-        =======================
+        =================

        Forward model 1D MT data.
        Setup and run a MT 1D inversion.
@@ -50,7 +50,7 @@ def run(plotIt=True):
    m_0 = np.log(sigma_0[active])

    # Set the mapping
-    actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
+    actMap = simpeg.Maps.InjectActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
    mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap

    ## Setup the layout of the survey, set the sources and the connected receivers
@@ -76,7 +76,7 @@ def run(plotIt=True):
    survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)

    if plotIt:
-        fig = MT.Utils.dataUtils.plotMT1DModelData(problem)
+        fig = MT.Utils.dataUtils.plotMT1DModelData(problem, [m_0])
        fig.suptitle('Target - smooth true')


@@ -12,7 +12,7 @@ except:
 def run(plotIt=True, nFreq=1):
    """
        MT: 3D: Forward
-        =======================
+        ===============

        Forward model 3D MT data.

@@ -46,16 +46,15 @@ def run(plotIt=True, nFreq=1):
    survey = MT.Survey(srcList)

    ## Setup the problem object
-    problem = MT.Problem3D.eForm_ps(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 = MT.Data(survey,dataVec)
+    mtData = MT.Data(survey, dataVec)
    # Add plots
    if plotIt:
        pass
@@ -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()
+
@@ -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()
+
@@ -1,22 +1,25 @@
 from SimPEG import Mesh, Utils, np, SolverLU

-## 2D DC forward modeling example with Tensor and Curvilinear Meshes
-
 def run(plotIt=True):
+
+    """
+        Mesh: Basic Forward 2D DC Resistivity
+        =====================================
+
+        2D DC forward modeling example with Tensor and Curvilinear Meshes
+    """
+
    # Step1: Generate Tensor and Curvilinear Mesh
    sz = [40,40]
-    # Tensor Mesh
    tM = Mesh.TensorMesh(sz)
-    # Curvilinear Mesh
    rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
+
    # Step2: Direct Current (DC) operator
    def DCfun(mesh, pts):
        D = mesh.faceDiv
-        G = D.T
        sigma = 1e-2*np.ones(mesh.nC)
-        Msigi = mesh.getFaceInnerProduct(1./sigma)
-        MsigI = Utils.sdInv(Msigi)
-        A = D*MsigI*G
+        MsigI = mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True)
+        A = -D*MsigI*D.T
        A[-1,-1] /= mesh.vol[-1] # Remove null space
        rhs = np.zeros(mesh.nC)
        txind = Utils.meshutils.closestPoints(mesh, pts)
@@ -37,39 +40,17 @@ def run(plotIt=True):
    if not plotIt: return

    import matplotlib.pyplot as plt
-    import matplotlib
-    from matplotlib.mlab import griddata

    #Step4: Making Figure
    fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-    label = ["(a)", "(b)"]
-    opts = {}
    vmin, vmax = phitM.min(), phitM.max()
    dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True)
-
-    #TODO: At the moment Curvilinear Mesh do not have plotimage
-
-    Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F')
-    Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F')
-    PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear')
-    axes[1].contourf(Xi, Yi, PHIrM, 100, vmin=vmin, vmax=vmax)
-
+    dat = rM.plotImage(phirM, ax=axes[1], clim=(vmin, vmax), grid=True)
    cb = plt.colorbar(dat[0], ax=axes[0]); cb.set_label("Voltage (V)")
    cb = plt.colorbar(dat[0], ax=axes[1]); cb.set_label("Voltage (V)")

-    tM.plotGrid(ax=axes[0], **opts)
    axes[0].set_title('TensorMesh')
-    rM.plotGrid(ax=axes[1], **opts)
    axes[1].set_title('CurvilinearMesh')
-    for i in range(2):
-        axes[i].set_xlim(0.025, 0.975)
-        axes[i].set_ylim(0.025, 0.975)
-        axes[i].text(0., 1.0, label[i], fontsize=20)
-        if i==0:
-            axes[i].set_ylabel("y")
-        else:
-            axes[i].set_ylabel(" ")
-        axes[i].set_xlabel("x")
    plt.show()


@@ -0,0 +1,43 @@
+from SimPEG import *
+from SimPEG.Utils import surface2ind_topo
+
+
+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.
+
+    """
+
+    mesh = Mesh.TensorMesh([nx,ny], x0='CC') # 2D mesh
+    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
+
+    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))
+        mesh.plotGrid(ax=ax, nodes=True, centers=True)
+        ax.plot(xtopo,topo,'k',linewidth=1)
+        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)
+        plt.show()
+
+
+if __name__ == '__main__':
+    run(plotIt=True)
@@ -8,9 +8,11 @@ import EM_FDEM_Analytic_MagDipoleWholespace
 import EM_Schenkel_Morrison_Casing
 import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
-import Forward_BasicDirectCurrent
 import Inversion_IRLS
 import Inversion_Linear
+import Maps_ComboMaps
+import Maps_Mesh2Mesh
+import Mesh_Basic_ForwardDC
 import Mesh_Basic_PlotImage
 import Mesh_Basic_Types
 import Mesh_Operators_CahnHilliard
@@ -20,8 +22,9 @@ import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
 import MT_1D_ForwardAndInversion
 import MT_3D_Foward
+import Utils_surface2ind_topo

-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
+__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 #####

@@ -37,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)

@@ -94,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')
@@ -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.
@@ -86,7 +86,7 @@ class polxy_1Dprimary(BaseMTSrc):
        Get the electrical field source
        """
        e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
+        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
        sigma_p = Map_sigma_p._transform(self.sigma1d)
        # Make mass matrix
        # Note: M(sig) - M(sig_p) = M(sig - sig_p)
@@ -163,7 +163,7 @@ class polxy_3Dprimary(BaseMTSrc):
        Get the electrical field source
        """
        e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
+        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
        sigma_p = Map_sigma_p._transform(self.sigma1d)
        # Make mass matrix
        # Note: M(sig) - M(sig_p) = M(sig - sig_p)
@@ -19,7 +19,7 @@ def getAppRes(MTdata):
        zList.append(zc)
    return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]

-def rotateData(MTdata,rotAngle):
+def rotateData(MTdata, rotAngle):
    '''
    Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
    '''
@@ -44,19 +44,19 @@ def rotateData(MTdata,rotAngle):
    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):
+def makeAnalyticSolution(mesh, model, elev, freqs):
    from SimPEG import MT
    data1D = []
    for freq in freqs:
@@ -70,7 +70,7 @@ 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):
+def plotMT1DModelData(problem, models, symList=None):
    from SimPEG import MT
    # Setup the figure
    fontSize = 15
@@ -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])
@@ -2,6 +2,7 @@ from SimPEG import Utils, np
 from BaseMesh import BaseRectangularMesh
 from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
+from View import CurvView

 # Some helper functions.
 length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
@@ -10,7 +11,7 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
 normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))


-class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
+class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView):
    """
    CurvilinearMesh is a mesh class that deals with curvilinear meshes.

@@ -330,102 +331,6 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):



-    #############################################
-    #            Plotting Functions             #
-    #############################################
-
-    def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True,  showIt=False):
-        """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
-
-
-        .. plot::
-            :include-source:
-
-            from SimPEG import Mesh, Utils
-            X, Y = Utils.exampleLrmGrid([3,3],'rotate')
-            M = Mesh.CurvilinearMesh([X, Y])
-            M.plotGrid(showIt=True)
-
-        """
-        import matplotlib.pyplot as plt
-        import matplotlib
-        from mpl_toolkits.mplot3d import Axes3D
-        mkvc = Utils.mkvc
-
-        axOpts = {'projection':'3d'} if self.dim == 3 else {}
-        if ax is None: ax = plt.subplot(111, **axOpts)
-
-        NN = self.r(self.gridN, 'N', 'N', 'M')
-        if self.dim == 2:
-
-            if lines:
-                X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
-                Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
-
-                X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
-                Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
-
-                X = np.r_[X1, X2]
-                Y = np.r_[Y1, Y2]
-
-                ax.plot(X, Y, 'b-')
-            if centers:
-                ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
-
-            # Nx = self.r(self.normals, 'F', 'Fx', 'V')
-            # Ny = self.r(self.normals, 'F', 'Fy', 'V')
-            # Tx = self.r(self.tangents, 'E', 'Ex', 'V')
-            # Ty = self.r(self.tangents, 'E', 'Ey', 'V')
-
-            # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
-
-            # nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
-            # ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
-            # ax.plot(nX, nY, 'r-')
-
-            # nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
-            # #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
-            # ax.plot(nX, nY, 'g-')
-
-            # tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
-            # tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
-            # ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
-            # ax.plot(tX, tY, 'r-')
-
-            # nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
-            # #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
-            # ax.plot(nX, nY, 'g-')
-
-        elif self.dim == 3:
-            X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
-            Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
-            Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
-
-            X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten()
-            Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten()
-            Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten()
-
-            X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten()
-            Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten()
-            Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten()
-
-            X = np.r_[X1, X2, X3]
-            Y = np.r_[Y1, Y2, Y3]
-            Z = np.r_[Z1, Z2, Z3]
-
-            ax.plot(X, Y, 'b', zs=Z)
-            ax.set_zlabel('x3')
-
-        ax.grid(True)
-        ax.set_xlabel('x1')
-        ax.set_ylabel('x2')
-
-        if showIt: plt.show()
-
-
 if __name__ == '__main__':
    nc = 5
    h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
@@ -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]
@@ -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"
@@ -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
@@ -162,12 +157,9 @@ class TensorMeshIO(object):
        """
            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()))
@@ -183,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
        """

@@ -201,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
@@ -231,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.
@@ -286,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
@@ -335,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:
@@ -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"
@@ -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::
@@ -552,7 +552,8 @@ class CurvView(object):
    def __init__(self):
        pass

-    def plotGrid(self, length=0.05, showIt=False):
+
+    def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True,  showIt=False):
        """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.


@@ -560,60 +561,63 @@ class CurvView(object):
            :include-source:

            from SimPEG import Mesh, Utils
-            X, Y = Utils.exampleCurvGird([3,3],'rotate')
+            X, Y = Utils.exampleLrmGrid([3,3],'rotate')
            M = Mesh.CurvilinearMesh([X, Y])
            M.plotGrid(showIt=True)

        """
+        import matplotlib.pyplot as plt
+        import matplotlib
+        from mpl_toolkits.mplot3d import Axes3D
+
+        axOpts = {'projection':'3d'} if self.dim == 3 else {}
+        if ax is None: ax = plt.subplot(111, **axOpts)
+
        NN = self.r(self.gridN, 'N', 'N', 'M')
        if self.dim == 2:
-            fig = plt.figure(2)
-            fig.clf()
-            ax = plt.subplot(111)
-            X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
-            Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()

-            X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
-            Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
+            if lines:
+                X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
+                Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()

-            X = np.r_[X1, X2]
-            Y = np.r_[Y1, Y2]
+                X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
+                Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()

-            plt.plot(X, Y)
+                X = np.r_[X1, X2]
+                Y = np.r_[Y1, Y2]

-            plt.hold(True)
-            Nx = self.r(self.normals, 'F', 'Fx', 'V')
-            Ny = self.r(self.normals, 'F', 'Fy', 'V')
-            Tx = self.r(self.tangents, 'E', 'Ex', 'V')
-            Ty = self.r(self.tangents, 'E', 'Ey', 'V')
+                ax.plot(X, Y, 'b-')
+            if centers:
+                ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')

-            plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
+            # Nx = self.r(self.normals, 'F', 'Fx', 'V')
+            # Ny = self.r(self.normals, 'F', 'Fy', 'V')
+            # Tx = self.r(self.tangents, 'E', 'Ex', 'V')
+            # Ty = self.r(self.tangents, 'E', 'Ey', 'V')

-            nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
-            nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
-            plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
-            plt.plot(nX, nY, 'r-')
+            # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')

-            nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
-            nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
-            #plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
-            plt.plot(nX, nY, 'g-')
+            # nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
+            # ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
+            # ax.plot(nX, nY, 'r-')

-            tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
-            tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
-            plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
-            plt.plot(tX, tY, 'r-')
+            # nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
+            # #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
+            # ax.plot(nX, nY, 'g-')

-            nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
-            nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
-            #plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
-            plt.plot(nX, nY, 'g-')
-            plt.axis('equal')
+            # tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
+            # tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
+            # ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
+            # ax.plot(tX, tY, 'r-')
+
+            # nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
+            # #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
+            # ax.plot(nX, nY, 'g-')

        elif self.dim == 3:
-            fig = plt.figure(3)
-            fig.clf()
-            ax = fig.add_subplot(111, projection='3d')
            X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
            Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
            Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
@@ -630,16 +634,50 @@ class CurvView(object):
            Y = np.r_[Y1, Y2, Y3]
            Z = np.r_[Z1, Z2, Z3]

-            plt.plot(X, Y, 'b', zs=Z)
+            ax.plot(X, Y, 'b', zs=Z)
            ax.set_zlabel('x3')

        ax.grid(True)
-        ax.hold(False)
        ax.set_xlabel('x1')
        ax.set_ylabel('x2')

        if showIt: plt.show()

+    def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
+        if self.dim == 3: raise NotImplementedError('This is not yet done!')
+
+        import matplotlib.pyplot as plt
+        import matplotlib
+        from mpl_toolkits.mplot3d import Axes3D
+        import matplotlib.colors as colors
+        import matplotlib.cm as cmx
+
+        if ax is None: ax = plt.subplot(111)
+        jet = cm = plt.get_cmap('jet')
+        cNorm  = colors.Normalize(
+            vmin=I.min() if clim is None else clim[0],
+            vmax=I.max() if clim is None else clim[1])
+
+        scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
+        # ax.set_xlim((self.x0[0], self.h[0].sum()))
+        # ax.set_ylim((self.x0[1], self.h[1].sum()))
+
+        Nx = self.r(self.gridN[:,0],'N','N','M')
+        Ny = self.r(self.gridN[:,1],'N','N','M')
+        cell = self.r(I,'CC','CC','M')
+
+        for ii in range(self.nCx):
+            for jj in range(self.nCy):
+                I = [ii,ii+1,ii+1,ii]
+                J = [jj,jj,jj+1,jj+1]
+                ax.add_patch(plt.Polygon(np.c_[Nx[I,J],Ny[I,J]], facecolor=scalarMap.to_rgba(cell[ii,jj]), edgecolor='k' if grid else 'none'))
+
+        scalarMap._A = []  # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
+        ax.set_xlabel('x')
+        ax.set_ylabel('y')
+        if showIt: plt.show()
+        return [scalarMap]
+

 if __name__ == '__main__':
    from SimPEG import *
@@ -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.'
@@ -1008,4 +1008,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
            indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
            delx[indx] = 0.

-            return delx
+            return delx
@@ -74,7 +74,7 @@ class Property(object):
                if linkedMap is None:
                    return None
                linkMap = linkMapClass(None) * linkedMap
-                m = getattr(self, '%s'%linkName)
+                m = getattr(self, '%sModel'%linkName)
                return linkMap.deriv( m )

            m = getattr(self, '%sModel'%prop.name)
@@ -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):
@@ -239,7 +239,7 @@ class PropMap(object):
            setattr(self, '%sMap'%name, mapping)
            setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP)))
            nP += mapping.nP
-        self.nP = nP 
+        self.nP = nP

    @property
    def defaultInvProp(self):
@@ -1,4 +1,6 @@
-import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
+import Utils, Maps, Mesh
+import numpy as np
+import scipy.sparse as sp

 class RegularizationMesh(object):
    """
@@ -8,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
    """

@@ -381,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:

@@ -403,7 +405,238 @@ class BaseRegularization(object):

        return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )

-class Tikhonov(BaseRegularization):
+class Simple(BaseRegularization):
+    """
+    Simple regularization that does not include length scales in the derivatives.
+    """
+
+    mrefInSmooth = False  #: include mref in the smoothness?
+    alpha_s      = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
+    alpha_x      = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'],     "Weight for the first derivative in the x direction")
+    alpha_y      = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'],     "Weight for the first derivative in the y direction")
+    alpha_z      = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'],     "Weight for the first derivative in the z direction")
+    cell_weights = 1.
+
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+
+        if isinstance(self.cell_weights,float):
+            self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
+
+    @property
+    def Wsmall(self):
+        """Regularization matrix Wsmall"""
+        if getattr(self,'_Wsmall', None) is None:
+            self._Wsmall = Utils.sdiag((self.alpha_s*self.cell_weights)**0.5)
+        return self._Wsmall
+
+    @property
+    def Wx(self):
+        """Regularization matrix Wx"""
+        if getattr(self, '_Wx', None) is None:
+            self._Wx = Utils.sdiag((self.alpha_x * (self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.regmesh.cellDiffxStencil
+        return self._Wx
+
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
+        if getattr(self, '_Wy', None) is None:
+            self._Wy = Utils.sdiag((self.alpha_y * (self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.regmesh.cellDiffyStencil
+        return self._Wy
+
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
+        if getattr(self, '_Wz', None) is None:
+            self._Wz = Utils.sdiag((self.alpha_z * (self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.regmesh.cellDiffzStencil
+        return self._Wz
+
+#    @property
+#    def Wsmooth(self):
+#        """Full smoothness regularization matrix W"""
+#        print 'wtf why are we using Wsmooth'
+#        raise NotImplementedError
+#        if getattr(self, '_Wsmooth', None) is None:
+#            wlist = (self.Wx,)
+#            if self.regmesh.dim > 1:
+#                wlist += (self.Wy,)
+#            if self.regmesh.dim > 2:
+#                wlist += (self.Wz,)
+#            self._Wsmooth = sp.vstack(wlist)
+#        return self._Wsmooth
+#
+#    @property
+#    def W(self):
+#        """Full regularization matrix W"""
+#        print 'wtf why are we using W'
+#        if getattr(self, '_W', None) is None:
+#            wlist = (self.Wsmall, self.Wx)
+#            if self.regmesh.dim > 1:
+#                wlist += (self.Wy,)
+#            if self.regmesh.dim > 2:
+#                wlist += (self.Wz,)
+#            self._W = sp.vstack(wlist)
+#        return self._W
+
+
+    @Utils.timeIt
+    def _evalSmall(self, m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmallDeriv(self, m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
+        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
+
+    @Utils.timeIt
+    def _evalSmall2Deriv(self, m, v = None):
+        rDeriv = self.Wsmall * ( self.mapping.deriv(m - self.mref) )
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothx(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wx * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wx * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothy(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wy * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wy * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothz(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wz * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wz * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmooth(self, m):
+        phiSmooth = self._evalSmoothx(m)
+        if self.regmesh.dim > 1:
+            phiSmooth += self._evalSmoothy(m)
+        if self.regmesh.dim > 2:
+            phiSmooth += self._evalSmoothz(m)
+        return phiSmooth
+
+    @Utils.timeIt
+    def _evalSmoothxDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wx * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wx * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wx * ( self.mapping * m )
+            return r.T * ( self.Wx * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothx2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wx * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wx * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothyDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wy * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wy * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wy * ( self.mapping * m )
+            return r.T * ( self.Wy * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothy2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wy * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wy * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothzDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wz * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wz * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wz * ( self.mapping * m )
+            return r.T * ( self.Wz * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothz2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wz * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wz * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothDeriv(self, m):
+        deriv = self._evalSmoothxDeriv(m)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyDeriv(m)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzDeriv(m)
+        return deriv
+
+    @Utils.timeIt
+    def _evalSmooth2Deriv(self, m, v=None):
+        deriv = self._evalSmoothx2Deriv(m, v)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothy2Deriv(m, v)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothz2Deriv(m, v)
+        return deriv
+
+
+    @Utils.timeIt
+    def eval(self, m):
+        return self._evalSmall(m) + self._evalSmooth(m)
+
+    @Utils.timeIt
+    def evalDeriv(self, m):
+        """
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+
+    @Utils.timeIt
+    def eval2Deriv(self, m, v=None):
+        return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v)
+
+
+
+class Tikhonov(Simple):
    """
    L2 Tikhonov regularization with both smallness and smoothness (first order
    derivative) contributions.
@@ -417,8 +650,8 @@ class Tikhonov(BaseRegularization):
    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
@@ -438,7 +671,7 @@ class Tikhonov(BaseRegularization):
    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
@@ -493,56 +726,131 @@ class Tikhonov(BaseRegularization):
            self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
        return self._Wzz

+
    @property
-    def Wsmooth(self):
+    def Wsmooth2(self):
        """Full smoothness regularization matrix W"""
        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx, self.Wxx)
+            wlist = (self.Wxx)
            if self.regmesh.dim > 1:
-                wlist += (self.Wy, self.Wyy)
+                wlist += (self.Wyy)
            if self.regmesh.dim > 2:
-                wlist += (self.Wz, self.Wzz)
+                wlist += (self.Wzz)
            self._Wsmooth = sp.vstack(wlist)
        return self._Wsmooth

-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
    @Utils.timeIt
-    def _evalSmall(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
+    def _evalSmoothxx(self, m):
        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
+            r = self.Wxx * ( self.mapping * (m - self.mref) )
        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * (m) )
+            r = self.Wxx * ( self.mapping * (m) )
        return 0.5 * r.dot(r)

+    @Utils.timeIt
+    def _evalSmoothyy(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wyy * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wyy * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothzz(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wzz * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wzz * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmooth2(self, m):
+        phiSmooth2 = self._evalSmoothxx(m)
+        if self.regmesh.dim > 1:
+            phiSmooth2 += self._evalSmoothyy(m)
+        if self.regmesh.dim > 2:
+            phiSmooth2 += self._evalSmoothzz(m)
+        return phiSmooth2
+
+    @Utils.timeIt
+    def _evalSmoothxxDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wxx * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wxx * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wxx * ( self.mapping * m )
+            return r.T * ( self.Wxx * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothyyDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wyy * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wyy * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wyy * ( self.mapping * m )
+            return r.T * ( self.Wyy * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothzzDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wzz * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wzz * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wzz * ( self.mapping * m )
+            return r.T * ( self.Wzz * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothxx2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wxx * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wxx * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothyy2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wyy * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wyy * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothzz2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wzz * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wzz * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothDeriv2(self, m):
+        deriv = self._evalSmoothxxDeriv(m)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyyDeriv(m)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzzDeriv(m)
+        return deriv
+
+    @Utils.timeIt
+    def _evalSmooth2Deriv2(self, m, v=None):
+        deriv = self._evalSmoothxx2Deriv(m, v)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyy2Deriv(m, v)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzz2Deriv(m, v)
+        return deriv
+
+
    @Utils.timeIt
    def eval(self, m):
-        return self._evalSmall(m) + self._evalSmooth(m)
-
-    @Utils.timeIt
-    def _evalSmallDeriv(self,m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
-
-    @Utils.timeIt
-    def _evalSmoothDeriv(self,m):
-        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * ( m - self.mref ) )
-            return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m )
-            return r.T * ( self.Wsmooth * self.mapping.deriv(m) )
+        return self._evalSmall(m) + self._evalSmooth(m) + self._evalSmooth2(m)

    @Utils.timeIt
    def evalDeriv(self, m):
@@ -560,184 +868,134 @@ class Tikhonov(BaseRegularization):
            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}

        """
-        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m)
+
+    def eval2Deriv(self, m, v=None):
+        """
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + self._evalSmooth2Deriv2(m, v)


-class Simple(Tikhonov):
+
+class Sparse(Simple):
    """
-    Simple regularization that does not include length scales in the derivatives.
+        The regularization is:
+    
+        .. math::
+    
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R W(m-m_\\text{ref})}
+    
+        where the IRLS weight
+    
+        .. math::
+    
+            R = \eta TO FINISH LATER!!!
+    
+        So the derivative is straight forward:
+    
+        .. math::
+    
+            R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})}
+    
+        The IRLS weights are recomputed after each beta solves.
+        It is strongly recommended to do a few Gauss-Newton iterations
+        before updating.
    """
-
-    mrefInSmooth = False  #: SMOOTH and SMOOTH_MOD_DIF options
-    alpha_s     = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
-    alpha_x     = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
-    alpha_y     = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
-    alpha_z     = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
-    wght        = 1.
+        
+    # set default values
+    eps_p = 1e-1        # Threshold value for the model norm
+    eps_q = 1e-1        # Threshold value for the model gradient norm
+    curModel = None     # Requires model to compute the weights
+    l2model = None
+    gamma = 1.          # Model norm scaling to smooth out convergence
+    norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz)
+    cell_weights = 1.        # Consider overwriting with sensitivity weights

    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)

-        if isinstance(self.wght,float):
-            self.wght = np.ones(self.regmesh.nC) * self.wght
+        if isinstance(self.cell_weights,float):
+            self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights

    @property
    def Wsmall(self):
        """Regularization matrix Wsmall"""
        if getattr(self,'_Wsmall', None) is None:
-            self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
+            if getattr(self, 'curModel', None) is None:
+                self.Rs = Utils.speye(self.regmesh.nC)
+
+            else:
+                f_m = self.mapping * (self.curModel - self.reg.mref)
+                self.rs = self.R(f_m , self.eps_p, self.norms[0])
+                self.Rs = Utils.sdiag( self.rs )
+
+            self._Wsmall = Utils.sdiag((self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs
+
        return self._Wsmall

    @property
    def Wx(self):
        """Regularization matrix Wx"""
-        if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil
+        if getattr(self,'_Wx', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel)
+                self.rx = self.R( f_m , self.eps_q, self.norms[1])
+                self.Rx = Utils.sdiag( self.rx )
+
+            self._Wx = Utils.sdiag(( self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
+
        return self._Wx

    @property
    def Wy(self):
        """Regularization matrix Wy"""
-        if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil
+        if getattr(self,'_Wy', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel)
+                self.ry = self.R( f_m , self.eps_q, self.norms[2])
+                self.Ry = Utils.sdiag( self.ry )
+
+            self._Wy = Utils.sdiag((self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
+
        return self._Wy

    @property
    def Wz(self):
        """Regularization matrix Wz"""
-        if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil
+        if getattr(self,'_Wz', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel)
+                self.rz = self.R( f_m , self.eps_q, self.norms[3])
+                self.Rz = Utils.sdiag( self.rz )
+
+            self._Wz = Utils.sdiag((self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
+
        return self._Wz

-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx,)
-            if self.regmesh.dim > 1:
-                wlist += (self.Wy,)
-            if self.regmesh.dim > 2:
-                wlist += (self.Wz,)
-            self._Wsmooth = sp.vstack(wlist)
-        return self._Wsmooth
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
-    @Utils.timeIt
-    def _evalSmall(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m)
-        return 0.5 * r.dot(r)
-
-
-class Sparse(Simple):
-
-    # set default values
-    eps_p      = 1e-1
-    eps_q      = 1e-1
-    curModel = None # use a model to compute the weights
-    gamma    = 1.
-    norms    = [0., 2., 2., 2.]
-    wght     = 1.
-
-    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
-
-        if isinstance(self.wght,float):
-            self.wght = np.ones(self.regmesh.nC) * self.wght
-
-    @property
-    def Wsmall(self):
-        """Regularization matrix Wsmall"""
-        if getattr(self, 'curModel', None) is None:
-            self.Rs = Utils.speye(self.regmesh.nC)
-
-        else:
-            f_m = self.curModel - self.reg.mref
-            self.rs = self.R(f_m , self.eps_p, self.norms[0])
-            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
-            self.Rs = Utils.sdiag( self.rs )
-
-        return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
-
-
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffxStencil * self.curModel
-            self.rx = self.R( f_m , self.eps_q, self.norms[1])
-            self.Rx = Utils.sdiag( self.rx )
-
-        return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
-
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffyStencil * self.curModel
-            self.ry = self.R( f_m , self.eps_q, self.norms[2])
-            self.Ry = Utils.sdiag( self.ry )
-
-        return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
-
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffzStencil * self.curModel
-            self.rz = self.R( f_m , self.eps_q, self.norms[3])
-            self.Rz = Utils.sdiag( self.rz )
-
-        return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
-
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        #if getattr(self, '_Wsmooth', None) is None:
-        wlist = (self.Wx,)
-        if self.regmesh.dim > 1:
-            wlist += (self.Wy,)
-        if self.regmesh.dim > 2:
-            wlist += (self.Wz,)
-        #self._Wsmooth = sp.vstack(wlist)
-        return sp.vstack(wlist)
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
    def R(self, f_m , eps, exponent):

+        # Eta scaling is important for mix-norms...do not mess with it
        eta = (eps**(1.-exponent/2.))**0.5
        r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)

@@ -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):
@@ -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
@@ -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
@@ -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):
@@ -7,3 +7,4 @@ from CounterUtils import *
 import ModelBuilder
 import SolverUtils
 from coordutils import *
+from modelutils import *
@@ -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::
@@ -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()

@@ -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) \
@@ -0,0 +1,63 @@
+from matutils import mkvc, ndgrid
+import numpy as np
+
+def surface2ind_topo(mesh, topo, gridLoc='CC'):
+# def genActiveindfromTopo(mesh, topo):
+    """
+    Get active indices from topography
+    """
+
+
+    if mesh.dim == 3:
+        from scipy.interpolate import NearestNDInterpolator
+        Ftopo = NearestNDInterpolator(topo[:,:2], topo[:,2])
+
+        if gridLoc == 'CC':
+            XY = ndgrid(mesh.vectorCCx, mesh.vectorCCy)
+            Zcc = mesh.gridCC[:,2].reshape((np.prod(mesh.vnC[:2]), mesh.nCz), order='F')
+
+            gridTopo = Ftopo(XY)
+            actind = [gridTopo[ixy] <= Zcc[ixy,:] for ixy in range(np.prod(mesh.vnC[0]))]
+            actind = np.hstack(actind)
+
+        elif gridLoc == 'N':
+
+            XY = ndgrid(mesh.vectorNx, mesh.vectorNy)
+            gridTopo = Ftopo(XY).reshape(mesh.vnN[:2], order='F')
+
+            if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']:
+                raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType)
+
+            Nz = mesh.vectorNz[1:] # TODO: this will only work for tensor meshes
+            actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F')
+
+            for ii in range(mesh.nCx):
+                for jj in range(mesh.nCy):
+                     actind[ii,jj,:] = [np.all(gridTopo[ii:ii+2, jj:jj+2] >= Nz[kk]) for kk in range(len(Nz)) ]
+
+    elif mesh.dim == 2:
+        from scipy.interpolate import interp1d
+        Ftopo = interp1d(topo[:,0], topo[:,1])
+
+        if gridLoc == 'CC':
+            gridTopo = Ftopo(mesh.gridCC[:,0])
+            actind = mesh.gridCC[:,1] <= gridTopo
+
+        elif gridLoc == 'N':
+
+            gridTopo = Ftopo(mesh.vectorNx)
+            if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']:
+                raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType)
+
+            Ny = mesh.vectorNy[1:] # TODO: this will only work for tensor meshes
+            actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F')
+
+            for ii in range(mesh.nCx):
+                actind[ii,:] = [np.all(gridTopo[ii:ii+2] > Ny[kk]) for kk in range(len(Ny)) ]
+
+    else:
+        raise NotImplementedError('surface2ind_topo not implemented for 1D mesh')
+
+    return mkvc(actind)
+
+
@@ -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'
@@ -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
@@ -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 %}
@@ -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
@@ -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:
@@ -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"
@@ -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'),
+]
@@ -7,7 +7,7 @@ Examples
    :maxdepth: 1
    :glob:

-    examples/*
+    ../examples/*


 External Notebooks
@@ -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
@@ -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:
@@ -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:
@@ -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:
@@ -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
 ---------
@@ -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:
@@ -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:
@@ -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:
+
@@ -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:
+


@@ -46,6 +46,8 @@ The API
 =======

 .. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapD
+    :noindex:

 .. autofunction:: SimPEG.Utils.SolverUtils.SolverWrapI
+    :noindex:

@@ -6,5 +6,6 @@ Utilities

   api_Solver
   api_Maps
+   api_PropMaps
   api_Utils
   api_Tests
@@ -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
 -------

@@ -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
@@ -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:
@@ -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
 -----------

@@ -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:
@@ -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:
+
+
@@ -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
@@ -16,6 +16,6 @@ DC Analytic Dipole
    from SimPEG import Examples
    Examples.DC_Analytic_Dipole.run()

-.. literalinclude:: ../../SimPEG/Examples/DC_Analytic_Dipole.py
+.. literalinclude:: ../../../SimPEG/Examples/DC_Analytic_Dipole.py
    :language: python
    :linenos:
@@ -20,8 +20,8 @@ INPUT:
 loc     = Location of spheres [[x1,y1,z1],[x2,y2,z2]]
 radi    = Radius of spheres [r1,r2]
 param   = Conductivity of background and two spheres [m0,m1,m2]
-stype   = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
-dtype   = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
+surveyType = survey type 'pole-dipole' or 'dipole-dipole'
+unitType = Data type "appResistivity" | "appConductivity"  | "volt"
 Created by @fourndo


@@ -31,6 +31,6 @@ Created by @fourndo
    from SimPEG import Examples
    Examples.DC_Forward_PseudoSection.run()

-.. literalinclude:: ../../SimPEG/Examples/DC_Forward_PseudoSection.py
+.. literalinclude:: ../../../SimPEG/Examples/DC_Forward_PseudoSection.py
    :language: python
    :linenos:
@@ -21,6 +21,6 @@ Here we will create and run a FDEM 1D inversion.
    from SimPEG import Examples
    Examples.EM_FDEM_1D_Inversion.run()

-.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py
+.. literalinclude:: ../../../SimPEG/Examples/EM_FDEM_1D_Inversion.py
    :language: python
    :linenos:
@@ -21,6 +21,6 @@ Here we plot the magnetic flux density from a harmonic dipole in a wholespace.
    from SimPEG import Examples
    Examples.EM_FDEM_Analytic_MagDipoleWholespace.run()

-.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py
+.. literalinclude:: ../../../SimPEG/Examples/EM_FDEM_Analytic_MagDipoleWholespace.py
    :language: python
    :linenos:
@@ -17,10 +17,13 @@ current inside a steel-cased. The model is based on the Schenkel and
 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
@@ -53,6 +56,6 @@ citation would be much appreciated!
    from SimPEG import Examples
    Examples.EM_Schenkel_Morrison_Casing.run()

-.. literalinclude:: ../../SimPEG/Examples/EM_Schenkel_Morrison_Casing.py
+.. literalinclude:: ../../../SimPEG/Examples/EM_Schenkel_Morrison_Casing.py
    :language: python
    :linenos:
@@ -21,6 +21,6 @@ Here we will create and run a TDEM 1D inversion.
    from SimPEG import Examples
    Examples.EM_TDEM_1D_Inversion.run()

-.. literalinclude:: ../../SimPEG/Examples/EM_TDEM_1D_Inversion.py
+.. literalinclude:: ../../../SimPEG/Examples/EM_TDEM_1D_Inversion.py
    :language: python
    :linenos:
@@ -47,6 +47,6 @@ Here we reproduce the results from Celia1990_ demonstrating the head-based formu
    from SimPEG import Examples
    Examples.FLOW_Richards_1D_Celia1990.run()

-.. literalinclude:: ../../SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
+.. literalinclude:: ../../../SimPEG/Examples/FLOW_Richards_1D_Celia1990.py
    :language: python
    :linenos:
@@ -21,6 +21,6 @@ Here we go over the basics of creating a linear problem and inversion.
    from SimPEG import Examples
    Examples.Inversion_IRLS.run()

-.. literalinclude:: ../../SimPEG/Examples/Inversion_IRLS.py
+.. literalinclude:: ../../../SimPEG/Examples/Inversion_IRLS.py
    :language: python
    :linenos:
@@ -21,6 +21,6 @@ Here we go over the basics of creating a linear problem and inversion.
    from SimPEG import Examples
    Examples.Inversion_Linear.run()

-.. literalinclude:: ../../SimPEG/Examples/Inversion_Linear.py
+.. literalinclude:: ../../../SimPEG/Examples/Inversion_Linear.py
    :language: python
    :linenos:
@@ -10,7 +10,7 @@


 MT: 1D: Inversion
-=======================
+=================

 Forward model 1D MT data.
 Setup and run a MT 1D inversion.
@@ -22,6 +22,6 @@ Setup and run a MT 1D inversion.
    from SimPEG import Examples
    Examples.MT_1D_ForwardAndInversion.run()

-.. literalinclude:: ../../SimPEG/Examples/MT_1D_ForwardAndInversion.py
+.. literalinclude:: ../../../SimPEG/Examples/MT_1D_ForwardAndInversion.py
    :language: python
    :linenos:
@@ -10,7 +10,7 @@


 MT: 3D: Forward
-=======================
+===============

 Forward model 3D MT data.

@@ -21,6 +21,6 @@ Forward model 3D MT data.
    from SimPEG import Examples
    Examples.MT_3D_Foward.run()

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

 .. --------------------------------- ..
 ..                                   ..
@@ -8,14 +8,20 @@
 ..                                   ..
 .. --------------------------------- ..

-Forward BasicDirectCurrent
-==========================
+
+
+Maps: Mesh2Mesh
+===============
+
+This mapping allows you to go from one mesh to another.
+
+

 .. plot::

    from SimPEG import Examples
-    Examples.Forward_BasicDirectCurrent.run()
+    Examples.Maps_Mesh2Mesh.run()

-.. literalinclude:: ../../SimPEG/Examples/Forward_BasicDirectCurrent.py
+.. literalinclude:: ../../../SimPEG/Examples/Maps_Mesh2Mesh.py
    :language: python
    :linenos:
@@ -0,0 +1,25 @@
+.. _examples_Mesh_Basic_ForwardDC:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+Mesh: Basic Forward 2D DC Resistivity
+=====================================
+
+2D DC forward modeling example with Tensor and Curvilinear Meshes
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.Mesh_Basic_ForwardDC.run()
+
+.. literalinclude:: ../../../SimPEG/Examples/Mesh_Basic_ForwardDC.py
+    :language: python
+    :linenos:
@@ -22,6 +22,6 @@ You can use M.PlotImage to plot images on all of the Meshes.
    from SimPEG import Examples
    Examples.Mesh_Basic_PlotImage.run()

-.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_PlotImage.py
+.. literalinclude:: ../../../SimPEG/Examples/Mesh_Basic_PlotImage.py
    :language: python
    :linenos:
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Lindsey Heagy	a3c0c10e48	cleanup in mapping tests	2016-07-20 09:24:32 -07:00
Lindsey Heagy	2c429f9dc2	Parametric maps renamed to ParametricXXXMap. Some Pep8-ing	2016-07-18 10:31:25 -07:00
Lindsey Heagy	62b3a3a218	light pep8 formatting in tests	2016-07-16 16:26:36 -07:00
Lindsey Heagy	6826df8989	test all new maps unless added to exclude list	2016-07-16 16:26:14 -07:00
Lindsey Heagy	78fce342cf	typo fix	2016-07-16 05:42:06 -07:00
Lindsey Heagy	94f799e7c4	fix assert	2016-07-16 05:38:04 -07:00
Lindsey Heagy	947c671c8c	Merge branch 'dev' into feat/mappingderivs	2016-07-15 14:29:46 -07:00
Lindsey Heagy	0367a172a2	light docs and pep8 cleanup	2016-07-15 13:50:02 -07:00
SEOGI KANG	4796b0f91f	Merge pull request #357 from simpeg/analytics Analytics	2016-06-30 00:25:28 -07:00
seogi_macbook	52b25e2dc5	Merge branch 'dev' of https://github.com/simpeg/simpeg into analytics	2016-06-30 00:23:12 -07:00
seogi_macbook	a289b656cd	Fixes for kwargs variables in FDEMDipolarfields.py	2016-06-29 13:09:11 -07:00
Lindsey Heagy	334cd8e454	Bump version: 0.1.11 → 0.1.12	2016-06-29 09:49:29 -07:00
Lindsey Heagy	ecbdd90f63	Merge pull request #354 from simpeg/dev Two new examples.	2016-06-29 09:46:33 -07:00
Lindsey Heagy	dcdcbd212a	Update Maps.py	2016-06-29 08:49:56 -07:00
Lindsey Heagy	17b1459b57	Update Maps.py	2016-06-29 08:48:35 -07:00
dfournier	394dc9106a	Merge pull request #332 from simpeg/ref/regularization Automate the epsilon picking based on percentile of model values for …	2016-06-29 08:39:29 -07:00
seogi_macbook	eda2394411	fix bug for omega.	2016-06-27 13:04:30 -07:00
Rowan Cockett	3deca9ed77	Merge pull request #351 from simpeg/example/mesh2mesh Mesh2Mesh and Combo Map examples.	2016-06-26 21:22:28 -06:00
Rowan Cockett	ba173674ec	Mesh2Mesh and Combo Map examples. Also fixed plotting codes to show the plots by default.	2016-06-26 17:07:07 -06:00
Rowan Cockett	303da372aa	Merged branch master into dev	2016-06-26 16:31:24 -06:00
Rowan Cockett	6d6e7fc8bd	Merge pull request #350 from simpeg/fix/docs-images Update index.rst	2016-06-26 16:30:29 -06:00
Rowan Cockett	8ed3ec18fa	Update README.rst	2016-06-26 16:29:20 -06:00
Rowan Cockett	3960cfc313	Update index.rst	2016-06-26 16:17:30 -06:00
Rowan Cockett	f4a8efab78	Whitespace.	2016-06-26 16:12:04 -06:00
Rowan Cockett	2eba0b841f	Merge pull request #338 from simpeg/dev Dev	2016-06-26 14:01:03 -06:00
sgkang	c79bb998cb	Merge pull request #348 from simpeg/analytics fix minor bugs in analytics (just for binder deploy)	2016-06-23 14:14:29 -07:00
seogi_macbook	0763925743	fix minor bugs in analytics	2016-06-23 14:12:49 -07:00
sgkang	1a0b81a206	Merge pull request #341 from simpeg/analytics Analytics	2016-06-23 11:33:33 -07:00
seogi_macbook	8b44f8d96b	change FDEM_fields.py to FDEMDipolarfield.py	2016-06-23 10:14:13 -07:00
seogi_macbook	2bfd01ed7c	Merge branch 'dev' of https://github.com/simpeg/simpeg into analytics	2016-06-23 10:09:56 -07:00
seogi_macbook	e1ba80883d	Incorporate Lindsey's suggestoins	2016-06-23 09:10:50 -07:00
Rowan Cockett	a54713f546	Change to NotImplementedError.	2016-06-22 11:32:54 -06:00
Lindsey Heagy	2f4c9a2a7a	testing for maps times a vec	2016-06-21 18:49:41 -07:00
Lindsey Heagy	b64d967e73	mapping derivs can take a vector and return a vector as per #342	2016-06-21 18:49:22 -07:00
Lindsey Heagy	ef382aed85	Merge branch 'master' into dev	2016-06-21 18:43:28 -06:00
seogi_macbook	1b33804e5a	Merge branch 'dev' of https://github.com/simpeg/simpeg into analytics	2016-06-21 11:15:23 -07:00
seogi_macbook	c161c5eab2	Merge branch 'master' of https://github.com/simpeg/simpeg into analytics	2016-06-21 11:14:56 -07:00
micmitch	f0944362c8	Silly mistake... needed zero arrays instead of scalars.	2016-06-14 15:22:36 -07:00
micmitch	e815ddaec7	Removed d typos from the end of function names.	2016-06-14 15:15:53 -07:00
Lindsey	64b0b4561f	Merge pull request #337 from simpeg/feat/docs-gae-travis-deploy Feat/docs gae travis deploy	2016-06-13 19:16:22 -06:00
micmitch	b9d30af4a8	Changed \sigma to \hat{\sigma} = \sigma + i \omega \epsilon in the E field calculations.	2016-06-13 17:47:46 -07:00
micmitch	093f441331	Added functions to split electric field into "galvanic" and "inductive" portions.	2016-06-13 17:35:52 -07:00
seogi_macbook	5e3c1da8e2	add place holder for galvanic and inductive electric fields...	2016-06-13 16:41:59 -07:00
seogi_macbook	a3a5c86008	Fix couple bugs in FDEM analytics	2016-06-11 21:58:03 +02:00
Lindsey Heagy	7a82f57367	- doc the simple regularization, sparse regularization and regularization mesh - typo fix in MT_3D_Forward example - Solver - use EM.Static.DC for DC example	2016-06-11 08:16:22 -07:00
Lindsey Heagy	845c3c10ab	version on travis commit	2016-06-11 07:48:24 -07:00
Lindsey Heagy	c4d86b4a29	Merge branch 'dev' into feat/docs-gae-travis-deploy	2016-06-11 07:14:00 -07:00
Lindsey Heagy	3fc855f3c9	only push on master	2016-06-10 20:50:41 -07:00
Lindsey Heagy	fea508a507	fix path to app.yaml	2016-06-10 19:39:02 -07:00
Lindsey Heagy	5ca52cf49f	linting travis deploy of docs	2016-06-10 19:14:24 -07:00
Lindsey Heagy	7a06453c42	editing travis deploy	2016-06-10 18:54:31 -07:00
Lindsey Heagy	291da78b97	try deploying only from docs branch	2016-06-10 15:16:56 -07:00
D Fournier	ef12a3674a	Automate the epsilon picking based on percentile of model values for DEFAULT mode. Fix example. Fix bug with Maps using array of values	2016-06-06 12:28:01 -07:00
Lindsey Heagy	504592c8de	travis typo fixes	2016-06-01 08:12:21 -07:00
Lindsey Heagy	75647f8fc3	first pass at gae deploy	2016-06-01 00:00:27 -07:00
Lindsey Heagy	e3462666bd	working on travis decrypt	2016-05-31 23:27:32 -07:00
Lindsey Heagy	fb37bf0fe2	move unpacking of credentials to after success, point credentials to docs folder	2016-05-31 23:15:50 -07:00
Lindsey Heagy	b023adbb33	add encrypted credentials for gae site and decrypt on travis	2016-05-31 23:03:55 -07:00
Lindsey Heagy	279bd49b4c	Bump version: 0.1.10 → 0.1.11	2016-05-31 15:38:05 -07:00
Lindsey Heagy	3398d5ab4d	use templates in conf.py, add google analytics to docs	2016-05-31 15:33:19 -07:00
Lindsey	3e26bb7de9	Merge pull request #277 from simpeg/dev Dev	2016-05-31 15:11:29 -07:00
Lindsey	382b31bd12	Merge pull request #310 from simpeg/ref/regularization modularizing regularization	2016-05-31 14:09:54 -07:00
Lindsey Heagy	01e19e4227	start of gae site	2016-05-30 22:01:15 -07:00
Lindsey Heagy	9fbdaaf0a5	fix typo in examples path	2016-05-30 21:12:31 -07:00
Lindsey Heagy	2a159e20b9	update docs for examples to point to the correct path	2016-05-30 20:34:10 -07:00
Lindsey Heagy	5e8d3fbc78	organizing the docs - put the content in a `content` folder. put the SimPEG core api docs in `core_api`	2016-05-30 17:06:29 -07:00
Lindsey Heagy	414418a996	Merge branch 'dev' into feat/docs-deploy # Conflicts: # SimPEG/Mesh/View.py	2016-05-30 15:56:24 -07:00
Lindsey	6cc509020a	Merge pull request #323 from simpeg/feat/plotImage-curvilinear plotImage for curvilinear mesh	2016-05-30 08:27:56 -07:00
Rowan Cockett	40f0874dfb	Doc testing, I think that is most of them!	2016-05-29 22:18:22 -07:00
Rowan Cockett	231e6dbc93	Suppress image warning.	2016-05-29 19:22:15 -07:00
Rowan Cockett	18f98b2ecd	Merge branch 'docs' of https://github.com/simpeg/simpeg into feat/docs-deploy # Conflicts: # SimPEG/Utils/meshutils.py # docs/api_Utils.rst # docs/conf.py # docs/flow/index.rst	2016-05-29 19:18:36 -07:00
Rowan Cockett	bc073e49b5	Updates to docs errors.	2016-05-29 18:57:38 -07:00
Rowan Cockett	a131383dae	Correct solver location.	2016-05-29 18:35:41 -07:00
Rowan Cockett	ad4a0240d1	Remove Vertical1DMap from tests.	2016-05-29 18:34:17 -07:00
Rowan Cockett	74f5395573	Surject1D updates.	2016-05-29 18:25:42 -07:00
Rowan Cockett	6f7a0b1279	Rename `Vertical1DMap` to `SurjectVertical1D` due to depreciation.	2016-05-29 18:14:42 -07:00
Rowan Cockett	e2bb9c8d8e	rename flow example.	2016-05-29 18:12:52 -07:00
Rowan Cockett	de693adaa7	Minor updates to get it 'working' There still seems to be a problem with this example: - The line search breaks. - The plots are not informative. - There are a lot of errors in the structured array codes.	2016-05-29 18:04:09 -07:00
Rowan Cockett	fc07993006	Spacings in functions.	2016-05-29 18:03:03 -07:00
Rowan Cockett	12a12c7b5a	updates to FDEM docs.	2016-05-29 17:51:41 -07:00
Rowan Cockett	0b4215f33e	Add DC and IP docs.	2016-05-29 17:45:49 -07:00
Rowan Cockett	5e2a8232a3	Minor updates to titles in examples.	2016-05-29 17:21:37 -07:00
Rowan Cockett	feba384911	Add solver parameter to the Casing example.	2016-05-29 17:08:11 -07:00
Rowan Cockett	4844b7230a	TOC updates for docs index.	2016-05-29 17:05:42 -07:00
Lindsey Heagy	f2e13182bf	Merge branch 'dev' into feat/plotImage-curvilinear # Conflicts: # SimPEG/Examples/__init__.py	2016-05-29 16:43:40 -07:00
D Fournier	09cd9c7fa3	Merge branch 'dev' into ref/regularization	2016-05-29 16:26:58 -07:00
Lindsey Heagy	62eb4541cb	update example name --> based on mesh	2016-05-29 15:48:15 -07:00
Lindsey	09eb2106ec	Merge pull request #307 from simpeg/ref/dev Ref/dev	2016-05-29 14:50:43 -07:00
D Fournier	f8b86abd5a	Merge branch 'ref/dev' into ref/regularization	2016-05-29 14:10:46 -07:00
D Fournier	e476bf0059	Cleanup Sparse Reg and Directives	2016-05-29 14:09:29 -07:00
Rowan Cockett	825511e9d3	Add a curvilinear plotImage function, update example.	2016-05-29 13:17:22 -07:00
D Fournier	3b4bec9c0b	Refactor IRLS iterations, full solves from l2->lp Adapt Example	2016-05-28 11:27:09 -07:00
micmitch	1960b52dfd	First stab at analytic functions for the fields from a harmonic electric dipole source. Not sure about the exception that I try to throw if multiple frequencies and multiple evaluation locations are both specified.	2016-05-27 14:56:48 -07:00
D Fournier	022e1f7660	Update IRLS directive to allow multiple GN iterations. Remove modifications to the ProjGN solver. Update IRLS example.	2016-05-27 13:11:31 -07:00
seogi_macbook	28d67e3112	Start of ED !!	2016-05-27 11:26:47 -07:00
D Fournier	406703f1c6	Merge branch 'dev' into ref/dev Conflicts: docs/examples/DC_Forward_PseudoSection.rst	2016-05-27 11:10:39 -07:00
D Fournier	7b72d3a92d	Merge branch 'dev' into ref/dev Conflicts: SimPEG/DCIP/DCIPUtils.py	2016-05-27 10:01:28 -07:00
Lindsey	cf89f5f6a2	Merge pull request #322 from simpeg/bug/propmap Bug/propmap	2016-05-26 20:52:42 -07:00
Lindsey Heagy	aa1086eba3	use fixed prop map in EM	2016-05-26 18:03:09 -07:00
Lindsey Heagy	1c53129da6	fix bug in prop map linked derivs	2016-05-26 17:58:30 -07:00
sgkang	6fd3be77de	Merge pull request #304 from simpeg/dcip/dev Dcip/dev	2016-05-26 13:27:44 -07:00
D Fournier	fd3bde787f	Propose change to the Projected_GNCG solver. Add inner GN iterations. Nice improvement to the convergence of IRLS	2016-05-12 14:58:16 -07:00
D Fournier	3cc46131a3	Temporary change ... comment out W and Wsmooth	2016-05-12 08:31:01 -07:00
D Fournier	cd2360b815	Stash the regularization between each beta	2016-05-11 23:04:14 -07:00
D Fournier	e10d6878fb	Remove Wsmooth from def W and replace by parts	2016-05-11 07:58:06 -07:00
Lindsey Heagy	3dd9ecc9cd	fix tikhonov 2Deriv	2016-05-10 22:10:19 -07:00
Lindsey Heagy	90a3030796	fixed 2 deriv	2016-05-10 21:39:15 -07:00
D Fournier	7964ebce50	Update directive to None the Wsmooth after iteration.	2016-05-10 17:20:46 -07:00
Lindsey Heagy	955bd54019	notation cleanup in Regularization	2016-05-10 16:46:11 -07:00
Lindsey Heagy	2a802c1aa3	weights --> cell_weights, removed vol term from simple regularization	2016-05-10 16:39:47 -07:00
Lindsey Heagy	3f0c89f10b	remove extra Ws	2016-05-10 14:53:43 -07:00
Lindsey Heagy	eaa37f42e4	remove duplicate evalSmall	2016-05-10 14:51:56 -07:00
D Fournier	fb5434695f	Alpha_s default to 1.0	2016-05-10 14:31:45 -07:00
D Fournier	e037597ecd	Merge branch 'feat/sparse-regularization' into ref/regularization	2016-05-10 13:37:16 -07:00
D Fournier	b4ab60c260	Add model mapping to sparse regularization	2016-05-05 11:55:56 -07:00
Lindsey Heagy	fbb8cf2731	modularizing regularization	2016-05-04 23:17:01 -07:00
Lindsey Heagy	0379df2bf2	attempt to clean up docs in DCIP utils	2016-05-04 22:27:02 -07:00
Lindsey Heagy	66440b0478	add depreciation warnings to DCIP utils for activeind from topo	2016-05-04 22:14:41 -07:00
Lindsey Heagy	dbdcc3cefb	use sigma in MfRhoDeriv - due to propmap bug	2016-05-04 22:06:32 -07:00
Lindsey Heagy	dd45a6a085	name updates in DC_Forward_PseudoSection, DC_Utils, example for Utils_surface2ind_topo	2016-05-02 11:40:02 -07:00
Lindsey Heagy	ba8f270b3a	start of surface2ind_topo	2016-05-01 13:17:16 -07:00
Lindsey Heagy	4257ea77b3	remove InjectActiveCellsTopo. you should use InjectActiveCells	2016-04-29 15:09:04 -07:00
Lindsey Heagy	a0174e4f30	kwarg name updates	2016-04-29 12:52:45 -07:00
Lindsey Heagy	756b738ef2	Merge branch 'dev' into docs	2016-04-05 17:47:17 -07:00
Lindsey Heagy	021e7c794c	Merge branch 'master' into docs # Conflicts: # SimPEG/Mesh/TensorMesh.py	2016-03-29 22:53:16 -07:00
Lindsey Heagy	b5b70390cb	Merge branch 'dev' into docs	2016-03-06 21:50:27 -08:00
Lindsey Heagy	6acaa81faf	fixed merge conflicts in FDEM docs that I missed	2016-02-09 09:03:20 -08:00
Lindsey Heagy	312b5d79c5	resolved merge conflicts in TensorMesh	2016-02-09 08:53:01 -08:00
Lindsey Heagy	999a37547e	Merge branch 'dev' into docs # Conflicts: # .travis.yml # SimPEG/EM/FDEM/FDEM.py # SimPEG/Mesh/TensorMesh.py	2016-02-09 08:32:41 -08:00
Lindsey Heagy	cbe8758465	corrected scipy.sparse.csr_matrix, move size descriptions to :return: instead of :type:	2016-02-01 08:22:00 -08:00
Lindsey Heagy	6b359f49b5	docs clean-up (using autoclass is more stable than automodule)	2016-01-31 15:21:46 -08:00
Lindsey Heagy	2254eedbac	indentations clean up in FDEM.py	2016-01-31 13:54:39 -08:00
Lindsey Heagy	012d2cadf1	use intersphinx mapping to get numpy, scipy, matplotlib, python inventories	2016-01-31 13:54:24 -08:00
Lindsey Heagy	841ba61006	clean up the html build	2016-01-31 13:22:36 -08:00
Lindsey Heagy	2874e204ee	exclude _static from warnings	2016-01-31 12:52:42 -08:00
Lindsey Heagy	adca273565	ignore nonlocal images in sphinx build	2016-01-31 12:46:24 -08:00
Lindsey Heagy	d9d6f70958	better description of paths in test_docs	2016-01-31 12:45:51 -08:00
Lindsey Heagy	f4ef767764	seperate out docs test so it runs independently (not on every test)	2016-01-31 12:08:29 -08:00
Lindsey Heagy	e314bdb740	add sphinx to travis conda install	2016-01-31 11:06:47 -08:00
Lindsey Heagy	e005ed8f5f	use cd to get into docs directories for testing	2016-01-31 10:51:55 -08:00
Lindsey Heagy	ac2e38e89d	test docs first	2016-01-31 10:47:00 -08:00
Lindsey Heagy	ade37fb493	add travis to docs. nit-picky testing on html, latex, link check	2016-01-31 09:51:57 -08:00