Merge branch 'dev' into em/dev

2026-06-27 19:48:52 +08:00 · 2016-05-09 11:32:59 -07:00
parent 6e00b4c2fe 906cca30f3
commit abd919e862
29 changed files with 1912 additions and 596 deletions
@@ -25,6 +25,10 @@ SimPEG
    :target: https://coveralls.io/r/simpeg/simpeg?branch=master
    :alt: Coverage status

+.. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
+    :alt: gitter chat room at https://gitter.im/simpeg/simpeg
+    :target: https://gitter.im/simpeg/simpeg
+    
 Simulation and Parameter Estimation in Geophysics  -  A python package for simulation and gradient based parameter estimation in the context of geophysical applications.

 The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:
@@ -169,7 +169,7 @@ def readUBC_DC2DModel(fileName):

    return model

-def plot_pseudoSection(DCsurvey, axs, stype):
+def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", 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.
@@ -179,7 +179,7 @@ def plot_pseudoSection(DCsurvey, axs, stype):
        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

@@ -192,9 +192,6 @@ def plot_pseudoSection(DCsurvey, axs, stype):
    from scipy.interpolate import griddata
    import pylab as plt

-    # Set depth to 0 for now
-    z0 = 0.
-
    # Pre-allocate
    midx = []
    midz = []
@@ -221,47 +218,97 @@ def plot_pseudoSection(DCsurvey, axs, stype):
        Cmid = (Tx[0][0] + Tx[1][0])/2
        Pmid = (Rx[0][:,0] + Rx[1][:,0])/2

-        # Compute pant leg of apparent rho
-        if stype == 'pdp':
-            leg =  data * 2*np.pi  * MA * ( MA + MN ) / MN
+        # Change output for dtype
+        if dtype == 'volt':

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

-        elif stype == 'dpdp':
-            leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
+        else:

+            # Compute pant leg of apparent rho
+            if stype == 'pdp':
+
+                leg =  data * 2*np.pi  * MA * ( MA + MN ) / MN
+
+            elif stype == 'dpdp':
+
+                leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
+
+            else:
+                print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
+                break
+
+
+            if dtype == 'appc':
+
+                leg = np.log10(abs(1./leg))
+                rho = np.hstack([rho,leg])
+
+            elif dtype == 'appr':
+
+                leg = np.log10(abs(leg))
+                rho = np.hstack([rho,leg])
+
+            else:
+                print """dtype must be 'appr' | 'appc' | 'volt' """
+                break

        midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
-        midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + z0 ])
-        rho = np.hstack([rho,leg])
-
-
-    ax = axs
+        midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ])

    # 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()
+    else:
+        vmin, vmax = clim[0], clim[1]
+
+    # 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':
+            cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal")
+
+        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)


-    plt.imshow(grid_rho.T, extent = (np.min(midx),np.max(midx),np.min(midz),np.max(midz)), origin='lower', alpha=0.8, vmin = np.min(rho), vmax = np.max(rho))
-    cbar = plt.colorbar(format = '%.2f',fraction=0.04,orientation="horizontal")

-    cmin,cmax = cbar.get_clim()
-    ticks = np.linspace(cmin,cmax,3)
-    cbar.set_ticks(ticks)
+    if not axlabel:
+        axs.set_xticklabels([])
+        axs.set_yticklabels([])

-    # Plot apparent resistivity
-    plt.scatter(midx,midz,s=50,c=rho.T)
-
-    ax.set_xticklabels([])
-
-    ax.set_ylabel('Z')
-    ax.yaxis.tick_right()
-    ax.yaxis.set_label_position('right')
    plt.gca().set_aspect('equal', adjustable='box')


-    return ax
+
+    return ph

 def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
    """
@@ -361,16 +408,6 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
                srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:])
            SrcList.append(srcClass)

-#==============================================================================
-#     elif re.match(stype,'dpdp'):
-#
-#         for ii in range(0, int(nstn)-2):
-#
-#             indx = np.min([ii+n+1,nstn])
-#             Tx.append(np.c_[M[ii,:],N[ii,:]])
-#             Rx.append(np.c_[M[ii+2:indx,:],N[ii+2:indx,:]])
-#==============================================================================
-
    elif stype == 'gradient':

        # Gradient survey only requires Tx at end of line and creates a square
@@ -423,15 +460,15 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
    survey = DC.SurveyDC(SrcList)
    return survey, Tx, Rx

-def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
+def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', 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'
+        :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
@@ -446,10 +483,16 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):

    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'"
-
+    
    fid = open(fileName,'w')
-    fid.write('! ' + stype + ' 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):
@@ -473,7 +516,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
                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')


@@ -481,18 +524,25 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):

                if stype == 'SURFACE':

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

                if stype == '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]

+                    # 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='%e',delimiter=' ',newline='\n')
+                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':

@@ -504,31 +554,32 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):

            if stype == 'GENERAL':

-                fid.writelines("%e " % ii for ii in mkvc(tx))
+                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):
+def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
    """
-        Read DC survey and data and change
-        coordinate system to distance along line assuming
-        all data is acquired along line.
-        First transmitter pole is assumed to be at the origin
+        Read DC survey and projects the coordinate system
+        according to the flag = 'Xloc' | 'Yloc' | 'local' (default)
+        In the 'local' system, station coordinates are referenced
+        to distance from the first srcLoc[0].loc[0]

-        Assumes flat topo for now...
+        The Z value is preserved, but Y coordinates zeroed.

        Input:
-        :param Tx, Rx
+        :param survey3D

        Output:
-        :figure Tx2d, Rx2d
+        :figure survey2D

-        Edited Feb 17th, 2016
+        Edited April 6th, 2016

        @author: dominiquef

@@ -570,25 +621,39 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID):
            Rx = DCsurvey.srcList[indx[ii]].rxList[0].locs
            nrx = Rx[0].shape[0]

-            # Find A electrode along line
-            vec, r = r_unit(x0,Tx[ii][0,0:2])
-            A = stn_id(vecTx,vec,r)
+            if flag == 'local':
+                # Find A electrode along line
+                vec, r = r_unit(x0,Tx[ii][0,0:2])
+                A = stn_id(vecTx,vec,r)

-            # Find B electrode along line
-            vec, r = r_unit(x0,Tx[ii][1,0:2])
-            B = stn_id(vecTx,vec,r)
+                # Find B electrode along line
+                vec, r = r_unit(x0,Tx[ii][1,0:2])
+                B = stn_id(vecTx,vec,r)

-            M = np.zeros(nrx)
-            N = np.zeros(nrx)
-            for kk in range(nrx):
+                M = np.zeros(nrx)
+                N = np.zeros(nrx)
+                for kk in range(nrx):

-                # Find all M electrodes along line
-                vec, r = r_unit(x0,Rx[0][kk,0:2])
-                M[kk] = stn_id(vecTx,vec,r)
+                    # Find all M electrodes along line
+                    vec, r = r_unit(x0,Rx[0][kk,0:2])
+                    M[kk] = stn_id(vecTx,vec,r)

-                # Find all N electrodes along line
-                vec, r = r_unit(x0,Rx[1][kk,0:2])
-                N[kk] = stn_id(vecTx,vec,r)
+                    # Find all N electrodes along line
+                    vec, r = r_unit(x0,Rx[1][kk,0:2])
+                    N[kk] = stn_id(vecTx,vec,r)
+            elif flag == 'Yloc':
+                """ Flip the XY axis locs"""
+                A = Tx[ii][0,1]
+                B = Tx[ii][1,1]
+                M = Rx[0][:,1]
+                N = Rx[1][:,1]
+
+            elif flag == 'Xloc':
+                """ Copy the rx-tx locs"""
+                A = Tx[ii][0,0]
+                B = Tx[ii][1,0]
+                M = Rx[0][:,0]
+                N = Rx[1][:,0]

            Rx = DC.RxDipole(np.c_[M,np.zeros(nrx),Rx[0][:,2]],np.c_[N,np.zeros(nrx),Rx[1][:,2]])

@@ -601,51 +666,59 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID):
    DCsurvey2D.std = np.asarray(DCsurvey.std)

    return DCsurvey2D
-
-def readUBC_DC3Dobs(fileName):
+    
+def readUBC_DC3Dobs(fileName, dtype = 'DC'):
    """
-        Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location
+        Read UBC GIF IP 3D observation file and generate survey

        Input:
        :param fileName, path to the UBC GIF 3D obs file

        Output:
-        :param rx, tx, d, wd
+        :param IPsurvey
        :return

-        Created on Mon December 7th, 2015
-
        @author: dominiquef

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

    # Countdown for number of obs/tx
    count = 0
    for ii in range(obsfile.shape[0]):

+        # Skip if blank line
        if not obsfile[ii]:
            continue

-        # First line is transmitter with number of receivers
+        # First line or end of a transmitter block, read transmitter info
        if count==0:
-
-            temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
+            # Read the line
+            temp = (np.fromstring(obsfile[ii], dtype=float, sep=' ').T)
            count = int(temp[-1])

            # Check if z value is provided, if False -> nan
            if len(temp)==5:
-                tx = np.r_[temp[0:2],np.nan,temp[0:2],np.nan]
-                zflag = False
+                tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan]
+                
+                zflag = False # Pass on the flag to the receiver loc

            else:
                tx = temp[:-1]
@@ -653,8 +726,16 @@ def readUBC_DC3Dobs(fileName):
            rx = []
            continue

-        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
+        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string

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

            rx.append(temp[:-2])
@@ -664,7 +745,7 @@ def readUBC_DC3Dobs(fileName):
                wd.append(temp[-1])

        else:
-            rx.append(np.r_[temp[0:2],np.nan,temp[0:2],np.nan] )
+            rx.append(np.r_[temp[0:2],np.nan,temp[2:4],np.nan] )
            # Check if there is data with the location
            if len(temp)==6:
                d.append(temp[-2])
@@ -672,7 +753,7 @@ def readUBC_DC3Dobs(fileName):

        count = count -1

-        # Reach the end of transmitter block
+        # Reach the end of transmitter block, append the src, rx and continue
        if count == 0:
            rx = np.asarray(rx)
            Rx = DC.RxDipole(rx[:,:3],rx[:,3:])
@@ -688,6 +769,7 @@ def readUBC_DC3Dobs(fileName):

 def readUBC_DC2Dobs(fileName):
    """
+        ------- NEEDS TO BE UPDATED ------
        Read UBC GIF 2D observation file and generate arrays for tx-rx location

        Input:
@@ -735,6 +817,73 @@ def readUBC_DC2Dobs(fileName):

    return tx, rx, d, wd

+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
+
+        Created on Mon March 9th, 2016 << Doug's 70th Birthday !! >>
+
+        @author: dominiquef
+
+    """
+
+    # Load file
+    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
+
+    # Pre-allocate
+    srcLists = []
+    Rx = []
+    d = []
+    zflag = True # Flag for z value provided
+
+    for ii in range(obsfile.shape[0]):
+
+        if not obsfile[ii]:
+            continue
+
+        # First line is transmitter with number of receivers
+
+
+        temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
+
+
+        # Check if z value is provided, if False -> nan
+        if len(temp)==5:
+            tx = np.r_[temp[0],np.nan,np.nan,temp[1],np.nan,np.nan]
+            zflag = False
+
+        else:
+            tx = np.r_[temp[0],np.nan,temp[1],temp[2],np.nan,temp[3]]
+
+
+        if zflag:
+            rx = np.c_[temp[4],np.nan,temp[5],temp[6],np.nan,temp[7]]
+
+
+        else:
+            rx = np.c_[temp[2],np.nan,np.nan,temp[3],np.nan,np.nan]
+            # Check if there is data with the location
+
+        d.append(temp[-1])
+
+
+        Rx = DC.RxDipole(rx[:,:3],rx[:,3:])
+        srcLists.append( DC.SrcDipole( [Rx], tx[:3],tx[3:]) )
+
+    # Create survey class
+    survey = DC.SurveyDC(srcLists)
+
+    survey.dobs = np.asarray(d)
+
+    return {'DCsurvey':survey}
+
 def readUBC_DC2DMesh(fileName):
    """
        Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg
@@ -928,7 +1077,6 @@ def getSrc_locs(DCsurvey):

    srcMat = np.zeros((DCsurvey.nSrc,2,3))
    for ii in range(DCsurvey.nSrc):
-        print np.asarray(DCsurvey.srcList[ii].loc).shape
        srcMat[ii,:,:] =  np.asarray(DCsurvey.srcList[ii].loc)

    return srcMat
@@ -216,13 +216,13 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
        # Save the data.
        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
        phi_ms = 0.5*ms.dot(ms)
-        if self.reg.smoothModel == True:
+        if self.reg.mrefInSmooth == True:
            mref = self.reg.mref
        else:
            mref = 0
        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
        phi_mx = 0.5 * mx.dot(mx)
-        if self.prob.mesh.dim==2:
+        if self.prob.mesh.dim >= 2:
            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
            phi_my = 0.5 * my.dot(my)
        else:
@@ -237,40 +237,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 SaveOutputDictEveryIteration(_SaveEveryIteration):
-    """SaveOutputDictEveryIteration
-    A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
-    """
-
-    def initialize(self):
-        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName
-
-    def endIter(self):
-        # Save the data.
-        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
-        phi_ms = 0.5*ms.dot(ms)
-        if self.reg.smoothModel == True:
-            mref = self.reg.mref
-        else:
-            mref = 0
-        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
-        phi_mx = 0.5 * mx.dot(mx)
-        if self.prob.mesh.dim==2:
-            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
-            phi_my = 0.5 * my.dot(my)
-        else:
-            phi_my = 'NaN'
-        if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType:
-            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
-            phi_mz = 0.5 * mz.dot(mz)
-        else:
-            phi_mz = 'NaN'
-
-
-        # Save the file as a npz
-        np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), 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):

@@ -283,3 +249,125 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 #             mref = self.mref0
 #         self.m_prev = self.invProb.m_current
 #         return mref
+
+class Update_IRLS(InversionDirective):
+
+    eps_min = None
+    factor = None
+    gamma = None
+    phi_m_last = None
+    phi_d_last = None
+
+
+    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
+
+    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])
+            else:
+                self.reg.eps = eps
+
+        # Get phi_m at the end of current iteration
+        self.phi_m_last = self.invProb.phi_m_last
+
+         # Update the model used for the IRLS weights
+        self.reg.curModel = self.invProb.curModel
+
+        # Temporarely set gamma to 1. to get raw phi_m
+        self.reg.gamma = 1.
+
+        # 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
+
+        # Set the weighting matrix to None so that it is recomputed next time
+        # it is called in the inversion
+        self.reg._W = None
+
+class Update_lin_PreCond(InversionDirective):
+    """
+    Create a Jacobi preconditioner for the linear problem
+    """
+    onlyOnStart=False
+    
+    def initialize(self):
+    
+        if getattr(self.opt, 'approxHinv', None) is None:
+            # Update the pre-conditioner
+            diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
+            self.opt.approxHinv = PC
+            
+    def endIter(self):
+        # Cool the threshold parameter
+        if self.onlyOnStart==True:
+            return
+            
+        if getattr(self.opt, 'approxHinv', None) is not None:
+            # Update the pre-conditioner
+            diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
+            self.opt.approxHinv = PC
+
+
+class Update_Wj(InversionDirective):
+    """
+        Create approx-sensitivity base weighting using the probing method
+    """
+    k = None # Number of probing cycles
+    itr = None # Iteration number to update Wj, or always update if None
+
+    def endIter(self):
+
+        if self.itr is None or self.itr == self.opt.iter:
+
+            m = self.invProb.curModel
+            if self.k is None:
+                self.k = int(self.survey.nD/10)
+
+            def JtJv(v):
+
+                Jv = self.prob.Jvec(m, v)
+
+                return self.prob.Jtvec(m,Jv)
+
+            JtJdiag = Utils.diagEst(JtJv,len(m),k=self.k)
+            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
@@ -2,19 +2,27 @@ 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', plotIt=True):
+def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', plotIt=True):
    """
        DC Forward Simulation
        =====================

-        Forward model conductive spheres in a half-space and plot a pseudo-section
+        Forward model two conductive spheres in a half-space and plot a
+        pseudo-section. Assumes an infinite line source and measures along the
+        center of the spheres.

-        Created by @fourndo on Mon Feb 01 19:28:06 2016
+        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)
+        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)"

    if loc is None:
        loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]]
@@ -27,7 +35,6 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):


    # First we need to create a mesh and a model.
-
    # This is our mesh
    dx    = 5.

@@ -52,14 +59,10 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
    # Get index of the center
    indy = int(mesh.nCy/2)

-
    # Plot the model for reference
    # Define core mesh extent
    xlim = 200
-    zlim = 125
-
-    # Specify the survey type: "pdp" | "dpdp"
-
+    zlim = 100

    # Then specify the end points of the survey. Let's keep it simple for now and survey above the anomalies, top of the mesh
    ends = [(-175,0),(175,0)]
@@ -77,12 +80,13 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
    dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )
    dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len
    dl_y = ( Tx[-1][1,1] - Tx[0][1,0]  ) / dl_len
-    azm =  np.arctan(dl_y/dl_x)
+    #azm =  np.arctan(dl_y/dl_x)

    #Set boundary conditions
    mesh.setCellGradBC('neumann')

-    # Define the differential operators needed for the DC problem
+    # Define the linear system needed for the DC problem. We assume an infitite
+    # line source for simplicity.
    Div = mesh.faceDiv
    Grad = mesh.cellGrad
    Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
@@ -145,16 +149,23 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
    print 'Forward completed'

    # Let's just convert the 3D format into 2D (distance along line) and plot
-    # [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc))
-    survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc))
+    survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc) , 'Xloc')
    survey2D.dobs =np.hstack(data)
-    # Here is an example for the first tx-rx array
+
    if plotIt:
        import matplotlib.pyplot as plt
-        fig = plt.figure()
+        fig = plt.figure(figsize=(7,7))
        ax = plt.subplot(2,1,1, aspect='equal')
-        mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True)
-        ax.set_title('E-W section at '+str(mesh.vectorCCy[indy])+' m')
+        # 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)
+        ax.add_artist(circle1)
+        ax.add_artist(circle2)
+
+        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')
        plt.gca().set_aspect('equal', adjustable='box')

        plt.scatter(Tx[0][0,:],Tx[0][2,:],s=40,c='g', marker='v')
@@ -163,22 +174,34 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])


-        ax = plt.subplot(2,1,2, aspect='equal')
+        pos =  ax.get_position()
+        ax.set_position([pos.x0 , pos.y0 + 0.025 ,  pos.width, pos.height])
+        pos =  ax.get_position()
+        cbarax = fig.add_axes([pos.x0 , pos.y0 + 0.025 ,  pos.width, pos.height * 0.04])  ## the parameters are the specified position you set
+        cb = fig.colorbar(dat[0],cax=cbarax, orientation="horizontal",
+                          ax = ax, ticks=np.linspace(np.log10(sig.min()),
+                         np.log10(sig.max()), 3), format="$10^{%.1f}$")
+        cb.set_label("Conductivity (S/m)",size=12)
+        cb.ax.tick_params(labelsize=12)
+
+        # Second plot for the predicted apparent resistivity data
+        ax2 = plt.subplot(2,1,2, aspect='equal')

        # Plot the location of the spheres for reference
-        circle1=plt.Circle((loc[0,0]-Tx[0][0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
-        circle2=plt.Circle((loc[0,1]-Tx[0][0,0],loc[2,1]),radi[1],color='k',fill=False, lw=3)
-        ax.add_artist(circle1)
-        ax.add_artist(circle2)
+        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
-        DC.plot_pseudoSection(survey2D,ax,stype)
+        dat = DC.plot_pseudoSection(survey2D,ax2,stype=stype, dtype = dtype)

        # 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')
-        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
+        ax2.set_title('Apparent Conductivity data')

+        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
        plt.show()

        return fig, ax
@@ -48,8 +48,7 @@ def run(plotIt=True):
    freqs = np.logspace(1,3,10)
    srcLoc = np.array([0., 0., 10.])

-    srcList = []
-    [srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
+    srcList = [EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z') for freq in freqs]

    survey = EM.FDEM.Survey(srcList)
    prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping)
@@ -0,0 +1,275 @@
+from SimPEG import *
+from SimPEG.EM import FDEM, Analytics, mu_0
+import time
+
+try:
+    from pymatsolver import MumpsSolver
+    solver = MumpsSolver
+except Exception:
+    solver = SolverLU
+    pass
+
+def run(plotIt=True):
+    """
+        EM: Schenkel and Morrison Casing Model
+        ======================================
+
+        Here we create and run a FDEM forward simulation to calculate the vertical
+        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.
+
+
+        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
+            - 300m long
+            - radius of 0.1m
+            - thickness of 6e-3m
+
+        Inside the casing, we take the same conductivity as the background.
+
+        We are using an EM code to simulate DC, so we use frequency low enough
+        that the skin depth inside the casing is longer than the casing length (f
+        = 1e-6 Hz). The plot produced is of the current inside the casing.
+
+        These results are shown in the SEG abstract by Yang et al., 2016: 3D DC
+        resistivity modeling of steel casing for reservoir monitoring using
+        equivalent resistor network. The solver used to produce these results and
+        achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_
+
+        .. _pymatsolver: https://github.com/rowanc1/pymatsolver
+
+        This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1
+
+        If you would use this example for a code comparison, or build upon it, a
+        citation would be much appreciated!
+
+    """
+
+    if plotIt:
+        import matplotlib.pylab as plt
+
+    # ------------------ MODEL ------------------
+    sigmaair = 1e-8 # air
+    sigmaback = 1e-2 # background
+    sigmacasing = 1e6 # casing
+    sigmainside = sigmaback # inside the casing
+
+
+    casing_t = 0.006  # 1cm thickness
+    casing_l = 300   # length of the casing
+
+    casing_r = 0.1
+    casing_a = casing_r - casing_t/2. # inner radius
+    casing_b = casing_r + casing_t/2. # outer radius
+    casing_z = np.r_[-casing_l,0.]
+
+
+    # ------------------ SURVEY PARAMETERS ------------------
+    freqs = np.r_[1e-6] #[1e-1, 1, 5] # frequencies
+    dsz = -300 # down-hole z source location
+    src_loc = np.r_[0.,0.,dsz]
+    inf_loc = np.r_[0.,0.,1e4]
+
+    print 'Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs]
+
+
+    # ------------------ MESH ------------------
+    # fine cells near well bore
+    csx1, csx2 = 2e-3, 60.
+    pfx1, pfx2 = 1.3, 1.3
+    ncx1 = np.ceil(casing_b/csx1+2)
+
+    # pad nicely to second cell size
+    npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1))
+    hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)])
+    dx1 = sum(hx1a)+sum(hx1b)
+    dx1 = np.floor(dx1/csx2)
+    hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b)
+
+    # second chunk of mesh
+    dx2 = 300. # uniform mesh out to here
+    ncx2 = np.ceil((dx2 - dx1)/csx2)
+    npadx2 = 45
+    hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)])
+    hx = np.hstack([hx1a,hx1b,hx2a,hx2b])
+
+    # z-direction
+    csz = 0.05
+    nza = 10
+    ncz, npadzu, npadzd = np.int(np.ceil(np.diff(casing_z)[0]/csz))+10, 68, 68 # cell size, number of core cells, number of padding cells in the x- direction
+    hz = Utils.meshTensor([(csz,npadzd,-1.3), (csz,ncz), (csz,npadzu,1.3)]) # vector of cell widths in the z-direction
+
+    # Mesh
+    mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])])
+
+    print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max())
+    print 'Number of cells', mesh.nC
+
+    if plotIt is True:
+        fig, ax = plt.subplots(1, 1, figsize=(6, 4))
+        ax.set_title('Simulation Mesh')
+        mesh.plotGrid(ax=ax)
+        plt.show()
+
+    # Put the model on the mesh
+    sigWholespace = sigmaback*np.ones((mesh.nC))
+
+    sigBack = sigWholespace.copy()
+    sigBack[mesh.gridCC[:,2] > 0.] = sigmaair
+
+    sigCasing = sigBack.copy()
+    iCasingZ = (mesh.gridCC[:,2] <= casing_z[1]) & (mesh.gridCC[:,2] >= casing_z[0])
+    iCasingX = (mesh.gridCC[:,0] >= casing_a) & (mesh.gridCC[:,0] <= casing_b)
+    iCasing = iCasingX & iCasingZ
+    sigCasing[iCasing] = sigmacasing
+
+
+    if plotIt is True:
+
+        # plotting parameters
+        xlim = np.r_[0., 0.2]
+        zlim = np.r_[-350., 10.]
+        clim_sig = np.r_[-8,6]
+
+        # plot models
+        fig, ax = plt.subplots(1,1,figsize=(4,4))
+
+        f = plt.colorbar(mesh.plotImage(np.log10(sigCasing),ax=ax)[0], ax=ax)
+        ax.grid(which='both')
+        ax.set_title('Log_10 (Sigma)')
+        ax.set_xlim(xlim)
+        ax.set_ylim(zlim)
+        f.set_clim(clim_sig)
+
+        plt.show()
+
+
+    # -------------- Sources --------------------
+    # Define Custom Current Sources
+
+    # surface source
+    sg_x = np.zeros(mesh.vnF[0],dtype=complex)
+    sg_y = np.zeros(mesh.vnF[1],dtype=complex)
+    sg_z = np.zeros(mesh.vnF[2],dtype=complex)
+
+    nza = 2 # put the wire two cells above the surface
+    ncin = 2
+
+    # vertically directed wire
+    sgv_indx = (mesh.gridFz[:,0] > casing_a) & (mesh.gridFz[:,0] < casing_a + csx1)  # hook it up to casing at the surface
+    sgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2)
+    sgv_ind = sgv_indx & sgv_indz
+    sg_z[sgv_ind] = -1.
+
+    # horizontally directed wire
+    sgh_indx = (mesh.gridFx[:,0] > casing_a)    & (mesh.gridFx[:,0] <= inf_loc[2])
+    sgh_indz = (mesh.gridFx[:,2] > csz*(nza-0.5)) & (mesh.gridFx[:,2] < csz*(nza+0.5))
+    sgh_ind = sgh_indx & sgh_indz
+    sg_x[sgh_ind] = -1.
+
+    sgv2_indx = (mesh.gridFz[:,0] >= mesh.gridFx[sgh_ind,0].max()) & (mesh.gridFz[:,0] <= inf_loc[2]*1.2)  # hook it up to casing at the surface
+    sgv2_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2)
+    sgv2_ind = sgv2_indx & sgv2_indz
+    sg_z[sgv2_ind] = 1.
+
+    # assemble the source
+    sg = np.hstack([sg_x,sg_y,sg_z])
+    sg_p = [FDEM.Src.RawVec_e([],_,sg/mesh.area) for _ in freqs]
+
+    # downhole source
+    dg_x = np.zeros(mesh.vnF[0],dtype=complex)
+    dg_y = np.zeros(mesh.vnF[1],dtype=complex)
+    dg_z = np.zeros(mesh.vnF[2],dtype=complex)
+
+    # vertically directed wire
+    dgv_indx = (mesh.gridFz[:,0] < csx1)  # go through the center of the well
+    dgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] > dsz + csz/2.)
+    dgv_ind = dgv_indx & dgv_indz
+    dg_z[dgv_ind] = -1.
+
+    # couple to the casing downhole
+    dgh_indx = mesh.gridFx[:,0] < casing_a + csx1
+    dgh_indz = (mesh.gridFx[:,2] < dsz + csz)  & (mesh.gridFx[:,2] >= dsz)
+    dgh_ind = dgh_indx & dgh_indz
+    dg_x[dgh_ind] = 1.
+
+    # horizontal part at surface
+    dgh2_indx = mesh.gridFx[:,0] <= inf_loc[2]*1.2
+    dgh2_indz = sgh_indz.copy()
+    dgh2_ind = dgh2_indx & dgh2_indz
+    dg_x[dgh2_ind] = -1.
+
+    # vertical part at surface
+    dgv2_ind = sgv2_ind.copy()
+    dg_z[dgv2_ind] = 1.
+
+    # assemble the source
+    dg = np.hstack([dg_x,dg_y,dg_z])
+    dg_p = [FDEM.Src.RawVec_e([],_,dg/mesh.area) for _ in freqs]
+
+    # ------------ Problem and Survey ---------------
+    survey = FDEM.Survey(sg_p + dg_p)
+    mapping = [('sigma', Maps.IdentityMap(mesh))]
+    problem = FDEM.Problem_h(mesh, mapping=mapping)
+    problem.pair(survey)
+
+    # ------------- Solve ---------------------------
+    t0 = time.time()
+    fieldsCasing = problem.fields(sigCasing)
+    print 'Time to solve 2 sources', time.time() - t0
+
+    # Plot current
+
+    # current density
+    jn0 = fieldsCasing[dg_p,'j']
+    jn1 = fieldsCasing[sg_p,'j']
+
+    # current
+    in0 = [mesh.area*fieldsCasing[dg_p,'j'][:,i] for i in range(len(freqs))]
+    in1 = [mesh.area*fieldsCasing[sg_p,'j'][:,i] for i in range(len(freqs))]
+
+    in0 = np.vstack(in0).T
+    in1 = np.vstack(in1).T
+
+    # integrate to get z-current inside casing
+    inds_inx = (mesh.gridFz[:,0] >= casing_a) & (mesh.gridFz[:,0] <= casing_b)
+    inds_inz = (mesh.gridFz[:,2] >= dsz )     & (mesh.gridFz[:,2] <= 0)
+    inds_fz = inds_inx & inds_inz
+
+    indsx = [False]*mesh.nFx
+    inds = list(indsx) + list(inds_fz)
+
+    in0_in = in0[np.r_[inds]]
+    in1_in = in1[np.r_[inds]]
+    z_in = mesh.gridFz[inds_fz,2]
+
+    in0_in = in0_in.reshape([in0_in.shape[0]/3,3])
+    in1_in = in1_in.reshape([in1_in.shape[0]/3,3])
+    z_in = z_in.reshape([z_in.shape[0]/3,3])
+
+    I0 = in0_in.sum(1).real
+    I1 = in1_in.sum(1).real
+    z_in = z_in[:,0]
+
+    if plotIt is True:
+        fig, ax = plt.subplots(1,2,figsize=(12,4))
+
+        ax[0].plot(z_in,np.absolute(I0), z_in,np.absolute(I1))
+        ax[0].legend(['top casing', 'bottom casing'],loc='best')
+        ax[0].set_title('Magnitude of Vertical Current in Casing')
+
+        ax[1].semilogy(z_in,np.absolute(I0), z_in,np.absolute(I1))
+        ax[1].legend(['top casing', 'bottom casing'],loc='best')
+        ax[1].set_title('Magnitude of Vertical Current in Casing')
+        ax[1].set_ylim([1e-2, 1.])
+
+        plt.show()
+
+if __name__ == '__main__':
+    run()
+
@@ -0,0 +1,132 @@
+from SimPEG import *
+
+
+def run(N=200, plotIt=True):
+    """
+        Inversion: Linear Problem
+        =========================
+
+        Here we go over the basics of creating a linear problem and inversion.
+
+    """
+
+
+    np.random.seed(1)
+
+    std_noise = 1e-2
+
+    mesh = Mesh.TensorMesh([N])
+
+    m0 = np.ones(mesh.nC) * 1e-4
+    nk = 10
+    jk = np.linspace(1.,nk,nk)
+    p = -2.
+    q = 1.
+
+    g = lambda k: np.exp(p*jk[k]*mesh.vectorCCx)*np.cos(np.pi*q*jk[k]*mesh.vectorCCx)
+
+    G = np.empty((nk, mesh.nC))
+
+    for i in range(nk):
+        G[i,:] = g(i)
+
+    mtrue = np.zeros(mesh.nC)
+    mtrue[mesh.vectorCCx > 0.3] = 1.
+    mtrue[mesh.vectorCCx > 0.45] = -0.5
+    mtrue[mesh.vectorCCx > 0.6] = 0
+
+
+    prob = Problem.LinearProblem(mesh, G)
+    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)
+
+#==============================================================================
+#     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 )
+
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
+
+    m0 = mrec
+
+    # Run inversion
+    mrec = inv.run(m0)
+
+    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
+
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+
+        fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
+        for i in range(prob.G.shape[0]):
+            axes[0].plot(prob.G[i,:])
+        axes[0].set_title('Columns of matrix G')
+
+        axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
+        axes[1].plot(mesh.vectorCCx, ml2, 'r-')
+        #axes[1].legend(('True Model', 'Recovered Model'))
+        axes[1].set_ylim(-1.0,1.25)
+
+        axes[1].plot(mesh.vectorCCx, mrec, 'k-',lw = 2)
+        axes[1].legend(('True Model', 'Smooth l2-l2',
+        'Sparse lp:' + str(reg.norms[0]) + ', lqx:' + str(reg.norms[1]) ), fontsize = 12)
+        plt.show()
+
+    return prob, survey, mesh, mrec
+
+if __name__ == '__main__':
+    run()
@@ -100,7 +100,7 @@ def run(plotIt=True):
    # Regularization - with a regularization mesh
    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
    reg = simpeg.Regularization.Tikhonov(regMesh)
-    reg.smoothModel = True
+    reg.mrefInSmooth = True
    reg.alpha_s = 1e-7
    reg.alpha_x = 1.
    # Inversion problem
@@ -5,9 +5,11 @@ import DC_Analytic_Dipole
 import DC_Forward_PseudoSection
 import EM_FDEM_1D_Inversion
 import EM_FDEM_Analytic_MagDipoleWholespace
+import EM_Schenkel_Morrison_Casing
 import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
 import Forward_BasicDirectCurrent
+import Inversion_IRLS
 import Inversion_Linear
 import Mesh_Basic_PlotImage
 import Mesh_Basic_Types
@@ -19,7 +21,7 @@ import Mesh_Tensor_Creation
 import MT_1D_ForwardAndInversion
 import MT_3D_Foward

-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "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", "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"]

 ##### AUTOIMPORTS #####

@@ -33,7 +33,9 @@ class BaseInversion(object):
        self._directiveList = value
        self._directiveList.inversion = self

-    def __init__(self, invProb, directiveList=[], **kwargs):
+    def __init__(self, invProb, directiveList=None, **kwargs):
+        if directiveList is None:
+            directiveList = []
        self.directiveList = directiveList
        Utils.setKwargs(self, **kwargs)

@@ -7,17 +7,16 @@ from SimPEG.MT.Utils.dataUtils import rec2ndarr
 # Import modules
 import numpy as np
 import os, sys, re
-try:
-    import osr
-except ImportError as e:
-    print 'Could not import osr, missing the gdal package'
-    pass
+

 class EDIimporter:
    """
    A class to import EDIfiles.

    """
+
+
+    # Define data converters
    _impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit
    _impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)

@@ -26,8 +25,8 @@ class EDIimporter:
    comps = None

    # Hidden properties
-    _outEPSG = None
-    _2out = None
+    _outEPSG = None # Project info
+    _2out = None    # The projection operator


    def __init__(self, EDIfilesList, compList=None, outEPSG=None):
@@ -113,6 +112,12 @@ class EDIimporter:
    # nOutData=length(obj.data);
    # obj.data(nOutData+1:nOutData+length(TEMP.data),:) = TEMP.data;
    def _transfromPoints(self,longD,latD):
+        # Import the coordinate projections
+        try:
+            import osr
+        except ImportError as e:
+            print 'Could not import osr, missing the gdal package\nCan not project coordinates'
+            raise e
        # Coordinates convertor
        if self._2out is None:
            src = osr.SpatialReference()
@@ -759,15 +759,29 @@ class PolyMap(IdentityMap):

            m = [\sigma_1, \sigma_2, c]

+        Can take in an actInd vector to account for topography.
+
    """
-    def __init__(self, mesh, order, logSigma=True, normal='X'):
+    def __init__(self, mesh, order, logSigma=True, normal='X', actInd = None):
        IdentityMap.__init__(self, mesh)
        self.logSigma = logSigma
        self.order = order
        self.normal = normal
+        self.actInd = actInd
+
+        if getattr(self, 'actInd', None) is None:
+            self.actInd = range(self.mesh.nC)
+            self.nC = self.mesh.nC
+
+        else:
+            self.nC = len(self.actInd)

    slope = 1e4

+    @property
+    def shape(self):
+        return (self.nC, self.nP)
+
    @property
    def nP(self):
        if np.isscalar(self.order):
@@ -785,8 +799,8 @@ class PolyMap(IdentityMap):
            sig1, sig2 = np.exp(sig1), np.exp(sig2)
        #2D
        if self.mesh.dim == 2:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
+            X = self.mesh.gridCC[self.actInd,0]
+            Y = self.mesh.gridCC[self.actInd,1]
            if self.normal =='X':
                f = polynomial.polyval(Y, c) - X
            elif self.normal =='Y':
@@ -795,9 +809,9 @@ class PolyMap(IdentityMap):
                raise(Exception("Input for normal = X or Y or Z"))
        #3D
        elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
-            Z = self.mesh.gridCC[:,2]
+            X = self.mesh.gridCC[self.actInd,0]
+            Y = self.mesh.gridCC[self.actInd,1]
+            Z = self.mesh.gridCC[self.actInd,2]
            if self.normal =='X':
                f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
            elif self.normal =='Y':
@@ -806,6 +820,7 @@ class PolyMap(IdentityMap):
                f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
            else:
                raise(Exception("Input for normal = X or Y or Z"))
+
        else:
            raise(Exception("Only supports 2D"))

@@ -819,8 +834,8 @@ class PolyMap(IdentityMap):
            sig1, sig2 = np.exp(sig1), np.exp(sig2)
        #2D
        if self.mesh.dim == 2:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
+            X = self.mesh.gridCC[self.actInd,0]
+            Y = self.mesh.gridCC[self.actInd,1]

            if self.normal =='X':
                f = polynomial.polyval(Y, c) - X
@@ -832,9 +847,9 @@ class PolyMap(IdentityMap):
                raise(Exception("Input for normal = X or Y or Z"))
        #3D
        elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]
-            Z = self.mesh.gridCC[:,2]
+            X = self.mesh.gridCC[self.actInd,0]
+            Y = self.mesh.gridCC[self.actInd,1]
+            Z = self.mesh.gridCC[self.actInd,2]

            if self.normal =='X':
                f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
@@ -330,7 +330,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
                raise NotImplementedError('wrapping in the averaging is not yet implemented')
        return self._aveF2CCV

-    def getInterpolationMatCartMesh(self, Mrect, locType='CC'):
+    def getInterpolationMatCartMesh(self, Mrect, locType='CC', locTypeTo=None):
        """
            Takes a cartesian mesh and returns a projection to translate onto the cartesian grid.
        """
@@ -338,19 +338,22 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
        assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes"


+        if locTypeTo is None:
+            locTypeTo = locType
+
        if locType == 'F':
            # do this three times for each component
-            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx')
-            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy')
-            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz')
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx', locTypeTo=locTypeTo+'x')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy', locTypeTo=locTypeTo+'y')
+            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz', locTypeTo=locTypeTo+'z')
            return sp.vstack((X,Y,Z))
        if locType == 'E':
-            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex')
-            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey')
-            Z = spzeros(Mrect.nEz, self.nE)
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex', locTypeTo=locTypeTo+'x')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey', locTypeTo=locTypeTo+'y')
+            Z = spzeros(getattr(Mrect, 'n' + locTypeTo + 'z'), self.nE)
            return sp.vstack((X,Y,Z))

-        grid = getattr(Mrect, 'grid' + locType)
+        grid = getattr(Mrect, 'grid' + locTypeTo)
        # This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice
        theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2
        theta[theta < 0] += np.pi*2.0
@@ -366,7 +369,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
                        'Ex': Mrect.tangents[:Mrect.nEx,:],
                        'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:],
                        'Ez': Mrect.tangents[-Mrect.nEz:,:],
-                    }[locType]
+                    }[locTypeTo]
            if 'F' in locType:
                normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)]
                proj = ( normals * dotMe ).sum(axis=1)
@@ -307,24 +307,28 @@ class DiffOperators(object):
        return BC
    _cellGradBC_list = 'neumann'

+    def _cellGradStencil(self):
+        BC = self.setCellGradBC(self._cellGradBC_list)
+        n = self.vnC
+        if(self.dim == 1):
+            G = ddxCellGrad(n[0], BC[0])
+        elif(self.dim == 2):
+            G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0]))
+            G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0]))
+            G = sp.vstack((G1, G2), format="csr")
+        elif(self.dim == 3):
+            G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0]))
+            G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0]))
+            G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0]))
+            G = sp.vstack((G1, G2, G3), format="csr")
+        return G
+
    def cellGrad():
        doc = "The cell centered Gradient, takes you to cell faces."

        def fget(self):
            if(self._cellGrad is None):
-                BC = self.setCellGradBC(self._cellGradBC_list)
-                n = self.vnC
-                if(self.dim == 1):
-                    G = ddxCellGrad(n[0], BC[0])
-                elif(self.dim == 2):
-                    G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0]))
-                    G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0]))
-                    G = sp.vstack((G1, G2), format="csr")
-                elif(self.dim == 3):
-                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0]))
-                    G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0]))
-                    G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0]))
-                    G = sp.vstack((G1, G2, G3), format="csr")
+                G = self._cellGradStencil()
                # Compute areas of cell faces & volumes
                S = self.area
                V = self.aveCC2F*self.vol  # Average volume between adjacent cells
@@ -361,19 +365,24 @@ class DiffOperators(object):
    _cellGradBC = None
    cellGradBC = property(**cellGradBC())

+    def _cellGradxStencil(self):
+        BC = ['neumann', 'neumann']
+        n = self.vnC
+        if(self.dim == 1):
+            G1 = ddxCellGrad(n[0], BC)
+        elif(self.dim == 2):
+            G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
+        elif(self.dim == 3):
+            G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
+        return G1
+
+
    def cellGradx():
        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."

        def fget(self):
            if getattr(self, '_cellGradx', None) is None:
-                BC = ['neumann', 'neumann']
-                n = self.vnC
-                if(self.dim == 1):
-                    G1 = ddxCellGrad(n[0], BC)
-                elif(self.dim == 2):
-                    G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
-                elif(self.dim == 3):
-                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
+                G1 = self._cellGradxStencil()
                # Compute areas of cell faces & volumes
                V = self.aveCC2F*self.vol
                L = self.r(self.area/V, 'F','Fx', 'V')
@@ -382,17 +391,22 @@ class DiffOperators(object):
        return locals()
    cellGradx = property(**cellGradx())

+    def _cellGradyStencil(self):
+        if self.dim < 2: return None
+        BC = ['neumann', 'neumann']
+        n = self.vnC
+        if(self.dim == 2):
+            G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
+        elif(self.dim == 3):
+            G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
+        return G2
+
    def cellGrady():
        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
        def fget(self):
            if self.dim < 2: return None
            if getattr(self, '_cellGrady', None) is None:
-                BC = ['neumann', 'neumann']
-                n = self.vnC
-                if(self.dim == 2):
-                    G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
-                elif(self.dim == 3):
-                    G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
+                G2 = self._cellGradyStencil()
                # Compute areas of cell faces & volumes
                V = self.aveCC2F*self.vol
                L = self.r(self.area/V, 'F','Fy', 'V')
@@ -401,14 +415,19 @@ class DiffOperators(object):
        return locals()
    cellGrady = property(**cellGrady())

+    def _cellGradzStencil(self):
+        if self.dim < 3: return None
+        BC = ['neumann', 'neumann']
+        n = self.vnC
+        G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
+        return G3
+
    def cellGradz():
        doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
        def fget(self):
            if self.dim < 3: return None
            if getattr(self, '_cellGradz', None) is None:
-                BC = ['neumann', 'neumann']
-                n = self.vnC
-                G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
+                G3 = self._cellGradzStencil()
                # Compute areas of cell faces & volumes
                V = self.aveCC2F*self.vol
                L = self.r(self.area/V, 'F','Fz', 'V')
@@ -21,10 +21,9 @@ class TensorMeshIO(object):
                if '*' in seg:
                    st = seg
                    sp = seg.split('*')
-                    re = np.array(sp[0],dtype=int)*(' ' + sp[1])
+                    re = int(sp[0])*(' ' + sp[1])
                    line = line.replace(st,re.strip())
            return np.array(line.split(),dtype=float)
-
        # Read the file as line strings, remove lines with comment = !
        msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!')

@@ -2131,10 +2131,16 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO):
    def plotSlice(self, v, vType='CC',
        normal='Z', ind=None, grid=True, view='real',
        ax=None, clim=None, showIt=False,
-        pcolorOpts={},
-        streamOpts={'color':'k'},
-        gridOpts={'color':'k', 'alpha':0.5}):
+        pcolorOpts=None,
+        streamOpts=None,
+        gridOpts=None):

+        if pcolorOpts is None:
+            pcolorOpts = {}
+        if streamOpts is None:
+            streamOpts = {'color':'k'}
+        if gridOpts is None:
+            gridOpts = {'color':'k', 'alpha':0.5}
        assert vType in ['CC','F','E']
        assert self.dim == 3

@@ -42,9 +42,9 @@ class TensorView(object):

    def plotImage(self, v, vType='CC', grid=False, view='real',
              ax=None, clim=None, showIt=False,
-              pcolorOpts={},
-              streamOpts={'color':'k'},
-              gridOpts={'color':'k'},
+              pcolorOpts=None,
+              streamOpts=None,
+              gridOpts=None,
              numbering=True, annotationColor='w'
              ):
        """
@@ -84,6 +84,12 @@ class TensorView(object):
            M.plotImage(v, annotationColor='k', showIt=True)

        """
+        if pcolorOpts is None:
+            pcolorOpts = {}
+        if streamOpts is None:
+            streamOpts = {'color':'k'}
+        if gridOpts is None:
+            gridOpts = {'color':'k'}

        if ax is None:
            fig = plt.figure()
@@ -174,9 +180,9 @@ class TensorView(object):
    def plotSlice(self, v, vType='CC',
                  normal='Z', ind=None, grid=False, view='real',
                  ax=None, clim=None, showIt=False,
-                  pcolorOpts={},
-                  streamOpts={'color':'k'},
-                  gridOpts={'color':'k', 'alpha':0.5}
+                  pcolorOpts=None,
+                  streamOpts=None,
+                  gridOpts=None
                  ):

        """
@@ -197,6 +203,12 @@ class TensorView(object):
            M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})

        """
+        if pcolorOpts is None:
+            pcolorOpts = {}
+        if streamOpts is None:
+            streamOpts = {'color':'k'}
+        if gridOpts is None:
+            gridOpts = {'color':'k', 'alpha':0.5}
        if type(vType) in [list, tuple]:
            assert ax is None, "cannot specify an axis to plot on with this function."
            fig, axs = plt.subplots(1,len(vType))
@@ -289,11 +301,17 @@ class TensorView(object):

    def _plotImage2D(self, v, vType='CC', grid=False, view='real',
              ax=None, clim=None, showIt=False,
-              pcolorOpts={},
-              streamOpts={'color':'k'},
-              gridOpts={'color':'k'}
+              pcolorOpts=None,
+              streamOpts=None,
+              gridOpts=None
              ):

+        if pcolorOpts is None:
+            pcolorOpts = {}
+        if streamOpts is None:
+            streamOpts = {'color':'k'}
+        if gridOpts is None:
+            gridOpts = {'color':'k'}
        vTypeOptsCC = ['N','CC','Fx','Fy','Ex','Ey']
        vTypeOptsV = ['CCv','F','E']
        vTypeOpts = vTypeOptsCC + vTypeOptsV
@@ -888,6 +888,8 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
    maxIterCG = 5
    tolCG = 1e-1

+    stepOffBoundsFact = 0.1 # perturbation of the inactive set off the bounds
+
    lower = -np.inf
    upper = np.inf

@@ -990,4 +992,20 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
                    cgFlag = 1
                # End CG Iterations

+            # Take a gradient step on the active cells if exist
+            if temp != self.xc.size:
+
+                rhs_a = (Active) * -self.g
+
+                dm_i = max( abs( delx ) )
+                dm_a = max( abs(rhs_a) )
+
+                # perturb inactive set off of bounds so that they are included in the step
+                delx = delx + self.stepOffBoundsFact * (rhs_a * dm_i / dm_a)
+
+
+            # Only keep gradients going in the right direction on the active set
+            indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
+            delx[indx] = 0.
+
            return delx
@@ -1,5 +1,289 @@
 import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp

+class RegularizationMesh(object):
+    """
+    **Regularization Mesh**
+
+    This contains the operators used in the regularization. Note that these
+    are not necessarily true differential operators, but are constructed from
+    a SimPEG Mesh.
+
+    :param Mesh 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
+    """
+
+    def __init__(self, mesh, indActive=None):
+        self.mesh = mesh
+        assert indActive is None or indActive.dtype == 'bool', 'indActive needs to be None or a bool'
+        self.indActive = indActive
+
+    @property
+    def vol(self):
+        """
+        reduced volume vector
+        :rtype: numpy.array
+        :return: reduced cell volume
+        """
+        if getattr(self, '_vol', None) is None:
+            self._vol = self._Pac.T * self.mesh.vol
+        return self._vol
+
+    @property
+    def nC(self):
+        """
+        reduced number of cells
+        :rtype: int
+        :return: number of cells being regularized
+        """
+        if getattr(self, '_nC', None) is None:
+            if self.indActive is None:
+                self._nC = self.mesh.nC
+            else:
+                self._nC = sum(self.indActive)
+        return self._nC
+
+    @property
+    def dim(self):
+        """
+        dimension of regularization mesh (1D, 2D, 3D)
+        :rtype: int
+        :return: dimension
+        """
+        if getattr(self, '_dim', None) is None:
+            self._dim = self.mesh.dim
+        return self._dim
+
+
+    @property
+    def _Pac(self):
+        """
+        projection matrix that takes from the reduced space of active cells to full modelling space (ie. nC x nindActive)
+        :rtype: scipy.sparse.csr_matrix
+        :return: active cell projection matrix
+        """
+        if getattr(self, '__Pac', None) is None:
+            if self.indActive is None:
+                self.__Pac = Utils.speye(self.mesh.nC)
+            else:
+                self.__Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
+        return self.__Pac
+
+    @property
+    def _Pafx(self):
+        """
+        projection matrix that takes from the reduced space of active x-faces to full modelling space (ie. nFx x nindActive_Fx )
+        :rtype: scipy.sparse.csr_matrix
+        :return: active face-x projection matrix
+        """
+        if getattr(self, '__Pafx', None) is None:
+            if self.indActive is None:
+                self.__Pafx = Utils.speye(self.mesh.nFx)
+            else:
+                indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1
+                self.__Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx]
+        return self.__Pafx
+
+    @property
+    def _Pafy(self):
+        """
+        projection matrix that takes from the reduced space of active y-faces to full modelling space (ie. nFy x nindActive_Fy )
+        :rtype: scipy.sparse.csr_matrix
+        :return: active face-y projection matrix
+        """
+        if getattr(self, '__Pafy', None) is None:
+            if self.indActive is None:
+                self.__Pafy = Utils.speye(self.mesh.nFy)
+            else:
+                indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1
+                self.__Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy]
+        return self.__Pafy
+
+    @property
+    def _Pafz(self):
+        """
+        projection matrix that takes from the reduced space of active z-faces to full modelling space (ie. nFz x nindActive_Fz )
+        :rtype: scipy.sparse.csr_matrix
+        :return: active face-z projection matrix
+        """
+        if getattr(self, '__Pafz', None) is None:
+            if self.indActive is None:
+                self.__Pafz = Utils.speye(self.mesh.nFz)
+            else:
+                indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1
+                self.__Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz]
+        return self.__Pafz
+
+    @property
+    def aveFx2CC(self):
+        """
+        averaging from active cell centers to active x-faces
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging from active cell centers to active x-faces
+        """
+        if getattr(self, '_aveFx2CC', None) is None:
+            self._aveFx2CC =  self._Pac.T * self.mesh.aveFx2CC * self._Pafx
+        return self._aveFx2CC
+
+    @property
+    def aveCC2Fx(self):
+        """
+        averaging from active x-faces to active cell centers
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging matrix from active x-faces to active cell centers
+        """
+        if getattr(self, '_aveCC2Fx', None) is None:
+            self._aveCC2Fx =  Utils.sdiag(1./(self.aveFx2CC.T).sum(1)) * self.aveFx2CC.T
+        return self._aveCC2Fx
+
+    @property
+    def aveFy2CC(self):
+        """
+        averaging from active cell centers to active y-faces
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging from active cell centers to active y-faces
+        """
+        if getattr(self, '_aveFy2CC', None) is None:
+            self._aveFy2CC = self._Pac.T * self.mesh.aveFy2CC * self._Pafy
+        return self._aveFy2CC
+
+    @property
+    def aveCC2Fy(self):
+        """
+        averaging from active y-faces to active cell centers
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging matrix from active y-faces to active cell centers
+        """
+        if getattr(self, '_aveCC2Fy', None) is None:
+            self._aveCC2Fy =  Utils.sdiag(1./(self.aveFy2CC.T).sum(1)) * self.aveFy2CC.T
+        return self._aveCC2Fy
+
+    @property
+    def aveFz2CC(self):
+        """
+        averaging from active cell centers to active z-faces
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging from active cell centers to active z-faces
+        """
+        if getattr(self, '_aveFz2CC', None) is None:
+            self._aveFz2CC = self._Pac.T * self.mesh.aveFz2CC * self._Pafz
+        return self._aveFz2CC
+
+    @property
+    def aveCC2Fz(self):
+        """
+        averaging from active z-faces to active cell centers
+        :rtype: scipy.sparse.csr_matrix
+        :return: averaging matrix from active z-faces to active cell centers
+        """
+        if getattr(self, '_aveCC2Fz', None) is None:
+            self._aveCC2Fz =  Utils.sdiag(1./(self.aveFz2CC.T).sum(1)) * self.aveFz2CC.T
+        return self._aveCC2Fz
+
+    @property
+    def cellDiffx(self):
+        """
+        cell centered difference in the x-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the x-direction
+        """
+        if getattr(self, '_cellDiffx', None) is None:
+            self._cellDiffx = self._Pafx.T * self.mesh.cellGradx * self._Pac
+        return self._cellDiffx
+
+    @property
+    def cellDiffy(self):
+        """
+        cell centered difference in the y-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the y-direction
+        """
+        if getattr(self, '_cellDiffy', None) is None:
+            self._cellDiffy = self._Pafy.T * self.mesh.cellGrady * self._Pac
+        return self._cellDiffy
+
+    @property
+    def cellDiffz(self):
+        """
+        cell centered difference in the z-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the z-direction
+        """
+        if getattr(self, '_cellDiffz', None) is None:
+            self._cellDiffz = self._Pafz.T * self.mesh.cellGradz * self._Pac
+        return self._cellDiffz
+
+    @property
+    def faceDiffx(self):
+        """
+        x-face differences
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active faces in the x-direction
+        """
+        if getattr(self, '_faceDiffx', None) is None:
+            self._faceDiffx = self._Pac.T * self.mesh.faceDivx * self._Pafx
+        return self._faceDiffx
+
+    @property
+    def faceDiffy(self):
+        """
+        y-face differences
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active faces in the y-direction
+        """
+        if getattr(self, '_faceDiffy', None) is None:
+            self._faceDiffy = self._Pac.T * self.mesh.faceDivy * self._Pafy
+        return self._faceDiffy
+
+    @property
+    def faceDiffz(self):
+        """
+        z-face differences
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active faces in the z-direction
+        """
+        if getattr(self, '_faceDiffz', None) is None:
+            self._faceDiffz = self._Pac.T * self.mesh.faceDivz * self._Pafz
+        return self._faceDiffz
+
+    @property
+    def cellDiffxStencil(self):
+        """
+        cell centered difference stencil (no cell lengths include) in the x-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the x-direction
+        """
+        if getattr(self, '_cellDiffxStencil', None) is None:
+
+            self._cellDiffxStencil = self._Pafx.T * self.mesh._cellGradxStencil() * self._Pac
+        return self._cellDiffxStencil
+
+    @property
+    def cellDiffyStencil(self):
+        """
+        cell centered difference stencil (no cell lengths include) in the y-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the y-direction
+        """
+        if self.dim < 2: return None
+        if getattr(self, '_cellDiffyStencil', None) is None:
+
+            self._cellDiffyStencil = self._Pafy.T * self.mesh._cellGradyStencil() * self._Pac
+        return self._cellDiffyStencil
+
+    @property
+    def cellDiffzStencil(self):
+        """
+        cell centered difference stencil (no cell lengths include) in the y-direction
+        :rtype: scipy.sparse.csr_matrix
+        :return: differencing matrix for active cells in the y-direction
+        """
+        if self.dim < 3: return None
+        if getattr(self, '_cellDiffzStencil', None) is None:
+
+            self._cellDiffzStencil = self._Pafz.T * self.mesh._cellGradzStencil() * self._Pac
+        return self._cellDiffzStencil
+
+
 class BaseRegularization(object):
    """
    **Base Regularization Class**
@@ -18,12 +302,19 @@ class BaseRegularization(object):

    mapping = None    #: A SimPEG.Map instance.
    mesh    = None    #: A SimPEG.Mesh instance.
-    mref = None       #: Reference model.
+    mref    = None    #: Reference model.

    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
        Utils.setKwargs(self, **kwargs)
-        self.mesh = mesh
        assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
+        if indActive is not None and indActive.dtype != 'bool':
+            tmp = indActive
+            indActive = np.zeros(mesh.nC, dtype=bool)
+            indActive[tmp] = True
+        if indActive is not None and mapping is None:
+            mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
+
+        self.regmesh = RegularizationMesh(mesh,indActive)
        self.mapping = mapping or self.mapPair(mesh)
        self.mapping._assertMatchesPair(self.mapPair)
        self.indActive = indActive
@@ -55,8 +346,7 @@ class BaseRegularization(object):
    @property
    def W(self):
        """Full regularization weighting matrix W."""
-        return sp.identity(self.mapping.nP)
-
+        return sp.identity(self.regmesh.nC)

    @Utils.timeIt
    def eval(self, m):
@@ -87,11 +377,12 @@ class BaseRegularization(object):
    @Utils.timeIt
    def eval2Deriv(self, m, v=None):
        """
+        Second derivative

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

        The regularization is:

@@ -112,112 +403,94 @@ class BaseRegularization(object):

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

-
 class Tikhonov(BaseRegularization):
    """
+    L2 Tikhonov regularization with both smallness and smoothness (first order
+    derivative) contributions.
+
+    .. math::
+        \phi_m(\mathbf{m}) = \\alpha_s \| W_s (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
+        + \\alpha_x \| W_x \\frac{\partial}{\partial x} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
+        + \\alpha_y \| W_y \\frac{\partial}{\partial y} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
+        + \\alpha_z \| W_z \\frac{\partial}{\partial z} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
+
+    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 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
+    :param float alpha_x: (default 1) smoothness weight for first derivative in the x-direction
+    :param float alpha_y: (default 1) smoothness weight for first derivative in the y-direction
+    :param float alpha_z: (default 1) smoothness weight for first derivative in the z-direction
+    :param float alpha_xx: (default 1) smoothness weight for second derivative in the x-direction
+    :param float alpha_yy: (default 1) smoothness weight for second derivative in the y-direction
+    :param float alpha_zz: (default 1) smoothness weight for second derivative in the z-direction
    """
-    smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
-    alpha_s  = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "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")
-    alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
-    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")
+    mrefInSmooth = False  # put mref in the smoothness contribution
+    alpha_s      = Utils.dependentProperty('_alpha_s', 1e-6, ['_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")
+    alpha_xx     = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'],    "Weight for the second derivative in the x direction")
+    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):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
-        self.indActive = indActive
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)

    @property
-    def Ws(self):
-        """Regularization matrix Ws"""
-        if getattr(self,'_Ws', None) is None:
-            self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Ws = Pac.T * self._Ws * Pac
-        return self._Ws
+    def Wsmall(self):
+        """Regularization matrix Wsmall"""
+        if getattr(self,'_Wsmall', None) is None:
+            self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5)
+        return self._Wsmall

    @property
    def Wx(self):
        """Regularization matrix Wx"""
        if getattr(self, '_Wx', None) is None:
-            Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol
-            self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx
-
-            if self.indActive is not None:
-                indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx]
-                self._Wx = Pafx.T*self._Wx*Pac
-
+            Ave_x_vol = self.regmesh.aveCC2Fx * self.regmesh.vol
+            self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellDiffx
        return self._Wx

    @property
    def Wy(self):
        """Regularization matrix Wy"""
        if getattr(self, '_Wy', None) is None:
-            Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol
-            self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
-
-            if self.indActive is not None:
-                indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy]
-                self._Wy = Pafy.T*self._Wy*Pac
-
+            Ave_y_vol = self.regmesh.aveCC2Fy * self.regmesh.vol
+            self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.regmesh.cellDiffy
        return self._Wy

    @property
    def Wz(self):
        """Regularization matrix Wz"""
        if getattr(self, '_Wz', None) is None:
-            Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol
-            self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
-
-            if self.indActive is not None:
-                indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz]
-                self._Wz = Pafz.T*self._Wz*Pac
-
+            Ave_z_vol = self.regmesh.aveCC2Fz * self.regmesh.vol
+            self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.regmesh.cellDiffz
        return self._Wz

    @property
    def Wxx(self):
        """Regularization matrix Wxx"""
        if getattr(self, '_Wxx', None) is None:
-            self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
-
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Wxx = Pac.T*self._Wxx*Pac
-
+            self._Wxx = Utils.sdiag((self.regmesh.vol*self.alpha_xx)**0.5)*self.regmesh.faceDiffx*self.regmesh.cellDiffx
        return self._Wxx

    @property
    def Wyy(self):
        """Regularization matrix Wyy"""
        if getattr(self, '_Wyy', None) is None:
-            self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
-
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Wyy = Pac.T*self._Wyy*Pac
-
+            self._Wyy = Utils.sdiag((self.regmesh.vol*self.alpha_yy)**0.5)*self.regmesh.faceDiffy*self.regmesh.cellDiffy
        return self._Wyy

    @property
    def Wzz(self):
        """Regularization matrix Wzz"""
        if getattr(self, '_Wzz', None) is None:
-            self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
-
-            if self.indActive is not None:
-                Pac = Utils.speye(self.mesh.nC)[:,self.indActive]
-                self._Wzz = Pac.T*self._Wzz*Pac
-
+            self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
        return self._Wzz

    @property
@@ -225,9 +498,9 @@ class Tikhonov(BaseRegularization):
        """Full smoothness regularization matrix W"""
        if getattr(self, '_Wsmooth', None) is None:
            wlist = (self.Wx, self.Wxx)
-            if self.mesh.dim > 1:
+            if self.regmesh.dim > 1:
                wlist += (self.Wy, self.Wyy)
-            if self.mesh.dim > 2:
+            if self.regmesh.dim > 2:
                wlist += (self.Wz, self.Wzz)
            self._Wsmooth = sp.vstack(wlist)
        return self._Wsmooth
@@ -236,25 +509,44 @@ class Tikhonov(BaseRegularization):
    def W(self):
        """Full regularization matrix W"""
        if getattr(self, '_W', None) is None:
-            wlist = (self.Ws, self.Wsmooth)
+            wlist = (self.Wsmall, self.Wsmooth)
            self._W = sp.vstack(wlist)
        return self._W

    @Utils.timeIt
-    def eval(self, m):
-        if self.smoothModel == True:
-            r1 = self.Wsmooth * ( self.mapping * (m) )
-            r2 = self.Ws * ( self.mapping * (m - self.mref) )
-            return 0.5*(r1.dot(r1)+r2.dot(r2))
-        elif self.smoothModel == False:
-            r = self.W * ( self.mapping * (m - self.mref) )
-            return 0.5*r.dot(r)
+    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)
+
+    @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) )

    @Utils.timeIt
    def evalDeriv(self, m):
        """
-
        The regularization is:

        .. math::
@@ -268,257 +560,185 @@ class Tikhonov(BaseRegularization):
            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}

        """
-        if self.smoothModel == True:
-            mD1 = self.mapping.deriv(m)
-            mD2 = self.mapping.deriv(m - self.mref)
-            r1 = self.Wsmooth * ( self.mapping * (m))
-            r2 = self.Ws * ( self.mapping * (m - self.mref) )
-            out1 = mD1.T * ( self.Wsmooth.T * r1 )
-            out2 = mD2.T * ( self.Ws.T * r2 )
-            out = out1+out2
-        elif self.smoothModel == False:
-            mD = self.mapping.deriv(m - self.mref)
-            r = self.W * ( self.mapping * (m - self.mref) )
-            out = mD.T * ( self.W.T * r )
-        return out
-
-# <<<<<<< HEAD
-
-# class Simple(BaseRegularization):
-#     """
-#         Only for tensor mesh
-#     """
-
-#     smoothModel = True  #: SMOOTH and SMOOTH_MOD_DIF options
-#     alpha_s  = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "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")
-#     alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
-#     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, **kwargs):
-#         BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
+        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)


+class Simple(Tikhonov):
+    """
+    Simple regularization that does not include length scales in the derivatives.
+    """

-#     @property
-#     def Ws(self):
-#         """Regularization matrix Ws"""
-#         if getattr(self,'_Ws', None) is None:
-#                 self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
-#         return self._Ws
+    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.

-#     @property
-#     def Wx(self):
-#         """Regularization matrix Wx"""
-#         if getattr(self, '_Wx', None) is None:
-#             self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.mesh.unitCellGradx
-#         return self._Wx
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)

-#     @property
-#     def Wy(self):
-#         """Regularization matrix Wy"""
-#         if getattr(self, '_Wy', None) is None:
-#             self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.mesh.unitCellGrady
-#         return self._Wy
+        if isinstance(self.wght,float):
+            self.wght = np.ones(self.regmesh.nC) * self.wght

-#     @property
-#     def Wz(self):
-#         """Regularization matrix Wz"""
-#         if getattr(self, '_Wz', None) is None:
-#             self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz
-#         return self._Wz
+    @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)
+        return self._Wsmall

-#     @property
-#     def Wxx(self):
-#         """Regularization matrix Wxx"""
-#         if getattr(self, '_Wxx', None) is None:
-#             self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
-#         return self._Wxx
+    @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
+        return self._Wx

-#     @property
-#     def Wyy(self):
-#         """Regularization matrix Wyy"""
-#         if getattr(self, '_Wyy', None) is None:
-#             self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
-#         return self._Wyy
+    @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
+        return self._Wy

-#     @property
-#     def Wzz(self):
-#         """Regularization matrix Wzz"""
-#         if getattr(self, '_Wzz', None) is None:
-#             self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
-#         return self._Wzz
+    @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
+        return self._Wz

-#     @property
-#     def Wsmooth(self):
-#         """Full smoothness regularization matrix W"""
-#         if getattr(self, '_Wsmooth', None) is None:
-#             wlist = (self.Wx, self.Wxx)
-#             if self.mesh.dim > 1:
-#                 wlist += (self.Wy, self.Wyy)
-#             if self.mesh.dim > 2:
-#                 wlist += (self.Wz, self.Wzz)
-#             self._Wsmooth = sp.vstack(wlist)
-#         return self._Wsmooth
+    @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.Ws, self.Wsmooth)
-#             self._W = sp.vstack(wlist)
-#         return self._W
+    @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 eval(self, m):
-#         if self.smoothModel == True:
-#             r1 = self.Wsmooth * ( self.mapping * (m) )
-#             r2 = self.Ws * ( self.mapping * (m - self.mref) )
-#             return 0.5*(r1.dot(r1)+r2.dot(r2))
-#         elif self.smoothModel == False:
-#             r = self.W * ( self.mapping * (m - self.mref) )
-#             return 0.5*r.dot(r)
+    @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)


-#     @Utils.timeIt
-#     def evalDeriv(self, m):
-#         """
+class Sparse(Simple):

-#         The regularization is:
+    # 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.

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

-#             R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+        if isinstance(self.wght,float):
+            self.wght = np.ones(self.regmesh.nC) * self.wght

-#         So the derivative is straight forward:
+    @property
+    def Wsmall(self):
+        """Regularization matrix Wsmall"""
+        if getattr(self, 'curModel', None) is None:
+            self.Rs = Utils.speye(self.regmesh.nC)

-#         .. math::
+        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 )

-#             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
-
-#         """
-#         if self.smoothModel == True:
-#             mD1 = self.mapping.deriv(m)
-#             mD2 = self.mapping.deriv(m - self.mref)
-#             r1 = self.Wsmooth * ( self.mapping * (m))
-#             r2 = self.Ws * ( self.mapping * (m - self.mref) )
-#             out1 = mD1.T * ( self.Wsmooth.T * r1 )
-#             out2 = mD2.T * ( self.Ws.T * r2 )
-#             out = out1+out2
-#         elif self.smoothModel == False:
-#             mD = self.mapping.deriv(m - self.mref)
-#             r = self.W * ( self.mapping * (m - self.mref) )
-#             out = mD.T * ( self.W.T * r )
-#         return out
-
-# class SparseRegularization(Simple):
-
-#     eps = 1e-1
-
-#     m = None
-#     gamma = 1.
-#     p = 0.
-#     qx = 2.
-#     qy = 2.
-#     qz = 2.
-
-#     def __init__(self, mesh, mapping=None, **kwargs):
-#         Simple.__init__(self, mesh, mapping=mapping, **kwargs)
+        return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs


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

-#     @property
-#     def W(self):
-#         """Full regularization matrix W"""
-#         if getattr(self, '_W', None) is None:
-#             wlist = (self.Ws, self.Wsmooth)
-#             self._W = sp.vstack(wlist)
-#         return self._W
+        if getattr(self, 'curModel', None) is None:
+            self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])

-#     @property
-#     def Ws(self):
-#         """Regularization matrix Ws"""
-#         if getattr(self, 'm', None) is None:
-#             self.Rs = Utils.speye(self.mesh.nC)
+        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 )

-#         else:
-#             f_m = self.m
-#             self.rs = self.R(f_m , self.p, self.eps)
-#             #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.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil

-#         self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""

-#         return self._Ws
+        if getattr(self, 'curModel', None) is None:
+            self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])

-#     @property
-#     def Wx(self):
-#         """Regularization matrix Wx"""
+        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 )

-#         if getattr(self, 'm', None) is None:
-#             self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0])
+        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

-#         else:
-#             f_m = self.mesh.unitCellGradx * self.m
-#             self.rx = self.R( f_m , self.qx, self.eps)
-#             self.Rx = Utils.sdiag( self.rx )
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""

-#         if getattr(self, '_Wx', None) is None:
-#             self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx
-#         return self._Wx
+        if getattr(self, 'curModel', None) is None:
+            self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])

-#     @property
-#     def Wy(self):
-#         """Regularization matrix Wy"""
+        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 )

-#         if getattr(self, 'm', None) is None:
-#             self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0])
+        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

-#         else:
-#             f_m = self.mesh.unitCellGrady * self.m
-#             self.ry = self.R( f_m , self.qy, self.eps)
-#             self.Ry = Utils.sdiag( self.ry )
+    @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)

-#         if getattr(self, '_Wy', None) is None:
-#             self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady
-#         return self._Wy
+    @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

-#     @property
-#     def Wz(self):
-#         """Regularization matrix Wz"""
+    def R(self, f_m , eps, exponent):

-#         if getattr(self, 'm', None) is None:
-#             self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0])
+        eta = (eps**(1.-exponent/2.))**0.5
+        r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)

-#         else:
-#             f_m = self.mesh.unitCellGradz * self.m
-#             self.rz = self.R( f_m , self.qz, self.eps)
-#             self.Rz = Utils.sdiag( self.rz )
-
-#         if getattr(self, '_Wz', None) is None:
-#             self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz
-#         return self._Wz
-
-
-#     def R(self, f_m , p, dec):
-
-#         eta = (self.eps**(1-p/2.))**0.5
-#         r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.)
-
-#         return r
-# =======
-# >>>>>>> 834de582844e8e1eac95819fbe03eed55dbeb001
+        return r
@@ -88,12 +88,14 @@ def getIndicesBlock(p0,p1,ccMesh):
    # Return a tuple
    return ind

-def defineBlock(ccMesh,p0,p1,vals=[0,1]):
+def defineBlock(ccMesh,p0,p1,vals=None):
    """
        Build a block with the conductivity specified by condVal.  Returns an array.
        vals[0]  conductivity of the block
        vals[1]  conductivity of the ground
    """
+    if vals is None:
+        vals = [0,1]
    sigma = np.zeros(ccMesh.shape[0]) + vals[1]
    ind   = getIndicesBlock(p0,p1,ccMesh)

@@ -101,7 +103,11 @@ def defineBlock(ccMesh,p0,p1,vals=[0,1]):

    return mkvc(sigma)

-def defineElipse(ccMesh, center=[0,0,0], anisotropy=[1,1,1], slope=10., theta=0.):
+def defineElipse(ccMesh, center=None, anisotropy=None, slope=10., theta=0.):
+    if center is None:
+        center = [0,0,0]
+    if anisotropy is None:
+        anisotropy = [1,1,1]
    G = ccMesh.copy()
    dim = ccMesh.shape[1]
    for i in range(dim):
@@ -156,7 +162,7 @@ def getIndicesSphere(center,radius,ccMesh):
    # Return a tuple
    return ind

-def defineTwoLayers(ccMesh,depth,vals=[0,1]):
+def defineTwoLayers(ccMesh,depth,vals=None):
    """
    Define a two layered model.  Depth of the first layer must be specified.
    CondVals vector with the conductivity values of the layers.  Eg:
@@ -167,6 +173,8 @@ def defineTwoLayers(ccMesh,depth,vals=[0,1]):
        0                          depth                                 zf
             1st layer                       2nd layer
    """
+    if vals is None:
+        vals = [0,1]
    sigma = np.zeros(ccMesh.shape[0]) + vals[1]

    dim = np.size(ccMesh[0,:])
@@ -252,7 +260,7 @@ def layeredModel(ccMesh, layerTops, layerValues):



-def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
+def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=None):
    """
        Create a random model by convolving a kernel with a
        uniformly distributed model.
@@ -276,6 +284,8 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):


    """
+    if bounds is None:
+        bounds = [0,1]

    if seed is None:
        seed = np.random.randint(1e3)
@@ -55,8 +55,10 @@ def hook(obj, method, name=None, overwrite=False, silent=False):
        print 'Method '+name+' was not overwritten.'


-def setKwargs(obj, ignore=[],  **kwargs):
+def setKwargs(obj, ignore=None,  **kwargs):
    """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
+    if ignore is None:
+        ignore = []
    for attr in kwargs:
        if attr in ignore:
            continue
@@ -0,0 +1,137 @@
+from SimPEG import np, Mesh
+import time as tm
+import vtk, vtk.util.numpy_support as npsup
+import re
+
+def read_GOCAD_ts(tsfile):
+    """
+
+    Read GOCAD triangulated surface (*.ts) file
+    INPUT:
+    tsfile: Triangulated surface
+
+    OUTPUT:
+    vrts : Array of vertices in XYZ coordinates [n x 3]
+    trgl : Array of index for triangles [m x 3]. The order of the vertices
+            is important and describes the normal
+            n = cross( (P2 - P1 ) , (P3 - P1) )
+
+    Author: @fourndo
+
+
+    .. note::
+
+        Remove all attributes from the GoCAD surface before exporting it!
+
+    """
+
+
+    fid = open(tsfile,'r')
+    line = fid.readline()
+
+    # Skip all the lines until the vertices
+    while re.match('TFACE',line)==None:
+        line = fid.readline()
+
+    line = fid.readline()
+    vrtx = []
+
+    # Run down all the vertices and save in array
+    while re.match('VRTX',line):
+        l_input  = re.split('[\s*]',line)
+        temp = np.array(l_input[2:5])
+        vrtx.append(temp.astype(np.float))
+
+        # Read next line
+        line = fid.readline()
+
+    vrtx = np.asarray(vrtx)
+
+    # Skip lines to the triangles
+    while re.match('TRGL',line)==None:
+        line = fid.readline()
+
+    # Run down the list of triangles
+    trgl = []
+
+    # Run down all the vertices and save in array
+    while re.match('TRGL',line):
+        l_input  = re.split('[\s*]',line)
+        temp = np.array(l_input[1:4])
+        trgl.append(temp.astype(np.int))
+
+        # Read next line
+        line = fid.readline()
+
+    trgl = np.asarray(trgl)
+
+    return vrtx, trgl
+
+def surface2inds(vrtx, trgl, mesh, boundaries=True, internal=True):
+    """"
+    Function to read gocad polystructure file and output indexes of mesh with in the structure.
+
+    """
+    # Adjust the index
+    trgl = trgl - 1
+
+    # Make vtk pts
+    ptsvtk = vtk.vtkPoints()
+    ptsvtk.SetData(npsup.numpy_to_vtk(vrtx,deep=1))
+
+    # Make the polygon connection
+    polys = vtk.vtkCellArray()
+    for face in trgl:
+        poly = vtk.vtkPolygon()
+        poly.GetPointIds().SetNumberOfIds(len(face))
+        for nrv, vert in enumerate(face):
+            poly.GetPointIds().SetId(nrv,vert)
+        polys.InsertNextCell(poly)
+
+    # Make the polydata, structure of connections and vrtx
+    polyData = vtk.vtkPolyData()
+    polyData.SetPoints(ptsvtk)
+    polyData.SetPolys(polys)
+
+    # Make implicit func
+    ImpDistFunc = vtk.vtkImplicitPolyDataDistance()
+    ImpDistFunc.SetInput(polyData)
+
+    # Convert the mesh
+    vtkMesh = vtk.vtkRectilinearGrid()
+    vtkMesh.SetDimensions(mesh.nNx,mesh.nNy,mesh.nNz)
+    vtkMesh.SetXCoordinates(npsup.numpy_to_vtk(mesh.vectorNx, deep=1))
+    vtkMesh.SetYCoordinates(npsup.numpy_to_vtk(mesh.vectorNy, deep=1))
+    vtkMesh.SetZCoordinates(npsup.numpy_to_vtk(mesh.vectorNz, deep=1))
+    # Add indexes
+    vtkInd = npsup.numpy_to_vtk(np.arange(mesh.nC), deep=1)
+    vtkInd.SetName('Index')
+    vtkMesh.GetCellData().AddArray(vtkInd)
+
+    extractImpDistRectGridFilt = vtk.vtkExtractGeometry() # Object constructor
+    extractImpDistRectGridFilt.SetImplicitFunction(ImpDistFunc) #
+    extractImpDistRectGridFilt.SetInputData(vtkMesh)
+
+    if boundaries is True:
+        extractImpDistRectGridFilt.ExtractBoundaryCellsOn()
+
+    else:
+        extractImpDistRectGridFilt.ExtractBoundaryCellsOff()
+
+    if internal is True:
+        extractImpDistRectGridFilt.ExtractInsideOn()
+
+    else:
+        extractImpDistRectGridFilt.ExtractInsideOff()
+
+    print "Extracting indices from grid..."
+    # Executing the pipe
+    extractImpDistRectGridFilt.Update()
+
+    # Get index inside
+    insideGrid = extractImpDistRectGridFilt.GetOutput()
+    insideGrid = npsup.vtk_to_numpy(insideGrid.GetCellData().GetArray('Index'))
+
+
+    # Return the indexes inside
+    return insideGrid
@@ -12,8 +12,16 @@
 DC Forward Simulation
 =====================

-Forward model conductive spheres in a half-space and plot a pseudo-section
+Forward model two conductive spheres in a half-space and plot a
+pseudo-section. Assumes an infinite line source and measures along the
+center of the spheres.

+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)
 Created by @fourndo on Mon Feb 01 19:28:06 2016


@@ -0,0 +1,58 @@
+.. _examples_EM_Schenkel_Morrison_Casing:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+EM: Schenkel and Morrison Casing Model
+======================================
+
+Here we create and run a FDEM forward simulation to calculate the vertical
+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.
+
+
+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
+    - 300m long
+    - radius of 0.1m
+    - thickness of 6e-3m
+
+Inside the casing, we take the same conductivity as the background.
+
+We are using an EM code to simulate DC, so we use frequency low enough
+that the skin depth inside the casing is longer than the casing length (f
+= 1e-6 Hz). The plot produced is of the current inside the casing.
+
+These results are shown in the SEG abstract by Yang et al., 2016: 3D DC
+resistivity modeling of steel casing for reservoir monitoring using
+equivalent resistor network. The solver used to produce these results and
+achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_
+
+.. _pymatsolver: https://github.com/rowanc1/pymatsolver
+
+This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1
+
+If you would use this example for a code comparison, or build upon it, a
+citation would be much appreciated!
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.EM_Schenkel_Morrison_Casing.run()
+
+.. literalinclude:: ../../SimPEG/Examples/EM_Schenkel_Morrison_Casing.py
+    :language: python
+    :linenos:
@@ -0,0 +1,26 @@
+.. _examples_Inversion_IRLS:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+Inversion: Linear Problem
+=========================
+
+Here we go over the basics of creating a linear problem and inversion.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.Inversion_IRLS.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Inversion_IRLS.py
+    :language: python
+    :linenos:
@@ -5,6 +5,8 @@ from scipy.sparse.linalg import dsolve
 import inspect

 TOL = 1e-20
+testReg = True
+testRegMesh = True

 class RegularizationTests(unittest.TestCase):

@@ -16,44 +18,80 @@ class RegularizationTests(unittest.TestCase):
        mesh3 = Mesh.TensorMesh([hx, hy, hz])
        self.meshlist = [mesh1,mesh2, mesh3]

-    def test_regularization(self):
-        for R in dir(Regularization):
-            r = getattr(Regularization, R)
-            if not inspect.isclass(r): continue
-            if not issubclass(r, Regularization.BaseRegularization):
-                continue
+    if testReg:
+        def test_regularization(self):
+            for R in dir(Regularization):
+                r = getattr(Regularization, R)
+                if not inspect.isclass(r): continue
+                if not issubclass(r, Regularization.BaseRegularization):
+                    continue
+
+                for i, mesh in enumerate(self.meshlist):
+
+                    print 'Testing %iD'%mesh.dim
+
+                    mapping = r.mapPair(mesh)
+                    reg = r(mesh, mapping=mapping)
+                    m = np.random.rand(mapping.nP)
+                    reg.mref = np.ones_like(m)*np.mean(m)
+
+                    print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                    passed = reg.eval(reg.mref) < TOL
+                    self.assertTrue(passed)
+
+                    print 'Check:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
+
+                    print 'Check 2 Deriv:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
+
+        def test_regularization_ActiveCells(self):
+            for R in dir(Regularization):
+                r = getattr(Regularization, R)
+                if not inspect.isclass(r): continue
+                if not issubclass(r, Regularization.BaseRegularization):
+                    continue
+
+                for i, mesh in enumerate(self.meshlist):
+
+                    print 'Testing Active Cells %iD'%(mesh.dim)
+
+                    if mesh.dim == 1:
+                        indActive = Utils.mkvc(mesh.gridCC <= 0.8)
+                    elif mesh.dim == 2:
+                        indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5)
+                    elif mesh.dim == 3:
+                        indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
+
+                    for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers
+                        reg = r(mesh, indActive=indAct)
+                        m = np.random.rand(mesh.nC)[indAct]
+                        reg.mref = np.ones_like(m)*np.mean(m)
+
+                    print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
+                    passed = reg.eval(reg.mref) < TOL
+                    self.assertTrue(passed)
+
+                    print 'Check:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
+
+                    print 'Check 2 Deriv:', R
+                    passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
+                    self.assertTrue(passed)
+
+    if testRegMesh:
+        def test_regularizationMesh(self):

            for i, mesh in enumerate(self.meshlist):

                print 'Testing %iD'%mesh.dim

-                mapping = r.mapPair(mesh)
-                reg = r(mesh, mapping=mapping)
-                m = np.random.rand(mapping.nP)
-                reg.mref = np.ones_like(m)*np.mean(m)
-
-                print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
-                passed = reg.eval(reg.mref) < TOL
-                self.assertTrue(passed)
-
-                print 'Check:', R
-                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
-
-                print 'Check 2 Deriv:', R
-                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
-
-    def test_regularization_ActiveCells(self):
-        for R in dir(Regularization):
-            r = getattr(Regularization, R)
-            if not inspect.isclass(r): continue
-            if not issubclass(r, Regularization.BaseRegularization):
-                continue
-
-            for i, mesh in enumerate(self.meshlist):
-
-                print 'Testing Active Cells %iD'%(mesh.dim)
+                # mapping = r.mapPair(mesh)
+                # reg = r(mesh, mapping=mapping)
+                # m = np.random.rand(mapping.nP)

                if mesh.dim == 1:
                    indAct = Utils.mkvc(mesh.gridCC <= 0.8)
@@ -62,23 +100,9 @@ class RegularizationTests(unittest.TestCase):
                elif mesh.dim == 3:
                    indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)

-                mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size)
+                regmesh = Regularization.RegularizationMesh(mesh, indActive=indAct)

-                reg = r(mesh, mapping=mapping, indActive=indAct)
-                m = np.random.rand(mesh.nC)[indAct]
-                reg.mref = np.ones_like(m)*np.mean(m)
-
-                print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref)
-                passed = reg.eval(reg.mref) < TOL
-                self.assertTrue(passed)
-
-                print 'Check:', R
-                passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
-
-                print 'Check 2 Deriv:', R
-                passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
-                self.assertTrue(passed)
+                assert (regmesh.vol == mesh.vol[indAct]).all()


 if __name__ == '__main__':
@@ -10,7 +10,9 @@ except ImportError, e:
    MumpsSolver = SolverLU


-def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,1e-3], showIt=False):
+def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None, showIt=False):
+    if bounds is None:
+        bounds = [1e-5,1e-3]
    if meshType == 'CYL':
        cs, ncx, ncz, npad = 5., 30, 10, 15
        hx = [(cs,ncx), (cs,npad,1.3)]
@@ -146,6 +146,20 @@ class TestCyl2DMesh(unittest.TestCase):

        assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3

+    def test_getInterpMatCartMesh_Cells2Nodes(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        mc = np.arange(Mc.nC)
+        xr = np.linspace(0,0.4,50)
+        xc = np.linspace(0,0.4,50) + 0.2
+        Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'N')
+        Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC')
+        Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC', locTypeTo='N')
+
+        assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
+
    def test_getInterpMatCartMesh_Faces(self):

        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
@@ -177,6 +191,37 @@ class TestCyl2DMesh(unittest.TestCase):
        assert np.abs(mag[dist > 0.1].min() - 1) < TOL


+    def test_getInterpMatCartMesh_Faces2Edges(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pf2e = Mc.getInterpolationMatCartMesh(Mr, 'F', locTypeTo='E')
+        mf = np.ones(Mc.nF)
+
+        ecart = Pf2e * mf
+
+        excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey')
+        ezcc = Mr.r(ecart, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 1) < TOL
+        assert np.abs(float(excc[indY]) - 0) < TOL
+        assert np.abs(float(eycc[indX]) - 0) < TOL
+        assert np.abs(float(eycc[indY]) - 1) < TOL
+        assert np.abs((ezcc - 1).sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
    def test_getInterpMatCartMesh_Edges(self):

        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
@@ -185,11 +230,42 @@ class TestCyl2DMesh(unittest.TestCase):
        Pe = Mc.getInterpolationMatCartMesh(Mr, 'E')
        me = np.ones(Mc.nE)

-        erect = Pe * me
+        ecart = Pe * me

-        excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex')
-        eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey')
-        ezcc = Mr.r(erect, 'E', 'Ez')
+        excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey')
+        ezcc = Mr.aveEz2CC*Mr.r(ecart, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 0) < TOL
+        assert np.abs(float(excc[indY]) + 1) < TOL
+        assert np.abs(float(eycc[indX]) - 1) < TOL
+        assert np.abs(float(eycc[indY]) - 0) < TOL
+        assert np.abs(ezcc.sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
+    def test_getInterpMatCartMesh_Edges2Faces(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pe2f = Mc.getInterpolationMatCartMesh(Mr, 'E', locTypeTo='F')
+        me = np.ones(Mc.nE)
+
+        frect = Pe2f * me
+
+        excc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx')
+        eycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy')
+        ezcc = Mr.r(frect, 'F', 'Fz')

        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
@@ -242,9 +242,6 @@ class TestAnalytics(unittest.TestCase):
    def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3))
    def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False))

-    # Do a derivative test
-    def test_derivProj1(self):self.assertTrue(DerivProjfieldsTest(halfSpace(1e-2)))
-
    # Do a derivative test of Jvec
    # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
    # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))