Added writeUBC_DCobs to StaticUtils from DCIPUtils

SimPEG/DCIP/DCIPUtils.py

@@ -1,16 +1,12 @@
-from SimPEG import np, Utils
+from SimPEG import np
 import BaseDC as DC
 import BaseDC as IP
-import warnings
 
 def getActiveindfromTopo(mesh, topo):
 # def genActiveindfromTopo(mesh, topo):
     """
         Get active indices from topography
     """
-    warnings.warn(
-            "`getActiveindfromTopo` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead",
-            FutureWarning)
     from scipy.interpolate import NearestNDInterpolator
     if mesh.dim==3:
         nCxy = mesh.nCx*mesh.nCy
@@ -32,9 +28,6 @@ def gettopoCC(mesh, airind):
     """
         Get topography from active indices of mesh.
     """
-    warnings.warn(
-            "`gettopoCC` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead",
-            FutureWarning)
     mesh2D = Mesh.TensorMesh([mesh.hx, mesh.hy], mesh.x0[:2])
     zc = mesh.gridCC[:,2]
     AIRIND = airind.reshape((mesh.vnC[0]*mesh.vnC[1],mesh.vnC[2]), order='F')
@@ -125,27 +118,34 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"):
 
 def readUBC_DC2DModel(fileName):
     """
-    Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg
+        Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg
 
-    :param string fileName: path to the UBC GIF 2D model file
-    :rtype: TensorMesh
-    :return: SimPEG TensorMesh 2D object
+        Input:
+        :param fileName, path to the UBC GIF 2D model file
+
+        Output:
+        :param SimPEG TensorMesh 2D object
+        :return
+
+        Created on Thu Nov 12 13:14:10 2015
+
+        @author: dominiquef
 
     """
     from SimPEG import np, mkvc
 
     # Open fileand skip header... assume that we know the mesh already
-    obsfile = np.genfromtxt(fileName, delimiter=' \n', dtype=np.str, comments='!')
+    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
 
-    dim = np.array(obsfile[0].split(), dtype=float)
+    dim = np.array(obsfile[0].split(),dtype=float)
 
-    temp = np.array(obsfile[1].split(), dtype=float)
+    temp = np.array(obsfile[1].split(),dtype=float)
 
     if len(temp) > 1:
         model = np.zeros(dim)
 
         for ii in range(len(obsfile)-1):
-            mm = np.array(obsfile[ii+1].split(), dtype=float)
+            mm = np.array(obsfile[ii+1].split(),dtype=float)
             model[:,ii] = mm
 
         model = model[:,::-1]
@@ -153,10 +153,10 @@ def readUBC_DC2DModel(fileName):
     else:
 
         if len(obsfile[1:])==1:
-            mm = np.array(obsfile[1:].split(), dtype=float)
+            mm = np.array(obsfile[1:].split(),dtype=float)
 
         else:
-            mm = np.array(obsfile[1:], dtype=float)
+            mm = np.array(obsfile[1:],dtype=float)
 
         # Permute the second dimension to flip the order
         model = mm.reshape(dim[1],dim[0])
@@ -169,19 +169,23 @@ def readUBC_DC2DModel(fileName):
 
     return model
 
-
-def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt', clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None):
+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.
+        Read list of 2D tx-rx location and plot a speudo-section of apparent
+        resistivity.
 
-    Assumes flat topo for now...
+        Assumes flat topo for now...
 
-    :param SurveyDC DCsurvey:
-    :param string surveyType: Either 'pole-dipole'  | 'dipole-dipole'
-    :param string unitType: Either 'appResistivity' | 'appConductivity'  | 'volt'
-    :rtype: matplotlib.plt
-    :return: figure scatter plot overlayed on image
+        Input:
+        :param d2D, z0
+        :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole)
+        :switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential)
+        Output:
+        :figure scatter plot overlayed on image
+
+        Edited Feb 17th, 2016
+
+        @author: dominiquef
 
     """
     from SimPEG import np
@@ -214,39 +218,39 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
         Cmid = (Tx[0][0] + Tx[1][0])/2
         Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
 
-        # Change output for unitType
-        if unitType == 'volt':
+        # Change output for dtype
+        if dtype == 'volt':
 
             rho = np.hstack([rho,data])
 
         else:
 
             # Compute pant leg of apparent rho
-            if surveyType == 'pole-dipole':
+            if stype == 'pdp':
 
                 leg =  data * 2*np.pi  * MA * ( MA + MN ) / MN
 
-            elif surveyType == 'dipole-dipole':
+            elif stype == 'dpdp':
 
                 leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
 
             else:
-                print """unitType must be 'pole-dipole' | 'dipole-dipole' """
+                print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """
                 break
 
 
-            if unitType == 'appConductivity':
+            if dtype == 'appc':
 
                 leg = np.log10(abs(1./leg))
                 rho = np.hstack([rho,leg])
 
-            elif unitType == 'appResistivity':
+            elif dtype == 'appr':
 
                 leg = np.log10(abs(leg))
                 rho = np.hstack([rho,leg])
 
             else:
-                print """unitType must be 'appResistivity' | 'appConductivity' | 'volt' """
+                print """dtype must be 'appr' | 'appc' | 'volt' """
                 break
 
         midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
@@ -255,7 +259,7 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
     # 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()
@@ -264,39 +268,38 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
 
     # 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 unitType == 'volt':
+        
+        if dtype == 'volt':
             cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal")
 
-        else:
+        else:               
             cbar = plt.colorbar(ph, ax = axs, format="$10^{%.1f}$",fraction=0.04,orientation="horizontal")
-
-            cmin,cmax = cbar.get_clim()
-            ticks = np.linspace(cmin,cmax,3)
-            cbar.set_ticks(ticks)
-            cbar.ax.tick_params(labelsize=10)
-
-            if unitType == 'appConductivity':
+        
+        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 unitType == 'appResistivity':
+            elif dtype == 'appr':
                 cbar.set_label("App.Res.",size=12)
-            elif unitType == 'volt':
+            elif dtype == 'volt':
                 cbar.set_label("Potential (V)",size=12)
 
 
+
     if not axlabel:
         axs.set_xticklabels([])
         axs.set_yticklabels([])
@@ -307,24 +310,27 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt
 
     return ph
 
-def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
+def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
     """
-    Load in endpoints and survey specifications to generate Tx, Rx location
-    stations.
+        Load in endpoints and survey specifications to generate Tx, Rx location
+        stations.
 
-    Assumes flat topo for now...
+        Assumes flat topo for now...
 
-    :param numpy.array endl: input endpoints [[x1, y1] , [x2, y2]]
-    :param Mesh mesh: SimPEG mesh object
-    :param string surveyType: 'dipole-dipole' | 'pole-dipole' | 'gradient'
-    :param float AM_sep: transmitter (A) - receiver (M) seperation
-    :param float b: receiver dipole seperation
-    :param float nrx: pole seperation, number of rx dipoles per tx
+        Input:
+        :param endl -> input endpoints [x1, y1, z1, x2, y2, z2]
+        :object mesh -> SimPEG mesh object
+        :switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient'
+        : param a, n -> pole seperation, number of rx dipoles per tx
 
-    :rtype: DC.Survey, Src, Rx
-    :returns: DC survey, Source
+        Output:
+        :param Tx, Rx -> List objects for each tx location
+            Lines: P1x, P1y, P1z, P2x, P2y, P2z
 
-    !! Require clean up to deal with DCsurvey
+        Created on Wed December 9th, 2015
+
+        @author: dominiquef
+        !! Require clean up to deal with DCsurvey
     """
 
     from SimPEG import np
@@ -340,17 +346,17 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
     dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
     dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
 
-    nstn = np.floor( dl_len / AM_sep )
+    nstn = np.floor( dl_len / a )
 
     # Compute discrete pole location along line
-    stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*AM_sep
-    stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*AM_sep
+    stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
+    stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
 
     # Create line of P1 locations
     M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
 
     # Create line of P2 locations
-    N = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
+    N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
 
     ## Build list of Tx-Rx locations depending on survey type
     # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn]
@@ -360,14 +366,14 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
     SrcList = []
 
 
-    if surveyType != 'gradient':
+    if stype != 'gradient':
 
         for ii in range(0, int(nstn)-1):
 
 
-            if surveyType == 'dipole-dipole':
+            if stype == 'dpdp':
                 tx = np.c_[M[ii,:],N[ii,:]]
-            elif surveyType == 'pole-dipole':
+            elif stype == 'pdp':
                 tx = np.c_[M[ii,:],M[ii,:]]
 
             # Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]])
@@ -376,33 +382,33 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
             AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
 
             # Number of receivers to fit
-            nstn = np.min([np.floor( (AB - MN_sep) / AM_sep ) , nrx])
+            nstn = np.min([np.floor( (AB - b) / a ) , n])
 
             # Check if there is enough space, else break the loop
             if nstn <= 0:
                 continue
 
             # Compute discrete pole location along line
-            stn_x = N[ii,0] + dl_x*MN_sep + np.array(range(int(nstn)))*dl_x*AM_sep
-            stn_y = N[ii,1] + dl_y*MN_sep + np.array(range(int(nstn)))*dl_y*AM_sep
+            stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
+            stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
 
             # Create receiver poles
             # Create line of P1 locations
             P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
 
             # Create line of P2 locations
-            P2 = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
+            P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
 
             Rx.append(np.c_[P1,P2])
             rxClass = DC.RxDipole(P1, P2)
             Tx.append(tx)
-            if surveyType == 'dipole-dipole':
+            if stype == 'dpdp':
                 srcClass = DC.SrcDipole([rxClass], M[ii,:],N[ii,:])
-            elif surveyType == 'pole-dipole':
+            elif stype == 'pdp':
                 srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:])
             SrcList.append(srcClass)
 
-    elif surveyType == 'gradient':
+    elif stype == 'gradient':
 
         # Gradient survey only requires Tx at end of line and creates a square
         # grid of receivers at in the middle at a pre-set minimum distance
@@ -410,23 +416,23 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
         Tx.append(np.c_[M[0,:],N[-1,:]])
 
         # Get the edge limit of survey area
-        min_x = endl[0,0] + dl_x * MN_sep
-        min_y = endl[0,1] + dl_y * MN_sep
+        min_x = endl[0,0] + dl_x * b
+        min_y = endl[0,1] + dl_y * b
 
-        max_x = endl[1,0] - dl_x * MN_sep
-        max_y = endl[1,1] - dl_y * MN_sep
+        max_x = endl[1,0] - dl_x * b
+        max_y = endl[1,1] - dl_y * b
 
         box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
         box_w = box_l/2.
 
-        nstn = np.floor( box_l / AM_sep )
+        nstn = np.floor( box_l / a )
 
         # Compute discrete pole location along line
-        stn_x = min_x + np.array(range(int(nstn)))*dl_x*AM_sep
-        stn_y = min_y + np.array(range(int(nstn)))*dl_y*AM_sep
+        stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
+        stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
 
         # Define number of cross lines
-        nlin = int(np.floor( box_w / AM_sep ))
+        nlin = int(np.floor( box_w / a ))
         lind = range(-nlin,nlin+1)
 
         ngrad = nstn * len(lind)
@@ -435,12 +441,12 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
         for ii in range( len(lind) ):
 
             # Move line in perpendicular direction by dipole spacing
-            lxx = stn_x - lind[ii]*AM_sep*dl_y
-            lyy = stn_y + lind[ii]*AM_sep*dl_x
+            lxx = stn_x - lind[ii]*a*dl_y
+            lyy = stn_y + lind[ii]*a*dl_x
 
 
             M = np.c_[ lxx, lyy , np.ones(nstn).T*mesh.vectorNz[-1]]
-            N = np.c_[ lxx+AM_sep*dl_x, lyy+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
+            N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
 
             rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N]
 
@@ -449,38 +455,44 @@ def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx):
         srcClass = DC.SrcDipole([rxClass], M[0,:], N[-1,:])
         SrcList.append(srcClass)
     else:
-        print """surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """
+        print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
 
     survey = DC.SurveyDC(SrcList)
     return survey, Tx, Rx
 
-
-def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
+def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
     """
         Write UBC GIF DCIP 2D or 3D observation file
 
-        :param string fileName: including path where the file is written out
-        :param Survey DCsurvey: DC survey class object
-        :param string dim:  either '2D' | '3D'
-        :param string surveyType:  either 'SURFACE' | 'GENERAL'
-        :rtype: file
-        :return: UBC2D-Data file
-    """
+        Input:
+        :string     fileName -> including path where the file is written out
+        :DCsurvey   DC survey class object
+        :string     dtype ->  either '2D' | '3D'
+        :string     stype ->  either 'SURFACE' | 'GENERAL'
 
+        Output:
+        :param UBC2D-Data file
+        :return
+
+        Last edit: February 16th, 2016
+
+        @author: dominiquef
+
+    """
     from SimPEG import mkvc
 
-    assert (dim=='2D') | (dim=='3D'), "Data must be either '2D' | '3D'"
-    assert (surveyType=='SURFACE') | (surveyType=='GENERAL') | (surveyType=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
-
+    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('! ' + surveyType + ' FORMAT\n')
-
+    
+    
     if iptype!=0:
         fid.write('IPTYPE=%i\n'%iptype)
-
+        
     else:
         fid.write('! ' + stype + ' FORMAT\n')
-
+        
     count = 0
 
     for ii in range(DCsurvey.nSrc):
@@ -494,33 +506,33 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
         M = rx[0]
         N = rx[1]
 
-        # Adapt source-receiver location for dim and surveyType
-        if dim=='2D':
+        # Adapt source-receiver location for dtype and stype
+        if dtype=='2D':
 
-            if surveyType == 'SIMPLE':
+            if stype == 'SIMPLE':
 
                 #fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
                 A = np.repeat(tx[0,0],M.shape[0],axis=0)
                 B = np.repeat(tx[0,1],M.shape[0],axis=0)
                 M = M[:,0]
                 N = N[:,0]
-
+                
                 np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
 
 
             else:
 
-                if surveyType == 'SURFACE':
+                if stype == 'SURFACE':
 
                     fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
                     M = M[:,0]
                     N = N[:,0]
 
-                if surveyType == 'GENERAL':
+                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]
@@ -528,31 +540,31 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType, iptype = 0):
                     # Flip sign for z-elevation to depth
                     M[:,1::2] = -M[:,1::2]
                     N[:,1::2] = -N[:,1::2]
-
+                    
                 fid.write('%i\n'% nD)
                 np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
 
-        if dim=='3D':
+        if dtype=='3D':
 
-            if surveyType == 'SURFACE':
+            if stype == 'SURFACE':
 
                 fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:]))
                 M = M[:,0:2]
                 N = N[:,0:2]
 
-            if surveyType == 'GENERAL':
+            if stype == 'GENERAL':
 
                 fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
 
             fid.write('%i\n'% nD)
             np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
             fid.write('\n')
-
+            
         count += nD
 
     fid.close()
 
-def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'):
+def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
     """
         Read DC survey and projects the coordinate system
         according to the flag = 'Xloc' | 'Yloc' | 'local' (default)
@@ -561,9 +573,15 @@ def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'):
 
         The Z value is preserved, but Y coordinates zeroed.
 
-        :param DC.Survey survey3D: 3D simpeg DC survey
-        :rtype: DC.Survey
-        :return: survey2D
+        Input:
+        :param survey3D
+
+        Output:
+        :figure survey2D
+
+        Edited April 6th, 2016
+
+        @author: dominiquef
 
     """
     from SimPEG import np
@@ -648,34 +666,39 @@ def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'):
     DCsurvey2D.std = np.asarray(DCsurvey.std)
 
     return DCsurvey2D
-
-def readUBC_DC3Dobs(fileName, rtype = 'DC'):
+    
+def readUBC_DC3Dobs(fileName, dtype = 'DC'):
     """
         Read UBC GIF IP 3D observation file and generate survey
 
-        :param string fileName:, path to the UBC GIF 3D obs file
-        :rtype: Survey
-        :return: DCIPsurvey
+        Input:
+        :param fileName, path to the UBC GIF 3D obs file
+
+        Output:
+        :param IPsurvey
+        :return
+
+        @author: dominiquef
 
     """
     zflag = True # Flag for z value provided
-
+    
     # Load file
-    if rtype == 'IP':
+    if dtype == 'IP':
         obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='IPTYPE')
-
-    elif rtype == 'DC':
+        
+    elif dtype == 'DC':
         obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
-
+        
     else:
-        print "rtype must be 'DC'(default) | 'IP'"
-
+        print "dtype must be 'DC'(default) | 'IP'"
+        
     # Pre-allocate
     srcLists = []
     Rx = []
     d = []
     wd = []
-
+    
 
     # Countdown for number of obs/tx
     count = 0
@@ -694,7 +717,7 @@ def readUBC_DC3Dobs(fileName, rtype = 'DC'):
             # Check if z value is provided, if False -> nan
             if len(temp)==5:
                 tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan]
-
+                
                 zflag = False # Pass on the flag to the receiver loc
 
             else:
@@ -706,12 +729,12 @@ def readUBC_DC3Dobs(fileName, rtype = 'DC'):
         temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string
 
         # Filter out negative IP
-#        if temp[-2] < 0:
+#        if temp[-2] < 0:       
 #            count = count -1
 #            print "Negative!"
-#
+#            
 #        else:
-
+            
         # If the Z-location is provided, otherwise put nan
         if zflag:
 
@@ -749,9 +772,17 @@ def readUBC_DC2Dobs(fileName):
         ------- NEEDS TO BE UPDATED ------
         Read UBC GIF 2D observation file and generate arrays for tx-rx location
 
-        :param string fileName: path to the UBC GIF 2D model file
-        :rtype: (DC.Src, DC.Rx, ??, ??)
-        :return: source_locs, rx_locs, ??, ??
+        Input:
+        :param fileName, path to the UBC GIF 2D model file
+
+        Output:
+        :param rx, tx
+        :return
+
+        Created on Thu Nov 12 13:14:10 2015
+
+        @author: dominiquef
+
     """
 
     from SimPEG import np
@@ -791,9 +822,11 @@ def readUBC_DC2Dpre(fileName):
         Read UBC GIF DCIP 2D observation file and generate arrays for tx-rx location
 
         Input:
-        :param string fileName: path to the UBC GIF 3D obs file
-        :rtype: DC.Survey
-        :return: DCsurvey
+        :param fileName, path to the UBC GIF 3D obs file
+
+        Output:
+        DCsurvey
+        :return
 
         Created on Mon March 9th, 2016 << Doug's 70th Birthday !! >>
 
@@ -855,9 +888,12 @@ def readUBC_DC2DMesh(fileName):
     """
         Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg
 
-        :param string fileName: path to the UBC GIF mesh file
-        :rtype: Mesh.TensorMesh
-        :return: SimPEG TensorMesh 2D object
+        Input:
+        :param fileName, path to the UBC GIF mesh file
+
+        Output:
+        :param SimPEG TensorMesh 2D object
+        :return
 
         Created on Thu Nov 12 13:14:10 2015
 
@@ -923,9 +959,12 @@ def xy_2_lineID(DCsurvey):
         they were collected. May need to generalize for random
         point locations, but will be more expensive
 
-        :param numpy.array DCdict: Vectors of station location
-        :rtype: numpy.array
-        :return: LineID Vector of integers
+        Input:
+        :param DCdict Vectors of station location
+
+        Output:
+        :param LineID Vector of integers
+        :return
 
         Created on Thu Feb 11, 2015
 

SimPEG/EM/Base.py

@@ -169,7 +169,9 @@ class BaseEMProblem(Problem.BaseProblem):
 
         dMeSigmaI_dI = -self.MeSigmaI**2
         dMe_dsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u)
-        return dMeSigmaI_dI * ( dMe_dsig * self.curModel.sigmaDeriv )
+        dsig_dm = self.curModel.sigmaDeriv
+        return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm))
+        # return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u)
 
     @property
     def MfRho(self):
@@ -185,7 +187,8 @@ class BaseEMProblem(Problem.BaseProblem):
         """
             Derivative of :code:`MfRho` with respect to the model.
         """
-        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv
+        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv)
+        # self.curModel.rhoDeriv
 
     @property
     def MfRhoI(self):
@@ -205,7 +208,9 @@ class BaseEMProblem(Problem.BaseProblem):
 
         dMfRhoI_dI = -self.MfRhoI**2
         dMf_drho = self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u)
-        return dMfRhoI_dI * ( dMf_drho * self.curModel.rhoDeriv )
+        return dMfRhoI_dI * ( dMf_drho * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) )
+
+        # return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
 
 class BaseEMSurvey(Survey.BaseSurvey):
 

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -257,7 +257,7 @@ class Fields3D_e(Fields):
         """
 
         # assuming primary does not depend on the model
-        return src.ePrimaryDeriv(self.prob, v, adjoint) #Zero()
+        return Zero()
 
     def _bPrimary(self, eSolution, srcList):
         """
@@ -600,8 +600,8 @@ class Fields3D_b(Fields):
 
 
         if adjoint:
-            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv + src.ePrimaryDeriv(self.prob, v, adjoint)
-        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv + src.ePrimaryDeriv(self.prob, v, adjoint)
+            return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv
+        return  self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv
 
     def _j(self, bSolution, srcList):
         """

SimPEG/EM/FDEM/ProblemFDEM.py

@@ -74,8 +74,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         self.curModel = m
 
-        # Jv = self.dataPair(self.survey)
-        Jv = []
+        Jv = self.dataPair(self.survey)
 
         for freq in self.survey.freqs:
             A = self.getA(freq)
@@ -90,9 +89,9 @@ class BaseFDEMProblem(BaseEMProblem):
                 for rx in src.rxList:
                     df_dmFun = getattr(f, '_{0}Deriv'.format(rx.projField), None)
                     df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
-                    Jv.append(rx.evalDeriv(src, self.mesh, f, df_dm_v))
+                    Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
             Ainv.clean()
-        return np.hstack(Jv)
+        return Utils.mkvc(Jv)
 
     def Jtvec(self, m, v, f=None):
         """
@@ -167,6 +166,7 @@ class BaseFDEMProblem(BaseEMProblem):
 
         for i, src in enumerate(Srcs):
             smi, sei = src.eval(self)
+            #Why are you adding?
             s_m[:,i] = s_m[:,i] + smi
             s_e[:,i] = s_e[:,i] + sei
 

SimPEG/EM/FDEM/SrcFDEM.py

@@ -60,18 +60,6 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._bPrimary
 
-    def bPrimaryDeriv(self, prob, v, adjoint=False):
-        """
-        Derivative of the primary magnetic flux density
-
-        :param Problem prob: FDEM Problem
-        :param numpy.ndarray v: vector
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: primary magnetic flux density
-        """
-        return Zero()
-
     def hPrimary(self, prob):
         """
         Primary magnetic field
@@ -84,18 +72,6 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._hPrimary
 
-    def hPrimaryDeriv(self, prob, v, adjoint=False):
-        """
-        Derivative of the primary magnetic field
-
-        :param Problem prob: FDEM Problem
-        :param numpy.ndarray v: vector
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: primary magnetic flux density
-        """
-        return Zero()
-
     def ePrimary(self, prob):
         """
         Primary electric field
@@ -108,18 +84,6 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._ePrimary
 
-    def ePrimaryDeriv(self, prob, v, adjoint=False):
-        """
-        Derivative of the primary electric field
-
-        :param Problem prob: FDEM Problem
-        :param numpy.ndarray v: vector
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: primary magnetic flux density
-        """
-        return Zero()
-
     def jPrimary(self, prob):
         """
         Primary current density
@@ -132,18 +96,6 @@ class BaseSrc(Survey.BaseSrc):
             return Zero()
         return self._jPrimary
 
-    def jPrimaryDeriv(self, prob, v, adjoint=False):
-        """
-        Derivative of the primary current density
-
-        :param Problem prob: FDEM Problem
-        :param numpy.ndarray v: vector
-        :param bool adjoint: adjoint?
-        :rtype: numpy.ndarray
-        :return: primary magnetic flux density
-        """
-        return Zero()
-
     def s_m(self, prob):
         """
         Magnetic source term
@@ -603,7 +555,7 @@ class CircularLoop(BaseSrc):
             a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
-            srcfct = MagneticLoopVectorPotential
+            srcfct = MagneticDipoleVectorPotential
             ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
             ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
             az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)
@@ -662,155 +614,5 @@ class CircularLoop(BaseSrc):
             return -C.T * (MMui_s * self.bPrimary(prob))
 
 
-class PrimSecSigma(BaseSrc):
-
-    def __init__(self, rxList, freq, sigBack, ePrimary, **kwargs):
-        self.sigBack = sigBack
-
-        BaseSrc.__init__(self, rxList, freq=freq, _ePrimary=ePrimary, **kwargs)
-
-    def s_e(self, prob):
-        return (prob.MeSigma -  prob.mesh.getEdgeInnerProduct(self.sigBack)) * self.ePrimary(prob)
-
-    def s_eDeriv(self, prob, v, adjoint=False):
-        if adjoint:
-            return prob.MeSigmaDeriv(self.ePrimary(prob)).T * v
-        return prob.MeSigmaDeriv(self.ePrimary(prob)) * v
-
-
-class PrimSecMappedSigma(BaseSrc):
-
-    """
-    Primary-Secondary Source in which a mapping is provided to put the current model
-    onto the primary mesh. This is solved on every model update.
-
-    There are a lot of layers to the derivatives here!
-
-    **Required**
-    :param list rxList: Receiver List
-    :param float freq: frequency
-    :param ProblemFDEM primaryProblem: FDEM primary problem
-    :param SurveyFDEM primarySurvey: FDEM primary survey
-
-    **Optional**
-    :param Mapping map2meshSecondary: mapping current model to act as primary model on the secondary mesh
-
-    """
-
-    def __init__(self, rxList, freq, primaryProblem, primarySurvey, map2meshSecondary = None ,**kwargs):
-
-        self.primaryProblem = primaryProblem
-        self.primarySurvey = primarySurvey
-
-        if self.primaryProblem.ispaired is False:
-            self.primaryProblem.pair(self.primarySurvey)
-
-        self.map2meshSecondary = map2meshSecondary
-
-        BaseSrc.__init__(self, rxList, freq=freq, **kwargs)
-
-    def _ProjPrimary(self, prob):
-        # if getattr(self, '__ProjPrimary', None) is None:
-        return self.primaryProblem.mesh.getInterpolationMatCartMesh(prob.mesh, locType='F', locTypeTo='E')
-        # return self.__ProjPrimary
-
-
-    def _primaryFields(self, prob, fieldType=None):
-
-        # TODO: cache and check if prob.curModel has changed
-        fields = self.primaryProblem.fields(prob.curModel.sigmaModel)
-
-        if fieldType is not None:
-            return fields[:,fieldType]
-        return fields
-
-    def _primaryFieldsDeriv(self, prob, v, adjoint=False, f=None):
-        if adjoint:
-            raise NotImplementedError
-
-        # TODO: this should not be hard-coded for j
-        # jp = self._primaryFields(prob)[:,'j']
-
-        # TODO: pull apart Jvec so that don't have to copy paste this code in
-        # A = self.primaryProblem.getA(self.freq)
-        # Ainv = self.primaryProblem.Solver(A, **self.primaryProblem.solverOpts) # create the concept of Ainv (actually a solve)
-
-        if f is None:
-           f = self._primaryFields(prob.curModel.sigmaModel)
-
-        freq = self.freq
-
-        A = self.primaryProblem.getA(freq)
-        Ainv = self.primaryProblem.Solver(A, **self.primaryProblem.solverOpts) # create the concept of Ainv (actually a solve)
-
-        src = self.primarySurvey.srcList[0]
-        # for src in self.survey.getSrcByFreq(freq):
-        u_src = Utils.mkvc(f[src, self.primaryProblem._solutionType])
-        dA_dm_v = self.primaryProblem.getADeriv(freq, u_src, v)
-        dRHS_dm_v = self.primaryProblem.getRHSDeriv(freq, src, v)
-        du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
-
-        df_dmFun = getattr(f, '_{0}Deriv'.format('j'), None)
-        df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
-        # Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
-        Ainv.clean()
-
-        return df_dm_v
-
-        # return self.primaryProblem.Jvec(prob.curModel, v, f=f)
-
-    def ePrimary(self, prob, f=None):
-        if f is None:
-            f = self._primaryFields(prob)
-
-        ep = self._ProjPrimary(prob) * (
-                self.primaryProblem.MfI * (
-                self.primaryProblem.MfRho * f[:,'j'])
-                )
-
-        return Utils.mkvc(ep)
-
-    def ePrimaryDeriv(self, prob, v, adjoint=False, f=None):
-
-        if adjoint is True:
-            raise NotImplementedError
-
-        if f is None:
-            f = self._primaryFields(prob)
-
-        epDeriv = self._ProjPrimary(prob) * (
-                    self.primaryProblem.MfI * (
-                    (self.primaryProblem.MfRhoDeriv(f[:,'j']) * v)
-                    +
-                    (self.primaryProblem.MfRho * self._primaryFieldsDeriv(prob, v, f=f))
-                    )
-                    )
-
-        return Utils.mkvc(epDeriv)
-
-
-    def s_e(self, prob):
-        sigmaPrimary = self.map2meshSecondary * prob.curModel.sigmaModel
-
-        return Utils.mkvc((prob.MeSigma -  prob.mesh.getEdgeInnerProduct(sigmaPrimary)) * self.ePrimary(prob))
-
-
-    def s_eDeriv(self, prob, v, adjoint=False):
-        if adjoint:
-            raise NotImplementedError
-            return prob.MeSigmaDeriv(self.ePrimary(prob)).T * v
-
-        sigmaPrimary = self.map2meshSecondary * prob.curModel.sigmaModel
-        sigmaPrimaryDeriv = self.map2meshSecondary.deriv(prob.curModel.sigmaModel)
-
-        f = self._primaryFields(prob)
-        ePrimary = self.ePrimary(prob,f=f)
-
-        return (prob.MeSigmaDeriv(ePrimary) * v
-                - prob.mesh.getEdgeInnerProductDeriv(sigmaPrimary)(ePrimary) * sigmaPrimaryDeriv * v
-                + (prob.MeSigma -  prob.mesh.getEdgeInnerProduct(sigmaPrimary)) * self.ePrimaryDeriv(prob, v, None, f=f)
-                )
-
-
 
 

SimPEG/EM/Static/DC/ProblemDC.py

@@ -35,10 +35,11 @@ class BaseDCProblem(BaseEMProblem):
 
         self.curModel = m
 
-        Jv = self.dataPair(self.survey) #same size as the data
-
+        # Jv = self.dataPair(self.survey) #same size as the data
         A = self.getA()
 
+        Jv = []
+
         for src in self.survey.srcList:
             u_src = f[src, self._solutionType] # solution vector
             dA_dm_v = self.getADeriv(u_src, v)
@@ -48,8 +49,10 @@ class BaseDCProblem(BaseEMProblem):
             for rx in src.rxList:
                 df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
                 df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
-                Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
-        return Utils.mkvc(Jv)
+                # Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
+                Jv.append(rx.evalDeriv(src, self.mesh, f, df_dm_v))
+        # return Utils.mkvc(Jv)
+        return np.hstack(Jv)
 
     def Jtvec(self, m, v, f=None):
         if f is None:
@@ -64,7 +67,6 @@ class BaseDCProblem(BaseEMProblem):
         Jtv = np.zeros(m.size)
         AT = self.getA()
 
-
         for src in self.survey.srcList:
             u_src = f[src, self._solutionType]
             for rx in src.rxList:

SimPEG/EM/Static/DC/RxDC.py

@@ -43,14 +43,7 @@ class BaseRx(SimPEG.Survey.BaseRx):
         elif adjoint:
             return P.T*v
 
-# DC.Rx.Pole(locs)
-class Pole(BaseRx):
-
-    def __init__(self, locs, rxType = 'phi', **kwargs):
-        BaseRx.__init__(self, locs, rxType)
-
-
-# DC.Rx.Dipole(locsM, locsN)
+# DC.Rx.Dipole(locs)
 class Dipole(BaseRx):
 
     def __init__(self, locsM, locsN, rxType = 'phi', **kwargs):

SimPEG/EM/Static/IP/ProblemIP.py

@@ -45,7 +45,8 @@ class BaseIPProblem(BaseEMProblem):
 
         self.curModel = m
 
-        Jv = self.dataPair(self.survey) #same size as the data
+        # Jv = self.dataPair(self.survey) #same size as the data
+        Jv = []
 
         A = self.getA()
 
@@ -58,13 +59,16 @@ class BaseIPProblem(BaseEMProblem):
             for rx in src.rxList:
                 df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
                 df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
-                Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
+                # Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
+                Jv.append(rx.evalDeriv(src, self.mesh, f, df_dm_v))
         # Conductivity (d u / d log sigma)
         if self._formulation is 'EB':
-            return -Utils.mkvc(Jv)
+            # return -Utils.mkvc(Jv)
+            return -np.hstack(Jv)
         # Conductivity (d u / d log rho)
         if self._formulation is 'HJ':
-            return Utils.mkvc(Jv)
+            # return Utils.mkvc(Jv)
+            return np.hstack(Jv)
 
     def Jtvec(self, m, v, f=None):
         if f is None:

SimPEG/EM/Static/Utils/StaticUtils.py

@@ -315,3 +315,107 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
 
     return SrcList
 
+
+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'
+
+        Output:
+        :param UBC2D-Data file
+        :return
+
+        Last edit: February 16th, 2016
+
+        @author: dominiquef
+
+    """
+    from SimPEG import mkvc
+
+    assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'"
+    assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
+
+    fid = open(fileName,'w')
+
+
+    if iptype!=0:
+        fid.write('IPTYPE=%i\n'%iptype)
+
+    else:
+        fid.write('! ' + stype + ' FORMAT\n')
+
+    count = 0
+
+    for ii in range(DCsurvey.nSrc):
+
+        tx = np.c_[DCsurvey.srcList[ii].loc]
+
+        rx = DCsurvey.srcList[ii].rxList[0].locs
+
+        nD = DCsurvey.srcList[ii].nD
+
+        M = rx[0]
+        N = rx[1]
+
+        # Adapt source-receiver location for dtype and stype
+        if dtype=='2D':
+
+            if stype == 'SIMPLE':
+
+                #fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
+                A = np.repeat(tx[0,0],M.shape[0],axis=0)
+                B = np.repeat(tx[0,1],M.shape[0],axis=0)
+                M = M[:,0]
+                N = N[:,0]
+
+                np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
+
+
+            else:
+
+                if stype == 'SURFACE':
+
+                    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='%f',delimiter=' ',newline='\n')
+
+        if dtype=='3D':
+
+            if stype == 'SURFACE':
+
+                fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:]))
+                M = M[:,0:2]
+                N = N[:,0:2]
+
+            if stype == 'GENERAL':
+
+                fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
+
+            fid.write('%i\n'% nD)
+            np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
+            fid.write('\n')
+
+        count += nD
+
+    fid.close()

SimPEG/Examples/DC_Forward_PseudoSection.py

@@ -2,7 +2,7 @@ from SimPEG import Mesh, Utils, np, sp
 import SimPEG.DCIP as DC
 import time
 
-def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', unitType='appConductivity', plotIt=True):
+def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', plotIt=True):
     """
         DC Forward Simulation
         =====================
@@ -15,14 +15,14 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
         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]
-        surveyType = survey type 'pole-dipole' or 'dipole-dipole'
-        unitType = Data type "appResistivity" | "appConductivity"  | "volt"
+        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 surveyType in ['pole-dipole', 'dipole-dipole'], "Source type (surveyType) must be pdp or dpdp (pole dipole or dipole dipole)"
-    assert unitType in ['appResistivity', 'appConductivity', 'volt'], "Unit type (unitType) must be appResistivity or appConductivity or volt (potential)"
+    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.]]
@@ -73,8 +73,8 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
     locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]]
 
     # We will handle the geometry of the survey for you and create all the combination of tx-rx along line
-    # [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])
-    survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2])
+    # [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2])
+    survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2])
 
     # Define some global geometry
     dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )
@@ -118,8 +118,8 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
         rxloc_N = np.asarray(Rx[ii][:,3:])
 
 
-        # For usual cases 'dipole-dipole' or "gradient"
-        if surveyType == 'pole-dipole':
+        # For usual cases "dpdp" or "gradient"
+        if stype == 'pdp':
             # Create an "inifinity" pole
             tx =  np.squeeze(Tx[ii][:,0:1])
             tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
@@ -157,12 +157,12 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
         fig = plt.figure(figsize=(7,7))
         ax = plt.subplot(2,1,1, aspect='equal')
         # Plot the location of the spheres for reference
-        circle1=plt.Circle((loc[0,0], loc[2,0]), radi[0], color='w', fill=False, lw=3)
-        circle2=plt.Circle((loc[0,1], loc[2,1]), radi[1], color='k', fill=False, lw=3)
+        circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
+        circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
         ax.add_artist(circle1)
         ax.add_artist(circle2)
 
-        dat = mesh.plotSlice(np.log10(model), ax = ax, normal = 'Y',
+        dat = mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y',
                              ind = indy,grid=True, clim = np.log10([sig.min(),sig.max()]))
 
         ax.set_title('3-D model')
@@ -188,13 +188,15 @@ def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', u
         ax2 = plt.subplot(2,1,2, aspect='equal')
 
         # Plot the location of the spheres for reference
-        circle1=plt.Circle((loc[0,0], loc[2,0]), radi[0], color='w', fill=False, lw=3)
-        circle2=plt.Circle((loc[0,1], loc[2,1]), radi[1], color='k', fill=False, lw=3)
+        circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
+        circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3)
         ax2.add_artist(circle1)
         ax2.add_artist(circle2)
 
         # Add the speudo section
-        dat = DC.plot_pseudoSection(survey2D, ax2, surveyType=surveyType, unitType=unitType)        # plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
+        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')
         ax2.set_title('Apparent Conductivity data')

SimPEG/Examples/Utils_surface2ind_topo.py

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

SimPEG/Examples/__init__.py

@@ -20,9 +20,8 @@ import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
 import MT_1D_ForwardAndInversion
 import MT_3D_Foward
-import Utils_surface2ind_topo
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Maps.py

@@ -1,4 +1,3 @@
-from __future__ import division
 import Utils, numpy as np, scipy.sparse as sp
 from scipy.sparse.linalg import LinearOperator
 from Tests import checkDerivative
@@ -6,7 +5,6 @@ from PropMaps import PropMap, Property
 from numpy.polynomial import polynomial
 from scipy.interpolate import UnivariateSpline
 import warnings
-from SimPEG.Utils import Zero
 
 class IdentityMap(object):
     """
@@ -19,7 +17,7 @@ class IdentityMap(object):
         Utils.setKwargs(self, **kwargs)
 
         if nP is not None:
-            assert type(nP) in [int, long, np.int64], ' Number of parameters must be an integer.'
+            assert type(nP) in [int, long], ' Number of parameters must be an integer.'
 
         self.mesh = mesh
         self._nP  = nP
@@ -131,15 +129,7 @@ class IdentityMap(object):
 
 
 class ComboMap(IdentityMap):
-    """
-        Combination of various maps.
-
-        The ComboMap holds the information for multiplying and combining
-        maps. It also uses the chain rule to create the derivative.
-        Remember, any time that you make your own combination of mappings
-        be sure to test that the derivative is correct.
-
-    """
+    """Combination of various maps."""
 
     def __init__(self, maps, **kwargs):
         IdentityMap.__init__(self, None, **kwargs)
@@ -188,12 +178,6 @@ class ComboMap(IdentityMap):
 
 class ExpMap(IdentityMap):
     """
-        Electrical conductivity varies over many orders of magnitude, so it is a common
-        technique when solving the inverse problem to parameterize and optimize in terms
-        of log conductivity. This makes sense not only because it ensures all conductivities
-        will be positive, but because this is fundamentally the space where conductivity
-        lives (i.e. it varies logarithmically).
-
         Changes the model into the physical property.
 
         A common example of this is to invert for electrical conductivity
@@ -465,32 +449,6 @@ class Mesh2Mesh(IdentityMap):
     """
         Takes a model on one mesh are translates it to another mesh.
 
-        .. plot::
-
-            from SimPEG import *
-            import matplotlib.pyplot as plt
-            M = Mesh.TensorMesh([100,100])
-            h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
-            h1 = h1/h1.sum()
-            M2 = Mesh.TensorMesh([h1,h1])
-            V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
-            v = Utils.mkvc(V)
-            modh = Maps.Mesh2Mesh([M,M2])
-            modH = Maps.Mesh2Mesh([M2,M])
-            H = modH * v
-            h = modh * H
-            ax = plt.subplot(131)
-            M.plotImage(v, ax=ax)
-            ax.set_title('Fine Mesh (Original)')
-            ax = plt.subplot(132)
-            M2.plotImage(H,clim=[0,1],ax=ax)
-            ax.set_title('Course Mesh')
-            ax = plt.subplot(133)
-            M.plotImage(h,clim=[0,1],ax=ax)
-            ax.set_title('Fine Mesh (Interpolated)')
-            plt.show()
-
-
     """
 
     def __init__(self, meshes, **kwargs):
@@ -543,19 +501,11 @@ class InjectActiveCells(IdentityMap):
         self.indInactive = np.logical_not(indActive)
         if Utils.isScalar(valInactive):
             self.valInactive = np.ones(self.nC)*float(valInactive)
-            self.valInactive[self.indActive] = 0.
         else:
-            if len(valInactive) == sum(self.indInactive):
-                self.valInactive = np.zeros(nC)
-                self.valInactive[self.indInactive] = valInactive.copy()
-            else:
-                assert len(self.valInactive) == self.nC, 'valInactive must be the size of nC or nInactive'
-                self.valInactive = valInactive.copy()
-                if any(self.valInactive[self.indActive] != 0.):
-                    warnings.warn('the inactive has non-zero values in the active set.')
+            self.valInactive = valInactive.copy()
+        self.valInactive[self.indActive] = 0
 
         inds = np.nonzero(self.indActive)[0]
-        # inds[self.indActive]
         self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
 
     @property
@@ -583,6 +533,83 @@ class ActiveCells(InjectActiveCells):
             FutureWarning)
         InjectActiveCells.__init__(self, mesh, indActive, valInactive, nC)
 
+class InjectActiveCellsTopo(IdentityMap):
+    """
+        Active model parameters. Extend for cells on topography to air cell (only works for tensor mesh)
+
+    """
+
+    indActive   = None #: Active Cells
+    valInactive = None #: Values of inactive Cells
+    nC          = None #: Number of cells in the full model
+
+    def __init__(self, mesh, indActive, nC=None):
+        self.mesh  = mesh
+
+        self.nC = nC or mesh.nC
+
+        if indActive.dtype is not bool:
+            z = np.zeros(self.nC,dtype=bool)
+            z[indActive] = True
+            indActive = z
+        self.indActive = indActive
+
+        self.indInactive = np.logical_not(indActive)
+        inds = np.nonzero(self.indActive)[0]
+        self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
+
+    @property
+    def shape(self):
+        return (self.nC, self.nP)
+
+    @property
+    def nP(self):
+        """Number of parameters in the model."""
+        return self.indActive.sum()
+
+    def _transform(self, m):
+        val_temp = np.zeros(self.mesh.nC)
+        val_temp[self.indActive] = m
+        valInactive = np.zeros(self.mesh.nC)
+        #1D
+        if self.mesh.dim == 1:
+            z_temp = self.mesh.gridCC
+            val_temp[~self.indActive] = val_temp[np.argmax(z_temp[self.indActive])]
+        #2D
+        elif self.mesh.dim == 2:
+            act_temp = self.indActive.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F')
+            val_temp = val_temp.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F')
+            y_temp = self.mesh.gridCC[:,1].reshape((self.mesh.nCx, self.mesh.nCy), order = 'F')
+            for i in range(self.mesh.nCx):
+                act_tempx = act_temp[i,:] == 1
+                val_temp[i,~act_tempx] = val_temp[i,np.argmax(y_temp[i,act_tempx])]
+            valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive]
+        #3D
+        elif self.mesh.dim == 3:
+            act_temp = self.indActive.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
+            val_temp = val_temp.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
+            z_temp = self.mesh.gridCC[:,2].reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F')
+            for i in range(self.mesh.nCx*self.mesh.nCy):
+                act_tempxy = act_temp[i,:] == 1
+                val_temp[i,~act_tempxy] = val_temp[i,np.argmax(z_temp[i,act_tempxy])]
+            valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive]
+
+        self.valInactive = valInactive
+
+        return self.P*m + self.valInactive
+
+    def inverse(self, D):
+        return self.P.T*D
+
+    def deriv(self, m):
+        return self.P
+
+class ActiveCellsTopo(InjectActiveCellsTopo):
+    def __init__(self, mesh, indActive, valInactive, nC=None):
+        warnings.warn(
+            "`ActiveCellsTopo` is deprecated and will be removed in future versions. Use `InjectActiveCellsTopo` instead",
+            FutureWarning)
+        InjectActiveCellsTopo.__init__(self, mesh, indActive, valInactive, nC)
 
 class Weighting(IdentityMap):
     """
@@ -624,37 +651,6 @@ class Weighting(IdentityMap):
     def deriv(self, m):
         return self.P
 
-class Projection(IdentityMap):
-    """
-        A map to rearrange parameters
-
-    """
-
-
-    def __init__(self, indTo, indFrom, shape, mesh=None, **kwargs):
-
-        assert len(indTo) == len(indFrom)
-
-        self.P = sp.csr_matrix((np.ones(len(indTo)), (indTo, indFrom)), shape=shape)
-        self._shape = shape
-
-        super(Projection, self).__init__(mesh, **kwargs)
-
-    @property
-    def shape(self):
-        return self._shape
-
-    @property
-    def nP(self):
-        """Number of parameters in the model."""
-        return self.shape[1]
-
-    def _transform(self, m):
-        return self.P*m
-
-    def deriv(self, m):
-        return self.P
-
 
 class ComplexMap(IdentityMap):
     """ComplexMap
@@ -697,13 +693,13 @@ class CircleMap(IdentityMap):
 
         Parameterize the model space using a circle in a wholespace.
 
-        .. math::
+        ..math::
 
             \sigma(m) = \sigma_1 + (\sigma_2 - \sigma_1)\left(\\arctan\left(100*\sqrt{(\\vec{x}-x_0)^2 + (\\vec{y}-y_0)}-r\\right) \pi^{-1} + 0.5\\right)
 
         Define the model as:
 
-        .. math::
+        ..math::
 
             m = [\sigma_1, \sigma_2, x_0, y_0, r]
 
@@ -1056,697 +1052,7 @@ class SplineMap(IdentityMap):
         return sp.csr_matrix(np.c_[g1,g2,g3])
 
 
-class ParametrizedLayer(IdentityMap):
-    """
-        Parametrized Layer Space
 
-        m = [val_background, val_layer, layer_center, layer_thickness]
 
 
-        .. plot::
-            :include-source:
-
-            from SimPEG import Mesh, Maps, np
-            import matplotlib.pyplot as plt
-
-            fig, ax = plt.subplots(1,1,figsize=(2,3))
-
-            mesh = Mesh.TensorMesh([50,50],x0='CC')
-            mapping = Maps.ParametrizedLayer(mesh)
-            m = np.hstack(np.r_[1., 2., -0.1, 0.2])
-            rho = mapping._transform(m)
-            mesh.plotImage(rho, ax=ax)
-
-        **Required**
-
-        :param Mesh mesh: SimPEG Mesh, 2D or 3D
-
-        **Optional**
-
-        :param float slopeFact: arctan slope factor - divided by the minimum h spacing to give the slope of the arctan functions
-        :param float slope: slope of the arctan function
-        :param numpy.ndarray indActive: bool vector with
-
-    """
-
-    slopeFact = 1e2 # will be scaled by the mesh.
-    slope = None
-    indActive = None
-
-    def __init__(self, mesh, **kwargs):
-
-        super(ParametrizedLayer, self).__init__(mesh, **kwargs)
-
-
-        if self.slope is None:
-            self.slope = self.slopeFact / np.hstack(self.mesh.h).min()
-
-        self.x = [self.mesh.gridCC[:,0] if self.indActive is None else self.mesh.gridCC[self.indActive,0]][0]
-
-        if self.mesh.dim > 1:
-            self.y = [self.mesh.gridCC[:,1] if self.indActive is None else self.mesh.gridCC[self.indActive,1]][0]
-
-        if self.mesh.dim > 2:
-            self.z = [self.mesh.gridCC[:,2] if self.indActive is None else self.mesh.gridCC[self.indActive,2]][0]
-
-    @property
-    def nP(self):
-        return 4
-
-    @property
-    def shape(self):
-        if self.indActive is not None:
-            return (sum(self.indActive), self.nP)
-        return (self.mesh.nC, self.nP)
-
-    def mDict(self, m):
-        return {
-                'val_background': m[0],
-                'val_layer': m[1],
-                'layer_center': m[2],
-                'layer_thickness': m[3],
-                }
-
-    def _atanfct(self, xyz, xyzi, slope):
-        return np.arctan(slope * (xyz - xyzi))/np.pi + 0.5
-
-    def _atanfctDeriv(self, xyz, xyzi, slope):
-        # d/dx(atan(x)) = 1/(1+x**2)
-        x = slope * (xyz - xyzi)
-        dx = - slope
-        return (1./(1 + x**2))/np.pi * dx
-
-    def _atanLayer(self, mDict):
-        if self.mesh.dim == 2:
-            z = self.y
-        elif self.mesh.dim == 3:
-            z = self.z
-
-        layer_bottom = mDict['layer_center'] - mDict['layer_thickness'] / 2.
-        layer_top = mDict['layer_center'] + mDict['layer_thickness'] / 2.
-        return self._atanfct(z, layer_bottom, self.slope)*self._atanfct(z, layer_top, -self.slope)
-
-    def _atanLayerDeriv_layer_center(self, mDict):
-        if self.mesh.dim == 2:
-            z = self.y
-        elif self.mesh.dim == 3:
-            z = self.z
-
-        layer_bottom = mDict['layer_center'] - mDict['layer_thickness'] / 2.
-        layer_top = mDict['layer_center'] + mDict['layer_thickness'] / 2.
-
-        return (self._atanfctDeriv(z, layer_bottom, self.slope)*self._atanfct(z, layer_top, -self.slope)
-                + self._atanfct(z, layer_bottom, self.slope)*self._atanfctDeriv(z, layer_top, -self.slope))
-
-    def _atanLayerDeriv_layer_thickness(self, mDict):
-        if self.mesh.dim == 2:
-            z = self.y
-        elif self.mesh.dim == 3:
-            z = self.z
-
-        layer_bottom = mDict['layer_center'] - mDict['layer_thickness'] / 2.
-        layer_top = mDict['layer_center'] + mDict['layer_thickness'] / 2.
-
-        return (-0.5*self._atanfctDeriv(z, layer_bottom, self.slope)*self._atanfct(z, layer_top, -self.slope)
-                + 0.5*self._atanfct(z, layer_bottom, self.slope)*self._atanfctDeriv(z, layer_top, -self.slope))
-
-    def layer_cont(self, mDict):
-        return mDict['val_background'] + (mDict['val_layer'] - mDict['val_background'])*self._atanLayer(mDict)
-
-    def _transform(self, m):
-        mDict = self.mDict(m)
-        return self.layer_cont(mDict)
-
-    def _deriv_val_background(self, mDict):
-        return np.ones_like(self.x) - self._atanLayer(mDict)
-
-    def _deriv_val_layer(self, mDict):
-        return self._atanLayer(mDict)
-
-    def _deriv_layer_center(self, mDict):
-        return (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_center(mDict)
-
-    def _deriv_layer_thickness(self, mDict):
-        return (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_thickness(mDict)
-
-    def deriv(self, m):
-
-        mDict = self.mDict(m)
-
-        return sp.csr_matrix(np.vstack([
-            self._deriv_val_background(mDict),
-            self._deriv_val_layer(mDict),
-            self._deriv_layer_center(mDict),
-            self._deriv_layer_thickness(mDict),
-            ]).T)
-
-
-class ParametrizedCasingAndLayer(ParametrizedLayer):
-    """
-        Parametrized layered space with casing.
-
-        m = [val_background, val_layer, val_casing, val_insideCasing, layer_center, layer_thickness, casing_radius, casing_thickness, casing_bottom, casing_top]
-    """
-
-    def __init__(self, mesh, **kwargs):
-
-        assert mesh._meshType == 'CYL', 'Parametrized Casing in a layer map only works for a cyl mesh.'
-
-        super(ParametrizedCasingAndLayer, self).__init__(mesh, **kwargs)
-
-
-    @property
-    def nP(self):
-        return 10
-
-    @property
-    def shape(self):
-        if self.indActive is not None:
-            return (sum(self.indActive), self.nP)
-        return (self.mesh.nC, self.nP)
-
-    def mDict(self, m):
-        #m = [val_background, val_layer, val_casing, val_insideCasing, layer_center, layer_thickness, casing_radius, casing_thickness, casing_bottom, casing_top]
-        return {
-                    'val_background': m[0],
-                    'val_layer': m[1],
-                    'val_casing': m[2],
-                    'val_insideCasing': m[3],
-                    'layer_center': m[4],
-                    'layer_thickness': m[5],
-                    'casing_radius': m[6],
-                    'casing_thickness': m[7],
-                    'casing_bottom': m[8],
-                    'casing_top': m[9]
-               }
-
-    def _atanCasingLength(self, mDict):
-        return (self._atanfct(self.z, mDict['casing_top'], -self.slope)
-                * self._atanfct(self.z, mDict['casing_bottom'], self.slope))
-
-    def _atanCasingLengthDeriv_casing_top(self, mDict):
-        return (self._atanfctDeriv(self.z, mDict['casing_top'], -self.slope)
-                * self._atanfct(self.z, mDict['casing_bottom'], self.slope))
-
-    def _atanCasingLengthDeriv_casing_bottom(self, mDict):
-        return (self._atanfct(self.z, mDict['casing_top'], -self.slope)
-                * self._atanfctDeriv(self.z, mDict['casing_bottom'], self.slope))
-
-    def _atanInsideCasing(self, mDict):
-        casing_a = mDict['casing_radius'] - 0.5*mDict['casing_thickness']
-        return (self._atanCasingLength(mDict)
-                * self._atanfct(self.x, casing_a, -self.slope))
-
-    def _atanInsideCasingDeriv_casing_radius(self, mDict):
-        casing_a = mDict['casing_radius'] - 0.5*mDict['casing_thickness']
-        return (self._atanCasingLength(mDict)
-                * self._atanfctDeriv(self.x, casing_a, -self.slope))
-
-    def _atanInsideCasingDeriv_casing_thickness(self, mDict):
-        casing_a = mDict['casing_radius'] - 0.5*mDict['casing_thickness']
-        return (self._atanCasingLength(mDict)
-                * - 0.5*self._atanfctDeriv(self.x, casing_a, -self.slope))
-
-    def _atanInsideCasingDeriv_casing_top(self, mDict):
-        casing_a = mDict['casing_radius'] - 0.5*mDict['casing_thickness']
-        return (self._atanCasingLengthDeriv_casing_top(mDict)
-                * self._atanfct(self.x, casing_a, -self.slope))
-
-    def _atanInsideCasingDeriv_casing_bottom(self, mDict):
-        casing_a = mDict['casing_radius'] - 0.5*mDict['casing_thickness']
-        return (self._atanCasingLengthDeriv_casing_bottom(mDict)
-                * self._atanfct(self.x, casing_a, -self.slope))
-
-    def _atanCasing(self, mDict):
-        casing_a, casing_b = mDict['casing_radius'] - 0.5*mDict['casing_thickness'], mDict['casing_radius'] + 0.5*mDict['casing_thickness']
-        return (self._atanCasingLength(mDict)
-                * self._atanfct(self.x, casing_a, self.slope)
-                * self._atanfct(self.x, casing_b, -self.slope))
-
-    def _atanCasingDeriv_casing_radius(self, mDict):
-        casing_a, casing_b = mDict['casing_radius'] - 0.5*mDict['casing_thickness'], mDict['casing_radius'] + 0.5*mDict['casing_thickness']
-        return (self._atanCasingLength(mDict) * (
-                self._atanfctDeriv(self.x, casing_a, self.slope)
-                * self._atanfct(self.x, casing_b, -self.slope)
-                +
-                self._atanfct(self.x, casing_a, self.slope)
-                * self._atanfctDeriv(self.x, casing_b, -self.slope)
-                ))
-
-    def _atanCasingDeriv_casing_thickness(self, mDict):
-        casing_a, casing_b = mDict['casing_radius'] - 0.5*mDict['casing_thickness'], mDict['casing_radius'] + 0.5*mDict['casing_thickness']
-        return (self._atanCasingLength(mDict) * (
-                - 0.5*self._atanfctDeriv(self.x, casing_a, self.slope)
-                * 0.5*self._atanfct(self.x, casing_b, -self.slope)
-                +
-                - 0.5*self._atanfct(self.x, casing_a, self.slope)
-                * 0.5*self._atanfctDeriv(self.x, casing_b, -self.slope)
-                ))
-
-    def _atanCasingDeriv_casing_bottom(self, mDict):
-        casing_a, casing_b = mDict['casing_radius'] - 0.5*mDict['casing_thickness'], mDict['casing_radius'] + 0.5*mDict['casing_thickness']
-        return (self._atanCasingLengthDeriv_casing_bottom(mDict)
-                * self._atanfct(self.x, casing_a, self.slope)
-                * self._atanfct(self.x, casing_b, -self.slope))
-
-    def _atanCasingDeriv_casing_top(self, mDict):
-        casing_a, casing_b = mDict['casing_radius'] - 0.5*mDict['casing_thickness'], mDict['casing_radius'] + 0.5*mDict['casing_thickness']
-        return (self._atanCasingLengthDeriv_casing_top(mDict)
-                * self._atanfct(self.x, casing_a, self.slope)
-                * self._atanfct(self.x, casing_b, -self.slope))
-
-    def layer_cont(self, mDict):
-        return mDict['val_background'] + (mDict['val_layer']-mDict['val_background']) * self._atanLayer(mDict) # contribution from the layered background
-
-
-    def _transform(self, m):
-
-        mDict = self.mDict(m)
-
-        # assemble the model
-        layer = self.layer_cont(mDict)
-        casing = (mDict['val_casing'] - layer) * self._atanCasing(mDict)
-        insideCasing = (mDict['val_insideCasing'] - layer) * self._atanInsideCasing(mDict)
-
-        return layer + casing + insideCasing
-
-
-    def _deriv_val_background(self, mDict):
-        d_layer_cont_dval_background = 1. - self._atanLayer(mDict) # contribution from the layered background
-        d_casing_cont_dval_background = -1. * d_layer_cont_dval_background * self._atanCasing(mDict)
-        d_insideCasing_cont_dval_background = -1. * d_layer_cont_dval_background * self._atanInsideCasing(mDict)
-        return d_layer_cont_dval_background + d_casing_cont_dval_background + d_insideCasing_cont_dval_background
-
-    def _deriv_val_layer(self, mDict):
-        d_layer_cont_dval_layer = self._atanLayer(mDict)
-        d_casing_cont_dval_layer = -1. * d_layer_cont_dval_layer * self._atanCasing(mDict)
-        d_insideCasing_cont_dval_layer = -1. * d_layer_cont_dval_layer * self._atanInsideCasing(mDict)
-        return d_layer_cont_dval_layer + d_casing_cont_dval_layer + d_insideCasing_cont_dval_layer
-
-    def _deriv_val_casing(self, mDict):
-        d_layer_cont_dval_casing = 0.
-        d_casing_cont_dval_casing = self._atanCasing(mDict)
-        d_insideCasing_cont_dval_casing = 0.
-        return d_layer_cont_dval_casing + d_casing_cont_dval_casing + d_insideCasing_cont_dval_casing
-
-    def _deriv_val_insideCasing(self, mDict):
-        d_layer_cont_dval_insideCasing = 0.
-        d_casing_cont_dval_insideCasing = 0.
-        d_insideCasing_cont_dval_insideCasing = self._atanInsideCasing(mDict)
-        return d_layer_cont_dval_insideCasing + d_casing_cont_dval_insideCasing + d_insideCasing_cont_dval_insideCasing
-
-    def _deriv_layer_center(self, mDict):
-        d_layer_cont_dlayer_center = (mDict['val_layer'] - mDict['val_background']) * self._atanLayerDeriv_layer_center(mDict)
-        d_casing_cont_dlayer_center =  - d_layer_cont_dlayer_center * self._atanCasing(mDict)
-        d_insideCasing_cont_dlayer_center = - d_layer_cont_dlayer_center * self._atanInsideCasing(mDict)
-        return d_layer_cont_dlayer_center + d_casing_cont_dlayer_center + d_insideCasing_cont_dlayer_center
-
-    def _deriv_layer_thickness(self, mDict):
-        d_layer_cont_dlayer_thickness = (mDict['val_layer']-mDict['val_background']) * self._atanLayerDeriv_layer_thickness(mDict)
-        d_casing_cont_dlayer_thickness =  - d_layer_cont_dlayer_thickness * self._atanCasing(mDict)
-        d_insideCasing_cont_dlayer_thickness = - d_layer_cont_dlayer_thickness * self._atanInsideCasing(mDict)
-        return d_layer_cont_dlayer_thickness + d_casing_cont_dlayer_thickness + d_insideCasing_cont_dlayer_thickness
-
-    def _deriv_casing_radius(self, mDict):
-        layer = self.layer_cont(mDict)
-        d_layer_cont_dcasing_radius = 0.
-        d_casing_cont_dcasing_radius = (mDict['val_casing'] - layer) * self._atanCasingDeriv_casing_radius(mDict)
-        d_insideCasing_cont_dcasing_radius = (mDict['val_insideCasing'] - layer) * self._atanInsideCasingDeriv_casing_radius(mDict)
-        return d_layer_cont_dcasing_radius + d_casing_cont_dcasing_radius + d_insideCasing_cont_dcasing_radius
-
-    def _deriv_casing_thickness(self, mDict):
-        d_layer_cont_dcasing_thickness = 0.
-        d_casing_cont_dcasing_thickness = (mDict['val_casing'] - self.layer_cont(mDict)) * self._atanCasingDeriv_casing_thickness(mDict)
-        d_insideCasing_cont_dcasing_thickness = (mDict['val_insideCasing'] - self.layer_cont(mDict)) * self._atanInsideCasingDeriv_casing_thickness(mDict)
-        return d_layer_cont_dcasing_thickness + d_casing_cont_dcasing_thickness + d_insideCasing_cont_dcasing_thickness
-
-    def _deriv_casing_bottom(self, mDict):
-        d_layer_cont_dcasing_bottom = 0.
-        d_casing_cont_dcasing_bottom = (mDict['val_casing'] - self.layer_cont(mDict)) * self._atanCasingDeriv_casing_bottom(mDict)
-        d_insideCasing_cont_dcasing_bottom = (mDict['val_insideCasing'] - self.layer_cont(mDict)) * self._atanInsideCasingDeriv_casing_bottom(mDict)
-        return d_layer_cont_dcasing_bottom + d_casing_cont_dcasing_bottom + d_insideCasing_cont_dcasing_bottom
-
-    def _deriv_casing_top(self, mDict):
-        d_layer_cont_dcasing_top = 0.
-        d_casing_cont_dcasing_top = (mDict['val_casing'] - self.layer_cont(mDict)) * self._atanCasingDeriv_casing_top(mDict)
-        d_insideCasing_cont_dcasing_top = (mDict['val_insideCasing'] - self.layer_cont(mDict)) * self._atanInsideCasingDeriv_casing_top(mDict)
-        return d_layer_cont_dcasing_top + d_casing_cont_dcasing_top + d_insideCasing_cont_dcasing_top
-
-
-    def deriv(self, m):
-
-        mDict = self.mDict(m)
-
-        return sp.csr_matrix(np.vstack([
-                self._deriv_val_background(mDict),
-                self._deriv_val_layer(mDict),
-                self._deriv_val_casing(mDict),
-                self._deriv_val_insideCasing(mDict),
-                self._deriv_layer_center(mDict),
-                self._deriv_layer_thickness(mDict),
-                self._deriv_casing_radius(mDict),
-                self._deriv_casing_thickness(mDict),
-                self._deriv_casing_bottom(mDict),
-                self._deriv_casing_top(mDict),
-                ]).T)
-
-
-
-class ParametrizedBlockInLayer(ParametrizedLayer):
-    """
-        Parametrized Block in a Layered Space
-
-        For 2D:
-        m = [val_background, val_layer, val_block, layer_center, layer_thickness, block_x0, block_dx]
-
-        For 3D:
-        m = [val_background, val_layer, val_block, layer_center, layer_thickness, block_x0, block_y0, block_dx, block_dy]
-
-        .. plot::
-            :include-source:
-
-            from SimPEG import Mesh, Maps, np
-            import matplotlib.pyplot as plt
-
-            fig, ax = plt.subplots(1,1,figsize=(2,3))
-
-            mesh = Mesh.TensorMesh([50,50],x0='CC')
-            mapping = Maps.ParametrizedBlockInLayer(mesh)
-            m = np.hstack(np.r_[1., 2., 3., -0.1, 0.2, 0.3, 0.2])
-            rho = mapping._transform(m)
-            mesh.plotImage(rho, ax=ax)
-
-        **Required**
-
-        :param Mesh mesh: SimPEG Mesh, 2D or 3D
-
-        **Optional**
-
-        :param float slopeFact: arctan slope factor - divided by the minimum h spacing to give the slope of the arctan functions
-        :param float slope: slope of the arctan function
-        :param numpy.ndarray indActive: bool vector with
-
-    """
-
-    def __init__(self, mesh, **kwargs):
-
-        super(ParametrizedBlockInLayer, self).__init__(mesh, **kwargs)
-
-    @property
-    def nP(self):
-        if self.mesh.dim == 2:
-            return 7
-        elif self.mesh.dim == 3:
-            return 9
-
-    @property
-    def shape(self):
-        if self.indActive is not None:
-            return (sum(self.indActive), self.nP)
-        return (self.mesh.nC, self.nP)
-
-    def _mDict2d(self, m):
-        return{
-                'val_background': m[0],
-                'val_layer': m[1],
-                'val_block': m[2],
-                'layer_center': m[3],
-                'layer_thickness': m[4],
-                'x0_block': m[5],
-                'dx_block': m[6]
-              }
-
-    def _mDict3d(self, m):
-        return{
-                'val_background': m[0],
-                'val_layer': m[1],
-                'val_block': m[2],
-                'layer_center': m[3],
-                'layer_thickness': m[4],
-                'x0_block': m[5],
-                'y0_block': m[6],
-                'dx_block': m[7],
-                'dy_block': m[8]
-              }
-
-    def mDict(self, m):
-        if self.mesh.dim == 2:
-            return self._mDict2d(m)
-        elif self.mesh.dim == 3:
-            return self._mDict3d(m)
-
-    def _atanBlock2d(self, mDict):
-        return (self._atanLayer(mDict)
-                          * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                          * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
-
-    def _atanBlock2dDeriv_layer_center(self, mDict):
-        return (self._atanLayerDeriv_layer_center(mDict)
-                          * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                          * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
-
-    def _atanBlock2dDeriv_layer_thickness(self, mDict):
-        return (self._atanLayerDeriv_layer_thickness(mDict)
-                          * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                          * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
-
-
-    def _atanBlock2dDeriv_x0(self, mDict):
-        return self._atanLayer(mDict) * (
-                    (self._atanfctDeriv(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
-                    +
-                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                    * self._atanfctDeriv(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
-                    )
-
-    def _atanBlock2dDeriv_dx(self, mDict):
-        return self._atanLayer(mDict) * (
-                    (self._atanfctDeriv(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope) * -0.5
-                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope))
-                    +
-                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                    * self._atanfctDeriv(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope) * 0.5)
-                    )
-
-    def _atanBlock3d(self, mDict):
-        return (self._atanLayer(mDict)
-                      * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                      * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
-                      * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
-                      * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
-
-
-    def _atanBlock3dDeriv_layer_center(self, mDict):
-        return (self._atanLayerDeriv_layer_center(mDict)
-                      * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                      * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
-                      * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
-                      * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
-
-    def _atanBlock3dDeriv_layer_thickness(self, mDict):
-        return (self._atanLayerDeriv_layer_thickness(mDict)
-                      * self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                      * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
-                      * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
-                      * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
-
-
-    def _atanBlock3dDeriv_x0(self, mDict):
-        return self._atanLayer(mDict) * (
-                    (self._atanfctDeriv(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
-                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
-                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
-                    +
-                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                    * self._atanfctDeriv(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
-                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
-                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
-                    )
-
-    def _atanBlock3dDeriv_y0(self, mDict):
-        return self._atanLayer(mDict) * (
-                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
-                    * self._atanfctDeriv(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
-                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
-                    +
-                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
-                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
-                    * self._atanfctDeriv(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
-                    )
-
-    def _atanBlock3dDeriv_dx(self, mDict):
-        return self._atanLayer(mDict) * (
-                    (self._atanfctDeriv(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope) * -0.5
-                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
-                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
-                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
-                    +
-                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                    * self._atanfctDeriv(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope) * 0.5
-                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
-                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
-                    )
-
-    def _atanBlock3dDeriv_dy(self, mDict):
-        return self._atanLayer(mDict) * (
-                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
-                    * self._atanfctDeriv(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope) * -0.5
-                    * self._atanfct(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope))
-                    +
-                    (self._atanfct(self.x, mDict['x0_block'] - 0.5*mDict['dx_block'], self.slope)
-                    * self._atanfct(self.x, mDict['x0_block'] + 0.5*mDict['dx_block'], -self.slope)
-                    * self._atanfct(self.y, mDict['y0_block'] - 0.5*mDict['dy_block'], self.slope)
-                    * self._atanfctDeriv(self.y, mDict['y0_block'] + 0.5*mDict['dy_block'], -self.slope) * 0.5)
-                    )
-
-
-    def _transform2d(self, m):
-        mDict = self.mDict(m)
-        # assemble the model
-        layer_cont = mDict['val_background'] + (mDict['val_layer']-mDict['val_background'])*self._atanLayer(mDict) # contribution from the layered background
-        block_cont = (mDict['val_block']-layer_cont)*self._atanBlock2d(mDict) # perturbation due to the block
-
-        return layer_cont + block_cont
-
-    def _deriv2d_val_background(self, mDict):
-        d_layer_dval_background = np.ones_like(self.x) - self._atanLayer(mDict)
-        d_block_dval_background = (-d_layer_dval_background)*self._atanBlock2d(mDict)
-        return d_layer_dval_background + d_block_dval_background
-
-    def _deriv2d_val_layer(self, mDict):
-        d_layer_dval_layer = self._atanLayer(mDict)
-        d_block_dval_layer = (-d_layer_dval_layer)*self._atanBlock2d(mDict)
-        return d_layer_dval_layer + d_block_dval_layer
-
-    def _deriv2d_val_block(self, mDict):
-        d_layer_dval_block = 0.
-        d_block_dval_block = (1.-d_layer_dval_block)*self._atanBlock2d(mDict)
-        return d_layer_dval_block + d_block_dval_block
-
-    def _deriv2d_layer_center(self, mDict):
-        d_layer_dlayer_center = (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_center(mDict)
-        d_block_dlayer_center = ((mDict['val_block']-self.layer_cont(mDict))*self._atanBlock2dDeriv_layer_center(mDict)
-                                    - d_layer_dlayer_center*self._atanBlock2d(mDict))
-        return d_layer_dlayer_center + d_block_dlayer_center
-
-    def _deriv2d_layer_thickness(self, mDict):
-        d_layer_dlayer_thickness = (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_thickness(mDict)
-        d_block_dlayer_thickness = ((mDict['val_block']-self.layer_cont(mDict))*self._atanBlock2dDeriv_layer_thickness(mDict)
-                                    - d_layer_dlayer_thickness*self._atanBlock2d(mDict))
-        return d_layer_dlayer_thickness + d_block_dlayer_thickness
-
-    def _deriv2d_x0_block(self, mDict):
-        d_layer_dx0 = 0.
-        d_block_dx0 = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock2dDeriv_x0(mDict)
-        return d_layer_dx0 + d_block_dx0
-
-    def _deriv2d_dx_block(self, mDict):
-        d_layer_ddx = 0.
-        d_block_ddx = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock2dDeriv_dx(mDict)
-        return d_layer_ddx + d_block_ddx
-
-    def _deriv2d(self, m):
-        mDict = self.mDict(m)
-
-        return np.vstack([
-                            self._deriv2d_val_background(mDict),
-                            self._deriv2d_val_layer(mDict),
-                            self._deriv2d_val_block(mDict),
-                            self._deriv2d_layer_center(mDict),
-                            self._deriv2d_layer_thickness(mDict),
-                            self._deriv2d_x0_block(mDict),
-                            self._deriv2d_dx_block(mDict)
-                        ]).T
-
-    def _transform3d(self, m):
-        # parse model
-        mDict = self.mDict(m)
-
-        # assemble the model
-        layer_cont = mDict['val_background'] + (mDict['val_layer']-mDict['val_background'])*self._atanLayer(mDict) # contribution from the layered background
-        block_cont = (mDict['val_block']-layer_cont)*self._atanBlock3d(mDict) # perturbation due to the block
-
-        return layer_cont + block_cont
-
-    def _deriv3d_val_background(self, mDict):
-        d_layer_dval_background = np.ones_like(self.x) - self._atanLayer(mDict)
-        d_block_dval_background = (-d_layer_dval_background)*self._atanBlock3d(mDict)
-        return d_layer_dval_background + d_block_dval_background
-
-    def _deriv3d_val_layer(self, mDict):
-        d_layer_dval_layer = self._atanLayer(mDict)
-        d_block_dval_layer = (-d_layer_dval_layer)*self._atanBlock3d(mDict)
-        return d_layer_dval_layer + d_block_dval_layer
-
-    def _deriv3d_val_block(self, mDict):
-        d_layer_dval_block = 0.
-        d_block_dval_block = (1.-d_layer_dval_block)*self._atanBlock3d(mDict)
-        return d_layer_dval_block + d_block_dval_block
-
-    def _deriv3d_layer_center(self, mDict):
-        d_layer_dlayer_center = (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_center(mDict)
-        d_block_dlayer_center = ((mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_layer_center(mDict)
-                                    - d_layer_dlayer_center*self._atanBlock3d(mDict))
-        return d_layer_dlayer_center + d_block_dlayer_center
-
-    def _deriv3d_layer_thickness(self, mDict):
-        d_layer_dlayer_thickness = (mDict['val_layer']-mDict['val_background'])*self._atanLayerDeriv_layer_thickness(mDict)
-        d_block_dlayer_thickness = ((mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_layer_thickness(mDict)
-                                - d_layer_dlayer_thickness*self._atanBlock3d(mDict))
-        return d_layer_dlayer_thickness + d_block_dlayer_thickness
-
-    def _deriv3d_x0_block(self, mDict):
-        d_layer_dx0 = 0.
-        d_block_dx0 = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_x0(mDict)
-        return d_layer_dx0 + d_block_dx0
-
-    def _deriv3d_y0_block(self, mDict):
-        d_layer_dy0 = 0.
-        d_block_dy0 = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_y0(mDict)
-        return d_layer_dy0 + d_block_dy0
-
-    def _deriv3d_dx_block(self, mDict):
-        d_layer_ddx = 0.
-        d_block_ddx = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_dx(mDict)
-        return d_layer_ddx + d_block_ddx
-
-    def _deriv3d_dy_block(self, mDict):
-        d_layer_ddy = 0.
-        d_block_ddy = (mDict['val_block']-self.layer_cont(mDict))*self._atanBlock3dDeriv_dy(mDict)
-        return d_layer_ddy + d_block_ddy
-
-    def _deriv3d(self, m):
-
-        mDict = self.mDict(m)
-
-        return np.vstack([
-                            self._deriv3d_val_background(mDict),
-                            self._deriv3d_val_layer(mDict),
-                            self._deriv3d_val_block(mDict),
-                            self._deriv3d_layer_center(mDict),
-                            self._deriv3d_layer_thickness(mDict),
-                            self._deriv3d_x0_block(mDict),
-                            self._deriv3d_y0_block(mDict),
-                            self._deriv3d_dx_block(mDict),
-                            self._deriv3d_dy_block(mDict),
-                        ]).T
-
-    def _transform(self, m):
-
-        if self.mesh.dim == 2:
-            return self._transform2d(m)
-        elif self.mesh.dim == 3:
-            return self._transform3d(m)
-
-    def deriv(self, m):
-
-        if self.mesh.dim == 2:
-            return sp.csr_matrix(self._deriv2d(m))
-        elif self.mesh.dim == 3:
-            return sp.csr_matrix(self._deriv3d(m))
 

SimPEG/Mesh/MeshIO.py

@@ -205,19 +205,30 @@ class TensorMeshIO(object):
             :param simpeg.Mesh.TensorMesh mesh: The mesh
 
         """
-        assert mesh.dim == 3
-        s = ''
-        s += '%i %i %i\n' %tuple(mesh.vnC)
-        origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated.
-        origin.dtype = float
+        if mesh.dim ==3:
+            s = ''
+            s += '%i %i %i\n' %tuple(mesh.vnC)
+            origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated.
+            origin.dtype = float
 
-        s += '%.2f %.2f %.2f\n' %tuple(origin)
-        s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
-        s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
-        s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1])
-        f = open(fileName, 'w')
-        f.write(s)
-        f.close()
+            s += '%.2f %.2f %.2f\n' %tuple(origin)
+            s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx)
+            s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy)
+            s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1])
+            f = open(fileName, 'w')
+            f.write(s)
+            f.close()
+
+        elif mesh.dim==2:
+            fid = open(fileName,'w')
+            fid.write('%i\n'% mesh.nCx)
+            fid.write('%f %f 1\n'% (mesh.vectorNx[0],mesh.vectorNx[1]))
+            np.savetxt(fid, np.c_[mesh.vectorNx[2:],np.ones(mesh.nCx-1)], fmt='\t %e %i',delimiter=' ',newline='\n')
+            fid.write('\n')
+            fid.write('%i\n'% mesh.nCy)
+            fid.write('%f %f 1\n'%( 0,mesh.hy[-1]))
+            np.savetxt(fid, np.c_[np.cumsum(mesh.hy[-2::-1])+mesh.hy[-1],np.ones(mesh.nCy-1)], fmt='\t %e %i',delimiter=' ',newline='\n')
+            fid.close()
 
         if models is None: return
         assert type(models) is dict, 'models must be a dict'

SimPEG/PropMaps.py

@@ -74,7 +74,7 @@ class Property(object):
                 if linkedMap is None:
                     return None
                 linkMap = linkMapClass(None) * linkedMap
-                m = getattr(self, '%sModel'%linkName)
+                m = getattr(self, '%s'%linkName)
                 return linkMap.deriv( m )
 
             m = getattr(self, '%sModel'%prop.name)
@@ -239,7 +239,7 @@ class PropMap(object):
             setattr(self, '%sMap'%name, mapping)
             setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP)))
             nP += mapping.nP
-        self.nP = nP
+        self.nP = nP 
 
     @property
     def defaultInvProp(self):

SimPEG/Regularization.py

@@ -39,7 +39,7 @@ class RegularizationMesh(object):
             if self.indActive is None:
                 self._nC = self.mesh.nC
             else:
-                self._nC = int(sum(self.indActive))
+                self._nC = sum(self.indActive)
         return self._nC
 
     @property
@@ -304,7 +304,7 @@ class BaseRegularization(object):
     mesh    = None    #: A SimPEG.Mesh instance.
     mref    = None    #: Reference model.
 
-    def __init__(self, mesh=None, nP=None, mapping=None, indActive=None, **kwargs):
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
         assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
         if indActive is not None and indActive.dtype != 'bool':
@@ -314,18 +314,10 @@ class BaseRegularization(object):
         if indActive is not None and mapping is None:
             mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
 
-        if mesh is None and nP is None:
-            raise Exception, 'either Mesh or number of parameters must be provided to the BaseRegularization'
-
         self.regmesh = RegularizationMesh(mesh,indActive)
-        self.indActive = indActive
-
-        if mesh is not None and nP is None:
-            nP = self.regmesh.nC
-        self.nP = nP
-
-        self.mapping = mapping or self.mapPair(nP=self.nP)
+        self.mapping = mapping or self.mapPair(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
+        self.indActive = indActive
 
     @property
     def parent(self):
@@ -354,7 +346,7 @@ class BaseRegularization(object):
     @property
     def W(self):
         """Full regularization weighting matrix W."""
-        return sp.identity(self.nP)
+        return sp.identity(self.regmesh.nC)
 
     @Utils.timeIt
     def eval(self, m):

SimPEG/Utils/__init__.py

@@ -7,4 +7,3 @@ from CounterUtils import *
 import ModelBuilder
 import SolverUtils
 from coordutils import *
-from modelutils import *

SimPEG/Utils/modelutils.py

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

docs/api_Maps.rst

@@ -122,10 +122,92 @@ When these are used in the inverse problem, this is extremely important!!
 The API
 =======
 
-.. automodule:: SimPEG.Maps
+.. autoclass:: SimPEG.Maps.IdentityMap
     :members:
     :undoc-members:
 
 
+Common Maps
+===========
 
 
+Exponential Map
+---------------
+
+Electrical conductivity varies over many orders of magnitude, so it is a common
+technique when solving the inverse problem to parameterize and optimize in terms
+of log conductivity. This makes sense not only because it ensures all conductivities
+will be positive, but because this is fundamentally the space where conductivity
+lives (i.e. it varies logarithmically).
+
+.. autoclass:: SimPEG.Maps.ExpMap
+    :members:
+    :undoc-members:
+
+
+Vertical 1D Map
+---------------
+
+.. autoclass:: SimPEG.Maps.Vertical1DMap
+    :members:
+    :undoc-members:
+
+
+Map 2D Cross-Section to 3D Model
+--------------------------------
+
+.. autoclass:: SimPEG.Maps.Map2Dto3D
+    :members:
+    :undoc-members:
+
+
+Mesh to Mesh Map
+----------------
+
+
+.. plot::
+
+    from SimPEG import *
+    import matplotlib.pyplot as plt
+    M = Mesh.TensorMesh([100,100])
+    h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
+    h1 = h1/h1.sum()
+    M2 = Mesh.TensorMesh([h1,h1])
+    V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
+    v = Utils.mkvc(V)
+    modh = Maps.Mesh2Mesh([M,M2])
+    modH = Maps.Mesh2Mesh([M2,M])
+    H = modH * v
+    h = modh * H
+    ax = plt.subplot(131)
+    M.plotImage(v, ax=ax)
+    ax.set_title('Fine Mesh (Original)')
+    ax = plt.subplot(132)
+    M2.plotImage(H,clim=[0,1],ax=ax)
+    ax.set_title('Course Mesh')
+    ax = plt.subplot(133)
+    M.plotImage(h,clim=[0,1],ax=ax)
+    ax.set_title('Fine Mesh (Interpolated)')
+    plt.show()
+
+
+.. autoclass:: SimPEG.Maps.Mesh2Mesh
+    :members:
+    :undoc-members:
+
+
+Some Extras
+===========
+
+Combo Map
+---------
+
+The ComboMap holds the information for multiplying and combining
+maps. It also uses the chain rule to create the derivative.
+Remember, any time that you make your own combination of mappings
+be sure to test that the derivative is correct.
+
+.. autoclass:: SimPEG.Maps.ComboMap
+    :members:
+    :undoc-members:
+

docs/examples/DC_Forward_PseudoSection.rst

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

docs/examples/Utils_surface2ind_topo.rst

@@ -1,24 +0,0 @@
-.. _examples_Utils_surface2ind_topo:
-
-.. --------------------------------- ..
-..                                   ..
-..    THIS FILE IS AUTO GENEREATED   ..
-..                                   ..
-..    SimPEG/Examples/__init__.py    ..
-..                                   ..
-.. --------------------------------- ..
-
-
-Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below
-a topographic surface.
-
-
-
-.. plot::
-
-    from SimPEG import Examples
-    Examples.Utils_surface2ind_topo.run()
-
-.. literalinclude:: ../../SimPEG/Examples/Utils_surface2ind_topo.py
-    :language: python
-    :linenos:

tests/base/test_PropMaps.py

@@ -1,7 +1,6 @@
 import unittest
 from SimPEG import *
 from scipy.constants import mu_0
-from SimPEG import Tests
 
 
 class MyPropMap(Maps.PropMap):
@@ -188,34 +187,6 @@ class TestPropMaps(unittest.TestCase):
 
         MyReciprocalPropMap([('sigma', iMap), ('mu', iMap)]) # This should be fine
 
-    def test_linked_derivs_sigma(self):
-        mesh = Mesh.TensorMesh([4,5], x0='CC')
-
-        mapping = Maps.ExpMap(mesh)
-        propmap = MyReciprocalPropMap([('rho', mapping)])
-
-        x0 = np.random.rand(mesh.nC)
-        m  = propmap(x0)
-
-        # test Sigma
-        testme = lambda v: [1./(m.rhoMap*v), m.sigmaDeriv]
-        print 'Testing Rho from Sigma'
-        Tests.checkDerivative(testme, x0, dx=0.01*x0, num=5, plotIt=False)
-
-    def test_linked_derivs_rho(self):
-        mesh = Mesh.TensorMesh([4,5], x0='CC')
-
-        mapping = Maps.ExpMap(mesh)
-        propmap = MyReciprocalPropMap([('sigma', mapping)])
-
-        x0 = np.random.rand(mesh.nC)
-        m  = propmap(x0)
-
-        # test Sigma
-        testme = lambda v: [1./(m.sigmaMap*v), m.rhoDeriv]
-        print 'Testing Rho from Sigma'
-        Tests.checkDerivative(testme, x0, dx=0.01*x0, num=5, plotIt=False)
-
 if __name__ == '__main__':
     unittest.main()
 

tests/base/test_maps.py

@@ -5,10 +5,8 @@ from scipy.sparse.linalg import dsolve
 
 TOL = 1e-14
 
-MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull", "FullMap", "Vertical1DMap", "ParametrizedLayer", "ParametrizedBlockInLayer"]
-MAPS_TO_TEST_3D = [             "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull", "FullMap", "Vertical1DMap", "ParametrizedLayer", "ParametrizedBlockInLayer"]
-MAPS_TO_TEST_CYL = [             "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull", "FullMap", "Vertical1DMap", "ParametrizedLayer"]
-
+MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull","FullMap","Vertical1DMap"]
+MAPS_TO_TEST_3D = [             "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull","FullMap","Vertical1DMap"]
 
 class MapTests(unittest.TestCase):
 
@@ -19,8 +17,6 @@ class MapTests(unittest.TestCase):
         self.mesh2 = Mesh.TensorMesh([a, b], x0=np.array([3, 5]))
         self.mesh3 = Mesh.TensorMesh([a, b, [3,4]], x0=np.array([3, 5, 2]))
         self.mesh22 = Mesh.TensorMesh([b, a], x0=np.array([3, 5]))
-        self.meshCyl = Mesh.CylMesh([10.,1.,10.], x0='00C')
-        print self.meshCyl._meshType
 
     def test_transforms2D(self):
         for M in MAPS_TO_TEST_2D:
@@ -32,15 +28,6 @@ class MapTests(unittest.TestCase):
             maps = getattr(Maps, M)(self.mesh3)
             self.assertTrue(maps.test())
 
-    def test_transformsCyl(self):
-        for M in MAPS_TO_TEST_CYL:
-            maps = getattr(Maps, M)(self.meshCyl)
-            self.assertTrue(maps.test())
-
-    def test_ParametricCasingAndLayer(self):
-        mapping = Maps.ParametrizedCasingAndLayer(self.meshCyl)
-        m = np.r_[-2., 1., 6., 2., -0.1, 0.2, 0.5, 0.2, -0.2, 0.2]
-        self.assertTrue(mapping.test(m))
 
     def test_transforms_logMap_reciprocalMap(self):
         # Note that log/reciprocal maps can be kinda finicky, so we are being explicit about the random seed.