Automate the epsilon picking based on percentile of model values for DEFAULT mode. Fix example.

Fix bug with Maps using array of values

SimPEG/DCIP/DCIPUtils.py

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

SimPEG/Directives.py

@@ -144,12 +144,18 @@ class BetaSchedule(InversionDirective):
             if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
             self.invProb.beta /= self.coolingFactor
 
+
 class TargetMisfit(InversionDirective):
 
+    chifact = 1.
+    phi_d_star = None
+
     @property
     def target(self):
         if getattr(self, '_target', None) is None:
-            self._target = self.survey.nD*0.5
+            if self.phi_d_star is None:
+                self.phi_d_star = 0.5 * self.survey.nD
+            self._target = self.chifact * self.phi_d_star # the factor of 0.5 is because we do phid = 0.5*|| dpred - dobs||^2
         return self._target
     @target.setter
     def target(self, val):
@@ -222,7 +228,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
             mref = 0
         mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
         phi_mx = 0.5 * mx.dot(mx)
-        if self.prob.mesh.dim==2:
+        if self.prob.mesh.dim >= 2:
             my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
             phi_my = 0.5 * my.dot(my)
         else:
@@ -237,47 +243,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         # Save the file as a npz
         np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
 
-#==============================================================================
-# class SaveOutputDictEveryIteration(_SaveEveryIteration):
-#     """SaveOutputDictEveryIteration
-#     A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
-#     """
-#
-#     def initialize(self):
-#         print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName
-#
-#     def endIter(self):
-#         # Save the data.
-#         ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
-#         phi_ms = 0.5*ms.dot(ms)
-#         if self.reg.mrefInSmooth == True:
-#             mref = self.reg.mref
-#         else:
-#             mref = 0
-#         mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
-#         phi_mx = 0.5 * mx.dot(mx)
-#         if self.prob.mesh.dim==2:
-#             my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
-#             phi_my = 0.5 * my.dot(my)
-#         else:
-#             phi_my = 'NaN'
-#         if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType:
-#             mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
-#             phi_mz = 0.5 * mz.dot(mz)
-#         else:
-#             phi_mz = 'NaN'
-#
-#
-#         # Save the file as a npz
-#         np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
-#
-#==============================================================================
-
-# class UpdateReferenceModel(Parameter):
-
-#     mref0 = None
-
-#     def nextIter(self):
 #         mref = getattr(self, 'm_prev', None)
 #         if mref is None:
 #             if self.debug: print 'UpdateReferenceModel is using mref0'
@@ -288,65 +253,177 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 class Update_IRLS(InversionDirective):
 
     eps_min = None
+    eps = None
+    norms = [2.,2.,2.,2.]
     factor = None
     gamma = None
     phi_m_last = None
     phi_d_last = None
+    f_old = None
+    f_min_change = 1e-2
+    beta_tol = 5e-2
+    prctile = 95
+
+    # Solving parameter for IRLS (mode:2)
+    IRLSiter   = 0
+    minGNiter = 5
+    maxIRLSiter = 10
+    iterStart = 0
+
+    # Beta schedule
+    coolingFactor = 2.
+    coolingRate = 1
+
+    mode = 1
+
+    @property
+    def target(self):
+        if getattr(self, '_target', None) is None:
+            self._target = self.survey.nD*0.5
+        return self._target
+    @target.setter
+    def target(self, val):
+        self._target = val
 
     def initialize(self):
 
-        # Scale the regularization for changes in norm
-        if getattr(self, 'phi_m_last', None) is not None:
-
-            self.reg.curModel = self.invProb.curModel
-            self.reg.gamma = 1.
-            phim_new = self.reg.eval(self.invProb.curModel)
-            self.gamma = self.phi_m_last / phim_new
-
-        self.reg.curModel = self.invProb.curModel
-        self.reg.gamma = self.gamma
-
-        if getattr(self, 'phi_d_last', None) is None:
-            self.phi_d_last = self.invProb.phi_d
+        if self.mode == 1:
+            self.reg.norms = [2., 2., 2., 2.]
 
     def endIter(self):
-        # Cool the threshold parameter
-        if getattr(self, 'factor', None) is not None:
-            eps = self.reg.eps / self.factor
 
-            if getattr(self, 'eps_min', None) is not None:
-                self.reg.eps = np.max([self.eps_min,eps])
+        # After reaching target misfit with l2-norm, switch to IRLS (mode:2)
+        if self.invProb.phi_d < self.target and self.mode == 1:
+            print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
+
+            self.mode = 2
+
+            # Either use the supplied epsilon, or fix base on distribution of
+            # model values
+            if getattr(self, 'reg.eps', None) is None:
+                self.reg.eps_p = np.percentile(np.abs(self.invProb.curModel),self.prctile)
+            else:                 
+                self.reg.eps_p = self.eps[0]
+                
+            if getattr(self, 'reg.eps', None) is None:
+                self.reg.eps_q = np.percentile(np.abs(self.reg.regmesh.cellDiffxStencil*(self.reg.mapping * self.invProb.curModel)),self.prctile)
+            else:                 
+                self.reg.eps_q = self.eps[1]
+            
+            print "L[p qx qy qz]-norm : " + str(self.reg.norms)
+            print "eps_p: " + str(self.reg.eps_p) + " eps_q: " + str(self.reg.eps_q)
+            
+            self.reg.norms = self.norms
+            self.coolingFactor = 1.
+            self.coolingRate = 1
+            self.iterStart = self.opt.iter
+            self.phi_d_last = self.invProb.phi_d
+            self.phi_m_last = self.invProb.phi_m_last
+
+            self.reg.l2model = self.invProb.curModel
+            self.reg.curModel = self.invProb.curModel
+
+            if getattr(self, 'f_old', None) is None:
+                self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False)
+
+        # Beta Schedule
+        if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
+            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
+            self.invProb.beta /= self.coolingFactor
+
+
+        # Only update after GN iterations
+        if (self.opt.iter-self.iterStart) % self.minGNiter == 0 and self.mode==2:
+
+            self.IRLSiter += 1
+
+            phim_new = self.reg.eval(self.invProb.curModel)
+            self.f_change = np.abs(self.f_old - phim_new) / self.f_old
+
+            print "Regularization decrease: %6.3e" % (self.f_change)
+
+            # Check for maximum number of IRLS cycles
+            if self.IRLSiter == self.maxIRLSiter:
+                print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
+                self.opt.stopNextIteration = True
+                return
+
+            # Check if the function has changed enough
+            if self.f_change < self.f_min_change and self.IRLSiter > 1:
+                print "Minimum decrease in regularization. End of IRLS"
+                self.opt.stopNextIteration = True
+                return
             else:
-                self.reg.eps = eps
+                self.f_old = phim_new
 
-        # Get phi_m at the end of current iteration
-        self.phi_m_last = self.invProb.phi_m_last
+#            # Cool the threshold parameter if required
+#            if getattr(self, 'factor', None) is not None:
+#                eps = self.reg.eps / self.factor
+#
+#                if getattr(self, 'eps_min', None) is not None:
+#                    self.reg.eps = np.max([self.eps_min,eps])
+#                else:
+#                    self.reg.eps = eps
 
-         # Update the model used for the IRLS weights
-        self.reg.curModel = self.invProb.curModel
+            # Get phi_m at the end of current iteration
+            self.phi_m_last = self.invProb.phi_m_last
 
-        # Temporarely set gamma to 1.
-        self.reg.gamma = 1.
+            # Reset the regularization matrices so that it is
+            # recalculated for current model
+            self.reg._Wsmall = None
+            self.reg._Wx = None
+            self.reg._Wy = None
+            self.reg._Wz = None
 
-        # Compute change in model objective function and update scaling
-        phim_new = self.reg.eval(self.invProb.curModel)
+             # Update the model used for the IRLS weights
+            self.reg.curModel = self.invProb.curModel
 
-        self.reg.gamma = self.phi_m_last / phim_new
+            # Temporarely set gamma to 1. to get raw phi_m
+            self.reg.gamma = 1.
 
-        self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
+            # Compute new model objective function value
+            phim_new = self.reg.eval(self.invProb.curModel)
+
+            # Update gamma to scale the regularization between IRLS iterations
+            self.reg.gamma = self.phi_m_last / phim_new
+
+            # Reset the regularization matrices again for new gamma
+            self.reg._Wsmall = None
+            self.reg._Wx = None
+            self.reg._Wy = None
+            self.reg._Wz = None
+
+            # Check if misfit is within the tolerance, otherwise scale beta
+            val = self.invProb.phi_d / (self.survey.nD*0.5)
+
+            if np.abs(1.-val) > self.beta_tol:
+                self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
 
 class Update_lin_PreCond(InversionDirective):
+    """
+    Create a Jacobi preconditioner for the linear problem
+    """
+    onlyOnStart=False
 
+    def initialize(self):
+
+        if getattr(self.opt, 'approxHinv', None) is None:
+            # Update the pre-conditioner
+            diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
+            self.opt.approxHinv = PC
 
     def endIter(self):
         # Cool the threshold parameter
+        if self.onlyOnStart==True:
+            return
 
         if getattr(self.opt, 'approxHinv', None) is not None:
             # Update the pre-conditioner
-            diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
-            PC     = Utils.sdiag(diagA**-1.)
+            diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
             self.opt.approxHinv = PC
-            print 'Updated pre-cond'
+
 
 class Update_Wj(InversionDirective):
     """

SimPEG/EM/Analytics/DC.py

@@ -2,20 +2,20 @@ import numpy as np
 from scipy.constants import mu_0, pi
 from scipy import special
 
-def DCAnalyticHalf(txloc, rxlocs, sigma, flag="wholespace"):
+def DCAnalyticHalf(txloc, rxlocs, sigma, earth_type="wholespace"):
     """
         Analytic solution for electric potential from a postive pole
 
-        Input variables:
-
-            txloc =  a xyz location of A (+) electrode (np.r_[xa, ya, za])
+        :param array txloc: a xyz location of A (+) electrode (np.r_[xa, ya, za])
+        :param list rxlocs: xyz locations of M (+) and N (-) electrodes [M, N]
 
+            e.g.
             rxlocs = [M, N]
-                M: xyz locations of M (+) electrode (np.c_[xmlocs, ymlocs, zmlocs])
-                N: xyz locations of N (-) electrode (np.c_[xnlocs, ynlocs, znlocs])
+            M: xyz locations of M (+) electrode (np.c_[xmlocs, ymlocs, zmlocs])
+            N: xyz locations of N (-) electrode (np.c_[xnlocs, ynlocs, znlocs])
 
-            sigma = conductivity (either float or complex)
-            flag = "wholsespace" or "halfspace"
+        :param float or complex sigma: values of conductivity
+        :param string earth_type: values of conductivity ("wholsespace" or "halfspace")
 
     """
     M = rxlocs[0]
@@ -28,7 +28,7 @@ def DCAnalyticHalf(txloc, rxlocs, sigma, flag="wholespace"):
     phiN = 1./(4*np.pi*rN*sigma)
     phi = phiM - phiN
 
-    if flag == "halfspace":
+    if earth_type == "halfspace":
         phi *= 2
 
     return phi
@@ -37,27 +37,26 @@ deg2rad  = lambda deg: deg/180.*np.pi
 rad2deg  = lambda rad: rad*180./np.pi
 
 def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
-                 flag = "sec", order=12, halfspace=False):
+                 field_type = "secondary", order=12, halfspace=False):
 # def DCSpherePointCurrent(txloc, rxloc, xc, radius, rho, rho1, \
-#                  flag = "sec", order=12):
+#                  field_type = "secondary", order=12):
     """
 
         Parameters:
 
-            txloc (array) : current electrode location (x,y,z)
-            xc (float)    : x center of depressed sphere
-            rxloc (array) : electrode locations
-                              (Nx3 array, # of electrodes)
-            radius (float): radius of the sphere (m)
-            rho (float)   : resistivity of the background (ohm-m)
-            rho1 (float)  : resistivity of the sphere
-            flag (string) : "sec", "total", "prim"
-                              (default="sec")
-                              "sec": secondary potential only due to sphere
-                              "prim": primary potential from the point source
-                              "total": "sec"+"prim"
-            order (float) : maximum order of Legendre polynomial
-                              (default=12)
+            :param array txloc: A (+) current electrode location (x,y,z)
+            :param array xc: x center of depressed sphere
+            :param array rxloc: M(+) electrode locations / (Nx3 array, # of electrodes)
+
+            :param float radius: radius (float): radius of the sphere (m)
+            :param float rho: resistivity of the background (ohm-m)
+            :param float rho1: resistivity of the sphere
+            :param string field_type: : "secondary", "total", "primary"
+                  (default="secondary")
+                  "secondary": secondary potential only due to sphere
+                  "primary": primary potential from the point source
+                  "total": "secondary"+"primary"
+            :param float order: maximum order of Legendre polynomial (default=12)
 
         Written by Seogi Kang (skang@eos.ubc.ca)
         Ph.D. Candidate of University of British Columbia, Canada
@@ -86,7 +85,7 @@ def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
     # primary potential in a whole space
     prim = rho*1./(4*np.pi*R)
 
-    if flag =="prim":
+    if field_type =="primary":
         return prim
 
     sphind = r < radius
@@ -105,9 +104,9 @@ def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
     else:
         scale = 1
 
-    if flag == "sec":
+    if field_type == "secondary":
         return scale*(out-prim)
-    elif flag == "total":
+    elif field_type == "total":
         return scale*out
 
 def AnBnfun(n, radius, x0, rho, rho1, I=1.):

SimPEG/EM/Base.py

@@ -62,6 +62,15 @@ class BaseEMProblem(Problem.BaseProblem):
             self._Me = self.mesh.getEdgeInnerProduct()
         return self._Me
 
+    @property
+    def MeI(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MeI', None) is None:
+            self._MeI = self.mesh.getEdgeInnerProduct(invMat=True)
+        return self._MeI
+
     @property
     def Mf(self):
         """
@@ -71,13 +80,21 @@ class BaseEMProblem(Problem.BaseProblem):
             self._Mf = self.mesh.getFaceInnerProduct()
         return self._Mf
 
+    @property
+    def MfI(self):
+        """
+            Face inner product matrix
+        """
+        if getattr(self, '_MfI', None) is None:
+            self._MfI = self.mesh.getFaceInnerProduct(invMat=True)
+        return self._MfI
+
     @property
     def Vol(self):
         if getattr(self, '_Vol', None) is None:
             self._Vol = Utils.sdiag(self.mesh.vol)
         return self._Vol
 
-
     # ----- Magnetic Permeability ----- #
     @property
     def MfMui(self):
@@ -152,9 +169,7 @@ class BaseEMProblem(Problem.BaseProblem):
 
         dMeSigmaI_dI = -self.MeSigmaI**2
         dMe_dsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u)
-        dsig_dm = self.curModel.sigmaDeriv
-        return dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm))
-        # return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u)
+        return dMeSigmaI_dI * ( dMe_dsig * self.curModel.sigmaDeriv )
 
     @property
     def MfRho(self):
@@ -170,8 +185,7 @@ class BaseEMProblem(Problem.BaseProblem):
         """
             Derivative of :code:`MfRho` with respect to the model.
         """
-        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv)
-        # self.curModel.rhoDeriv
+        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv
 
     @property
     def MfRhoI(self):
@@ -191,9 +205,7 @@ class BaseEMProblem(Problem.BaseProblem):
 
         dMfRhoI_dI = -self.MfRhoI**2
         dMf_drho = self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u)
-        return dMfRhoI_dI * ( dMf_drho * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) )
-
-        # return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
+        return dMfRhoI_dI * ( dMf_drho * self.curModel.rhoDeriv )
 
 class BaseEMSurvey(Survey.BaseSurvey):
 

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -160,9 +160,9 @@ class Fields(SimPEG.Problem.Fields):
             return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
         return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex)
 
-class Fields_e(Fields):
+class Fields3D_e(Fields):
     """
-    Fields object for Problem_e.
+    Fields object for Problem3D_e.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey
@@ -426,9 +426,9 @@ class Fields_e(Fields):
 
 
 
-class Fields_b(Fields):
+class Fields3D_b(Fields):
     """
-    Fields object for Problem_b.
+    Fields object for Problem3D_b.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey
@@ -693,9 +693,9 @@ class Fields_b(Fields):
         return Zero()
 
 
-class Fields_j(Fields):
+class Fields3D_j(Fields):
     """
-    Fields object for Problem_j.
+    Fields object for Problem3D_j.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey
@@ -988,9 +988,9 @@ class Fields_j(Fields):
         return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) )
 
 
-class Fields_h(Fields):
+class Fields3D_h(Fields):
     """
-    Fields object for Problem_h.
+    Fields object for Problem3D_h.
 
     :param Mesh mesh: mesh
     :param Survey survey: survey

SimPEG/EM/FDEM/ProblemFDEM.py

@@ -1,7 +1,7 @@
 from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
 from SurveyFDEM import Survey as SurveyFDEM
-from FieldsFDEM import Fields, Fields_e, Fields_b, Fields_h, Fields_j
+from FieldsFDEM import Fields, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j
 from SimPEG.EM.Base import BaseEMProblem
 from SimPEG.EM.Utils import omega
 
@@ -17,8 +17,8 @@ class BaseFDEMProblem(BaseEMProblem):
             \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
             {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
 
-        if using the E-B formulation (:code:`Problem_e`
-        or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
+        if using the E-B formulation (:code:`Problem3D_e`
+        or :code:`Problem3D_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity.
 
         If we write Maxwell's equations in terms of
         \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\)
@@ -28,7 +28,7 @@ class BaseFDEMProblem(BaseEMProblem):
             \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\
             \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e}
 
-        if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
+        if using the H-J formulation (:code:`Problem3D_j` or :code:`Problem3D_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity.
 
         The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\)
     """
@@ -87,7 +87,7 @@ class BaseFDEMProblem(BaseEMProblem):
                 du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
 
                 for rx in src.rxList:
-                    df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                    df_dmFun = getattr(f, '_{0}Deriv'.format(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)
             Ainv.clean()
@@ -125,7 +125,7 @@ class BaseFDEMProblem(BaseEMProblem):
                 for rx in src.rxList:
                     PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m
 
-                    df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                    df_duTFun = getattr(f, '_{0}Deriv'.format(rx.projField), None)
                     df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
 
                     ATinvdf_duT = ATinv * df_duT
@@ -137,10 +137,9 @@ class BaseFDEMProblem(BaseEMProblem):
                     df_dmT = df_dmT + du_dmT
 
                     # TODO: this should be taken care of by the reciever?
-                    real_or_imag = rx.projComp
-                    if real_or_imag is 'real':
+                    if rx.component is 'real':
                         Jtv +=   np.array(df_dmT, dtype=complex).real
-                    elif real_or_imag is 'imag':
+                    elif rx.component is 'imag':
                         Jtv += - np.array(df_dmT, dtype=complex).real
                     else:
                         raise Exception('Must be real or imag')
@@ -178,7 +177,7 @@ class BaseFDEMProblem(BaseEMProblem):
 ################################ E-B Formulation #########################################
 ##########################################################################################
 
-class Problem_e(BaseFDEMProblem):
+class Problem3D_e(BaseFDEMProblem):
     """
     By eliminating the magnetic flux density using
 
@@ -200,7 +199,7 @@ class Problem_e(BaseFDEMProblem):
 
     _solutionType = 'eSolution'
     _formulation  = 'EB'
-    fieldsPair    = Fields_e
+    fieldsPair    = Fields3D_e
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -289,7 +288,7 @@ class Problem_e(BaseFDEMProblem):
             return C.T * (MfMui * s_mDeriv(v)) -1j * omega(freq) * s_eDeriv(v)
 
 
-class Problem_b(BaseFDEMProblem):
+class Problem3D_b(BaseFDEMProblem):
     """
     We eliminate :math:`\mathbf{e}` using
 
@@ -311,7 +310,7 @@ class Problem_b(BaseFDEMProblem):
 
     _solutionType = 'bSolution'
     _formulation  = 'EB'
-    fieldsPair    = Fields_b
+    fieldsPair    = Fields3D_b
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -437,7 +436,7 @@ class Problem_b(BaseFDEMProblem):
 ##########################################################################################
 
 
-class Problem_j(BaseFDEMProblem):
+class Problem3D_j(BaseFDEMProblem):
     """
     We eliminate \\\(\\\mathbf{h}\\\) using
 
@@ -459,7 +458,7 @@ class Problem_j(BaseFDEMProblem):
 
     _solutionType = 'jSolution'
     _formulation  = 'HJ'
-    fieldsPair    = Fields_j
+    fieldsPair    = Fields3D_j
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -578,7 +577,7 @@ class Problem_j(BaseFDEMProblem):
 
 
 
-class Problem_h(BaseFDEMProblem):
+class Problem3D_h(BaseFDEMProblem):
     """
     We eliminate \\\(\\\mathbf{j}\\\) using
 
@@ -597,7 +596,7 @@ class Problem_h(BaseFDEMProblem):
 
     _solutionType = 'hSolution'
     _formulation  = 'HJ'
-    fieldsPair    = Fields_h
+    fieldsPair    = Fields3D_h
 
     def __init__(self, mesh, **kwargs):
         BaseFDEMProblem.__init__(self, mesh, **kwargs)

SimPEG/EM/FDEM/RxFDEM.py

@@ -0,0 +1,126 @@
+import SimPEG
+from SimPEG import sp
+
+class BaseRx(SimPEG.Survey.BaseRx):
+    """
+    Frequency domain receiver base class
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        assert(orientation in ['x','y','z']), "Orientation %s not known. Orientation must be in 'x', 'y', 'z'. Arbitrary orientations have not yet been implemented."%orientation
+        assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not %s"%component
+
+        self.projComp = orientation
+        self.component = component
+
+        SimPEG.Survey.BaseRx.__init__(self, locs, rxType=None) #TODO: remove rxType from baseRx
+
+    def projGLoc(self, u):
+        """Grid Location projection (e.g. Ex Fy ...)"""
+        return u._GLoc(self.projField) + self.projComp
+
+    def eval(self, src, mesh, f):
+        """
+        Project fields to recievers to get data.
+
+        :param Source src: FDEM source
+        :param Mesh mesh: mesh used
+        :param Fields f: fields object
+        :rtype: numpy.ndarray
+        :return: fields projected to recievers
+        """
+
+        P = self.getP(mesh, self.projGLoc(f))
+        f_part_complex = f[src, self.projField]
+        f_part = getattr(f_part_complex, self.component) # get the real or imag component
+
+        return P*f_part
+
+    def evalDeriv(self, src, mesh, f, v, adjoint=False):
+        """
+        Derivative of projected fields with respect to the inversion model times a vector.
+
+        :param Source src: FDEM source
+        :param Mesh mesh: mesh used
+        :param Fields f: fields object
+        :param numpy.ndarray v: vector to multiply
+        :rtype: numpy.ndarray
+        :return: fields projected to recievers
+        """
+
+        P = self.getP(mesh, self.projGLoc(f))
+
+        if not adjoint:
+            Pv_complex = P * v
+            Pv = getattr(Pv_complex, self.component)
+        elif adjoint:
+            Pv_real = P.T * v
+
+            if self.component == 'imag':
+                Pv = 1j*Pv_real
+            elif self.component == 'real':
+                Pv = Pv_real.astype(complex)
+            else:
+                raise NotImplementedError('must be real or imag')
+
+        return Pv
+
+
+class Point_e(BaseRx):
+    """
+    Electric field FDEM receiver
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        self.projField = 'e'
+        super(Point_e, self).__init__(locs, orientation, component)
+
+
+class Point_b(BaseRx):
+    """
+    Magnetic flux FDEM receiver
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        self.projField = 'b'
+        super(Point_b, self).__init__(locs, orientation, component)
+
+
+class Point_h(BaseRx):
+    """
+    Magnetic field FDEM receiver
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        self.projField = 'h'
+        super(Point_h, self).__init__(locs, orientation, component)
+
+
+class Point_j(BaseRx):
+    """
+    Current density FDEM receiver
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        self.projField = 'j'
+        super(Point_j, self).__init__(locs, orientation, component)

SimPEG/EM/FDEM/SrcFDEM.py

@@ -9,8 +9,14 @@ class BaseSrc(Survey.BaseSrc):
     """
 
     freq = None
-    # rxPair = RxFDEM
-    integrate = True
+    integrate = False
+    _ePrimary = None
+    _bPrimary = None
+    _hPrimary = None
+    _jPrimary = None
+
+    def __init__(self, rxList, **kwargs):
+        Survey.BaseSrc.__init__(self, rxList, **kwargs)
 
     def eval(self, prob):
         """
@@ -50,7 +56,9 @@ class BaseSrc(Survey.BaseSrc):
         :rtype: numpy.ndarray
         :return: primary magnetic flux density
         """
-        return Zero()
+        if self._bPrimary is None:
+            return Zero()
+        return self._bPrimary
 
     def hPrimary(self, prob):
         """
@@ -60,7 +68,9 @@ class BaseSrc(Survey.BaseSrc):
         :rtype: numpy.ndarray
         :return: primary magnetic field
         """
-        return Zero()
+        if self._hPrimary is None:
+            return Zero()
+        return self._hPrimary
 
     def ePrimary(self, prob):
         """
@@ -70,7 +80,9 @@ class BaseSrc(Survey.BaseSrc):
         :rtype: numpy.ndarray
         :return: primary electric field
         """
-        return Zero()
+        if self._ePrimary is None:
+            return Zero()
+        return self._ePrimary
 
     def jPrimary(self, prob):
         """
@@ -80,7 +92,9 @@ class BaseSrc(Survey.BaseSrc):
         :rtype: numpy.ndarray
         :return: primary current density
         """
-        return Zero()
+        if self._jPrimary is None:
+            return Zero()
+        return self._jPrimary
 
     def s_m(self, prob):
         """
@@ -135,15 +149,14 @@ class RawVec_e(BaseSrc):
     :param list rxList: receiver list
     :param float freq: frequency
     :param numpy.array s_e: electric source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
 
-    def __init__(self, rxList, freq, s_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
+    def __init__(self, rxList, freq, s_e, **kwargs):
         self._s_e = np.array(s_e, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
 
-        BaseSrc.__init__(self, rxList)
+        BaseSrc.__init__(self, rxList, **kwargs)
 
     def s_e(self, prob):
         """
@@ -165,15 +178,14 @@ class RawVec_m(BaseSrc):
     :param float freq: frequency
     :param rxList: receiver list
     :param numpy.array s_m: magnetic source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
 
-    def __init__(self, rxList, freq, s_m, integrate=True):  #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
+    def __init__(self, rxList, freq, s_m, **kwargs):  #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
         self._s_m = np.array(s_m, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
 
-        BaseSrc.__init__(self, rxList)
+        BaseSrc.__init__(self, rxList, **kwargs)
 
     def s_m(self, prob):
         """
@@ -196,14 +208,13 @@ class RawVec(BaseSrc):
     :param float freq: frequency
     :param numpy.array s_m: magnetic source term
     :param numpy.array s_e: electric source term
-    :param bool integrate: Integrate the source term (multiply by Me) [True]
+    :param bool integrate: Integrate the source term (multiply by Me) [False]
     """
-    def __init__(self, rxList, freq, s_m, s_e, integrate=True):
+    def __init__(self, rxList, freq, s_m, s_e, **kwargs):
         self._s_m = np.array(s_m, dtype=complex)
         self._s_e = np.array(s_e, dtype=complex)
         self.freq = float(freq)
-        self.integrate = integrate
-        BaseSrc.__init__(self, rxList)
+        BaseSrc.__init__(self, rxList, **kwargs)
 
     def s_m(self, prob):
         """
@@ -277,14 +288,13 @@ class MagDipole(BaseSrc):
     :param float mu: background magnetic permeability
     """
 
-    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0):
+    def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0, **kwargs):
         self.freq = float(freq)
         self.loc = loc
         self.orientation = orientation
         assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
         self.moment = moment
         self.mu = mu
-        self.integrate = False
         BaseSrc.__init__(self, rxList)
 
     def bPrimary(self, prob):
@@ -542,7 +552,7 @@ class CircularLoop(BaseSrc):
             if not prob.mesh.isSymmetric:
                 # TODO ?
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
+            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
             srcfct = MagneticDipoleVectorPotential

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -4,126 +4,9 @@ from SimPEG.EM.Base import BaseEMSurvey
 from scipy.constants import mu_0
 from SimPEG.Utils import Zero, Identity
 import SrcFDEM as Src
+import RxFDEM as Rx
 from SimPEG import sp
 
-
-####################################################
-# Receivers
-####################################################
-
-class Rx(SimPEG.Survey.BaseRx):
-    """
-    Frequency domain receivers
-
-    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
-    :param string rxType: reciever type from knownRxTypes
-    """
-
-    knownRxTypes = {
-                    'exr':['e', 'x', 'real'],
-                    'eyr':['e', 'y', 'real'],
-                    'ezr':['e', 'z', 'real'],
-                    'exi':['e', 'x', 'imag'],
-                    'eyi':['e', 'y', 'imag'],
-                    'ezi':['e', 'z', 'imag'],
-
-                    'bxr':['b', 'x', 'real'],
-                    'byr':['b', 'y', 'real'],
-                    'bzr':['b', 'z', 'real'],
-                    'bxi':['b', 'x', 'imag'],
-                    'byi':['b', 'y', 'imag'],
-                    'bzi':['b', 'z', 'imag'],
-
-                    'jxr':['j', 'x', 'real'],
-                    'jyr':['j', 'y', 'real'],
-                    'jzr':['j', 'z', 'real'],
-                    'jxi':['j', 'x', 'imag'],
-                    'jyi':['j', 'y', 'imag'],
-                    'jzi':['j', 'z', 'imag'],
-
-                    'hxr':['h', 'x', 'real'],
-                    'hyr':['h', 'y', 'real'],
-                    'hzr':['h', 'z', 'real'],
-                    'hxi':['h', 'x', 'imag'],
-                    'hyi':['h', 'y', 'imag'],
-                    'hzi':['h', 'z', 'imag'],
-                   }
-    radius = None
-
-    def __init__(self, locs, rxType):
-        SimPEG.Survey.BaseRx.__init__(self, locs, rxType)
-
-    @property
-    def projField(self):
-        """Field Type projection (e.g. e b ...)"""
-        return self.knownRxTypes[self.rxType][0]
-
-    @property
-    def projComp(self):
-        """Component projection (real/imag)"""
-        return self.knownRxTypes[self.rxType][2]
-
-    def projGLoc(self, f):
-        """Grid Location projection (e.g. Ex Fy ...)"""
-        return f._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
-
-    def eval(self, src, mesh, f):
-        """
-        Project fields to recievers to get data.
-
-        :param Source src: FDEM source
-        :param Mesh mesh: mesh used
-        :param Fields f: fields object
-        :rtype: numpy.ndarray
-        :return: fields projected to recievers
-        """
-        # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
-        # projGLoc += self.knownRxTypes[self.rxType][1]
-
-        P = self.getP(mesh, self.projGLoc(f))
-        f_part_complex = f[src, self.projField]
-        # get the real or imag component
-        real_or_imag = self.projComp
-        f_part = getattr(f_part_complex, real_or_imag)
-
-        return P*f_part
-
-    def evalDeriv(self, src, mesh, f, v, adjoint=False):
-        """
-        Derivative of projected fields with respect to the inversion model times a vector.
-
-        :param Source src: FDEM source
-        :param Mesh mesh: mesh used
-        :param Fields f: fields object
-        :param numpy.ndarray v: vector to multiply
-        :rtype: numpy.ndarray
-        :return: fields projected to recievers
-        """
-
-        P = self.getP(mesh, self.projGLoc(f))
-
-        if not adjoint:
-            Pv_complex = P * v
-            real_or_imag = self.projComp
-            Pv = getattr(Pv_complex, real_or_imag)
-        elif adjoint:
-            Pv_real = P.T * v
-
-            real_or_imag = self.projComp
-            if real_or_imag == 'imag':
-                Pv = 1j*Pv_real
-            elif real_or_imag == 'real':
-                Pv = Pv_real.astype(complex)
-            else:
-                raise NotImplementedError('must be real or imag')
-
-        return Pv
-
-
-####################################################
-# Survey
-####################################################
-
 class Survey(BaseEMSurvey):
     """
     Frequency domain electromagnetic survey
@@ -132,7 +15,7 @@ class Survey(BaseEMSurvey):
     """
 
     srcPair = Src.BaseSrc
-    rxPair = Rx
+    rxPair = Rx.BaseRx
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency

SimPEG/EM/FDEM/__init__.py

@@ -1,3 +1,5 @@
-from SurveyFDEM import Rx, Src, Survey
-from FDEM import BaseFDEMProblem, Problem_e, Problem_b, Problem_j, Problem_h
-from FieldsFDEM import *
+from SurveyFDEM import Survey
+import SrcFDEM as Src
+import RxFDEM as Rx
+from ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h
+from FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h

SimPEG/EM/Static/DC/ProblemDC.py

@@ -122,13 +122,12 @@ class Problem3D_CC(BaseDCProblem):
 
         Make the A matrix for the cell centered DC resistivity problem
 
-        A = D MfRhoI D^\\top V
+        A = D MfRhoI G
 
         """
 
         D = self.Div
         G = self.Grad
-        # TODO: this won't work for full anisotropy
         MfRhoI = self.MfRhoI
         A = D * MfRhoI * G
 
@@ -144,13 +143,8 @@ class Problem3D_CC(BaseDCProblem):
         MfRhoIDeriv = self.MfRhoIDeriv
 
         if adjoint:
-            # if self._makeASymmetric is True:
-            #     v = V * v
             return(MfRhoIDeriv( G * u ).T) * ( D.T * v)
 
-        # I think we should deprecate this for DC problem.
-        # if self._makeASymmetric is True:
-        #     return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) )
         return D * (MfRhoIDeriv( G * u ) * v)
 
     def getRHS(self):
@@ -162,10 +156,6 @@ class Problem3D_CC(BaseDCProblem):
 
         RHS = self.getSourceTerm()
 
-        # I think we should deprecate this for DC problem.
-        # if self._makeASymmetric is True:
-        #     return self.Vol.T * RHS
-
         return RHS
 
     def getRHSDeriv(self, src, v, adjoint=False):
@@ -255,11 +245,10 @@ class Problem3D_N(BaseDCProblem):
 
         Make the A matrix for the cell centered DC resistivity problem
 
-        A = D MfRhoI D^\\top V
+        A = G.T MeSigma G
 
         """
 
-        # TODO: this won't work for full anisotropy
         MeSigma = self.MeSigma
         Grad = self.mesh.nodalGrad
         A = Grad.T * MeSigma * Grad

SimPEG/EM/Static/DC/ProblemDC_2D.py

@@ -161,14 +161,13 @@ class Problem2D_CC(BaseDCProblem_2D):
 
         Make the A matrix for the cell centered DC resistivity problem
 
-        A = D MfRhoI D^\\top V
+        A = D MfRhoI G
 
         """
 
         D = self.Div
         G = self.Grad
         vol = self.mesh.vol
-        # TODO: this won't work for full anisotropy
         MfRhoI = self.MfRhoI
         # Get resistivity rho
         rho = self.curModel.rho
@@ -304,11 +303,10 @@ class Problem2D_N(BaseDCProblem_2D):
 
         Make the A matrix for the cell centered DC resistivity problem
 
-        A = D MfRhoI D^\\top V
+        A = D MfRhoI G
 
         """
 
-        # TODO: this won't work for full anisotropy
         MeSigma = self.MeSigma
         MnSigma = self.MnSigma
         Grad = self.mesh.nodalGrad

SimPEG/EM/Static/IP/ProblemIP.py

@@ -180,13 +180,12 @@ class Problem3D_CC(BaseIPProblem):
 
         Make the A matrix for the cell centered DC resistivity problem
 
-        A = D MfRhoI D^\\top V
+        A = D MfRhoI G
 
         """
 
         D = self.Div
         G = self.Grad
-        # TODO: this won't work for full anisotropy
         MfRhoI = self.MfRhoI
         A = D * MfRhoI * G
 
@@ -313,11 +312,10 @@ class Problem3D_N(BaseIPProblem):
 
         Make the A matrix for the cell centered DC resistivity problem
 
-        A = D MfRhoI D^\\top V
+        A = G.T MeSigma G
 
         """
 
-        # TODO: this won't work for full anisotropy
         MeSigma = self.MeSigma
         Grad = self.mesh.nodalGrad
         A = Grad.T * MeSigma * Grad

SimPEG/EM/Static/SIP/ProblemSIP.py

@@ -251,7 +251,7 @@ class Problem3D_CC(BaseSIPProblem):
 
         Make the A matrix for the cell centered DC resistivity problem
 
-        A = D MfRhoI D^\\top V
+        A = D MfRhoI G
 
         """
 
@@ -384,7 +384,7 @@ class Problem3D_N(BaseSIPProblem):
 
         Make the A matrix for the cell centered DC resistivity problem
 
-        A = D MfRhoI D^\\top V
+        A = G.T MeSigma G
 
         """
 

SimPEG/EM/Utils/testingUtils.py

@@ -20,56 +20,61 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
     mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
 
     if useMu is True:
-        mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] 
+        mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))]
     else:
         mapping = Maps.ExpMap(mesh)
 
     x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid
     XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5))
-    Rx0 = EM.FDEM.Rx(XYZ, comp)
+    Rx0 = getattr(EM.FDEM.Rx, 'Point_' + comp[0])
+    if comp[2] == 'r':
+        real_or_imag = 'real'
+    elif comp[2] == 'i':
+        real_or_imag = 'imag'
+    rx0 = Rx0(XYZ, comp[1], 'imag')
 
     Src = []
 
     for SrcType in SrcList:
         if SrcType is 'MagDipole':
-            Src.append(EM.FDEM.Src.MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
+            Src.append(EM.FDEM.Src.MagDipole([rx0], freq=freq, loc=np.r_[0.,0.,0.]))
         elif SrcType is 'MagDipole_Bfield':
-            Src.append(EM.FDEM.Src.MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
+            Src.append(EM.FDEM.Src.MagDipole_Bfield([rx0], freq=freq, loc=np.r_[0.,0.,0.]))
         elif SrcType is 'CircularLoop':
-            Src.append(EM.FDEM.Src.CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]))
+            Src.append(EM.FDEM.Src.CircularLoop([rx0], freq=freq, loc=np.r_[0.,0.,0.]))
         elif SrcType is 'RawVec':
             if fdemType is 'e' or fdemType is 'b':
                 S_m = np.zeros(mesh.nF)
                 S_e = np.zeros(mesh.nE)
                 S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
                 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
-                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
+                Src.append(EM.FDEM.Src.RawVec([rx0], freq, S_m, mesh.getEdgeInnerProduct()*S_e))
 
             elif fdemType is 'h' or fdemType is 'j':
                 S_m = np.zeros(mesh.nE)
                 S_e = np.zeros(mesh.nF)
                 S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
                 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
-                Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
+                Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*S_m, S_e))
 
     if verbose:
         print '  Fetching %s problem' % (fdemType)
 
     if fdemType == 'e':
         survey = EM.FDEM.Survey(Src)
-        prb = EM.FDEM.Problem_e(mesh, mapping=mapping)
+        prb = EM.FDEM.Problem3D_e(mesh, mapping=mapping)
 
     elif fdemType == 'b':
         survey = EM.FDEM.Survey(Src)
-        prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
+        prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping)
 
     elif fdemType == 'j':
         survey = EM.FDEM.Survey(Src)
-        prb = EM.FDEM.Problem_j(mesh, mapping=mapping)
+        prb = EM.FDEM.Problem3D_j(mesh, mapping=mapping)
 
     elif fdemType == 'h':
         survey = EM.FDEM.Survey(Src)
-        prb = EM.FDEM.Problem_h(mesh, mapping=mapping)
+        prb = EM.FDEM.Problem3D_h(mesh, mapping=mapping)
 
     else:
         raise NotImplementedError()
@@ -90,7 +95,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM
     prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
     mesh = prb1.mesh
     print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
-    
+
     logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
     mu = np.ones(mesh.nC)*MU
 

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

SimPEG/Examples/EM_FDEM_1D_Inversion.py

@@ -42,8 +42,8 @@ def run(plotIt=True):
         ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
 
 
-    rxOffset=10. 
-    bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi')
+    rxOffset=10.
+    bzi = EM.FDEM.Rx.Point_b(np.array([[rxOffset, 0., 1e-3]]), orientation='z', component='imag')
 
     freqs = np.logspace(1,3,10)
     srcLoc = np.array([0., 0., 10.])
@@ -51,7 +51,7 @@ def run(plotIt=True):
     srcList = [EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z') for freq in freqs]
 
     survey = EM.FDEM.Survey(srcList)
-    prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
+    prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping)
 
     try:
         from pymatsolver import MumpsSolver

SimPEG/Examples/EM_Schenkel_Morrison_Casing.py

@@ -215,7 +215,7 @@ def run(plotIt=True):
     # ------------ Problem and Survey ---------------
     survey = FDEM.Survey(sg_p + dg_p)
     mapping = [('sigma', Maps.IdentityMap(mesh))]
-    problem = FDEM.Problem_h(mesh, mapping=mapping)
+    problem = FDEM.Problem3D_h(mesh, mapping=mapping)
     problem.pair(survey)
 
     # ------------- Solve ---------------------------

SimPEG/Examples/Inversion_IRLS.py

@@ -1,7 +1,7 @@
 from SimPEG import *
 
 
-def run(N=200, plotIt=True):
+def run(N=100, plotIt=True):
     """
         Inversion: Linear Problem
         =========================
@@ -18,6 +18,8 @@ def run(N=200, plotIt=True):
     mesh = Mesh.TensorMesh([N])
 
     m0 = np.ones(mesh.nC) * 1e-4
+    mref = np.zeros(mesh.nC)
+
     nk = 10
     jk = np.linspace(1.,nk,nk)
     p = -2.
@@ -40,67 +42,35 @@ def run(N=200, plotIt=True):
     survey = Survey.LinearSurvey()
     survey.pair(prob)
     survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk)
-    #survey.makeSyntheticData(mtrue, std=std_noise)
 
     wd = np.ones(nk) * std_noise
 
-    #print survey.std[0]
-    #M = prob.mesh
     # Distance weighting
     wr = np.sum(prob.G**2.,axis=0)**0.5
     wr = ( wr/np.max(wr) )
 
-    reg = Regularization.Simple(mesh)
-    reg.wght = wr
-
     dmis = DataMisfit.l2_DataMisfit(survey)
     dmis.Wd = 1./wd
 
-    opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
-    invProb.curModel = m0
-
-    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
-    target = Directives.TargetMisfit()
-
     betaest = Directives.BetaEstimate_ByEig()
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
-
-
-    mrec = inv.run(m0)
-    ml2 = mrec
-    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
-
-    # Switch regularization to sparse
-    phim = invProb.phi_m_last
-    phid =  invProb.phi_d
 
     reg = Regularization.Sparse(mesh)
+    reg.mref = mref
+    reg.cell_weights = wr
 
-#==============================================================================
-#     fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-#     dmdx = reg.mesh.cellDiffxStencil * mrec
-#     plt.plot(np.sort(dmdx))
-#==============================================================================
-
-    #reg.recModel = mrec
-    reg.wght = np.ones(mesh.nC)
     reg.mref = np.zeros(mesh.nC)
-    reg.eps_p = 2e-3
-    reg.eps_q = 2e-3
-    reg.norms   = [0., 0., 2., 2.]
-    reg.wght = wr
+    
 
-    opt = Optimization.ProjectedGNCG(maxIter=5 ,lower=-2.,upper=2., maxIterCG= 100, tolCG = 1e-3)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.)
-    beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
-    #betaest = Directives.BetaEstimate_ByEig()
-    target = Directives.TargetMisfit()
-    IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
+    opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
+    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
+    update_Jacobi = Directives.Update_lin_PreCond()
+    
+    # Set the IRLS directive, penalize the lowest 25 percentile of model values
+    # Start with an l2-l2, then switch to lp-norms
+    norms   = [0., 0., 2., 2.]    
+    IRLS = Directives.Update_IRLS( norms=norms, prctile = 25, maxIRLSiter = 15, minGNiter=3)
 
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
-
-    m0 = mrec
+    inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
 
     # Run inversion
     mrec = inv.run(m0)
@@ -117,7 +87,7 @@ def run(N=200, plotIt=True):
         axes[0].set_title('Columns of matrix G')
 
         axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
-        axes[1].plot(mesh.vectorCCx, ml2, 'r-')
+        axes[1].plot(mesh.vectorCCx, reg.l2model, 'r-')
         #axes[1].legend(('True Model', 'Recovered Model'))
         axes[1].set_ylim(-1.0,1.25)
 

SimPEG/Examples/Mesh_Basic_ForwardDC.py

@@ -1,22 +1,25 @@
 from SimPEG import Mesh, Utils, np, SolverLU
 
-## 2D DC forward modeling example with Tensor and Curvilinear Meshes
-
 def run(plotIt=True):
+
+    """
+        Mesh: Basic Forward 2D DC Resistivity
+        =====================================
+
+        2D DC forward modeling example with Tensor and Curvilinear Meshes
+    """
+
     # Step1: Generate Tensor and Curvilinear Mesh
     sz = [40,40]
-    # Tensor Mesh
     tM = Mesh.TensorMesh(sz)
-    # Curvilinear Mesh
     rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate'))
+
     # Step2: Direct Current (DC) operator
     def DCfun(mesh, pts):
         D = mesh.faceDiv
-        G = D.T
         sigma = 1e-2*np.ones(mesh.nC)
-        Msigi = mesh.getFaceInnerProduct(1./sigma)
-        MsigI = Utils.sdInv(Msigi)
-        A = D*MsigI*G
+        MsigI = mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True)
+        A = -D*MsigI*D.T
         A[-1,-1] /= mesh.vol[-1] # Remove null space
         rhs = np.zeros(mesh.nC)
         txind = Utils.meshutils.closestPoints(mesh, pts)
@@ -37,39 +40,17 @@ def run(plotIt=True):
     if not plotIt: return
 
     import matplotlib.pyplot as plt
-    import matplotlib
-    from matplotlib.mlab import griddata
 
     #Step4: Making Figure
     fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-    label = ["(a)", "(b)"]
-    opts = {}
     vmin, vmax = phitM.min(), phitM.max()
     dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True)
-
-    #TODO: At the moment Curvilinear Mesh do not have plotimage
-
-    Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F')
-    Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F')
-    PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear')
-    axes[1].contourf(Xi, Yi, PHIrM, 100, vmin=vmin, vmax=vmax)
-
+    dat = rM.plotImage(phirM, ax=axes[1], clim=(vmin, vmax), grid=True)
     cb = plt.colorbar(dat[0], ax=axes[0]); cb.set_label("Voltage (V)")
     cb = plt.colorbar(dat[0], ax=axes[1]); cb.set_label("Voltage (V)")
 
-    tM.plotGrid(ax=axes[0], **opts)
     axes[0].set_title('TensorMesh')
-    rM.plotGrid(ax=axes[1], **opts)
     axes[1].set_title('CurvilinearMesh')
-    for i in range(2):
-        axes[i].set_xlim(0.025, 0.975)
-        axes[i].set_ylim(0.025, 0.975)
-        axes[i].text(0., 1.0, label[i], fontsize=20)
-        if i==0:
-            axes[i].set_ylabel("y")
-        else:
-            axes[i].set_ylabel(" ")
-        axes[i].set_xlabel("x")
     plt.show()
 
 

SimPEG/Examples/Utils_surface2ind_topo.py

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

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

SimPEG/MT/BaseMT.py

@@ -1,5 +1,5 @@
 from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
-from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem
+from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
 from SurveyMT import Survey, Data
 from FieldsMT import BaseMTFields
 

SimPEG/MT/Utils/ediFilesUtils.py

@@ -7,17 +7,16 @@ from SimPEG.MT.Utils.dataUtils import rec2ndarr
 # Import modules
 import numpy as np
 import os, sys, re
-try:
-    import osr
-except ImportError as e:
-    print 'Could not import osr, missing the gdal package'
-    pass
+
 
 class EDIimporter:
     """
     A class to import EDIfiles.
 
     """
+
+
+    # Define data converters
     _impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit
     _impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)
 
@@ -26,8 +25,8 @@ class EDIimporter:
     comps = None
 
     # Hidden properties
-    _outEPSG = None
-    _2out = None
+    _outEPSG = None # Project info
+    _2out = None    # The projection operator
 
 
     def __init__(self, EDIfilesList, compList=None, outEPSG=None):
@@ -113,6 +112,12 @@ class EDIimporter:
     # nOutData=length(obj.data);
     # obj.data(nOutData+1:nOutData+length(TEMP.data),:) = TEMP.data;
     def _transfromPoints(self,longD,latD):
+        # Import the coordinate projections
+        try:
+            import osr
+        except ImportError as e:
+            print 'Could not import osr, missing the gdal package\nCan not project coordinates'
+            raise e
         # Coordinates convertor
         if self._2out is None:
             src = osr.SpatialReference()

SimPEG/Maps.py

@@ -502,7 +502,9 @@ class InjectActiveCells(IdentityMap):
         if Utils.isScalar(valInactive):
             self.valInactive = np.ones(self.nC)*float(valInactive)
         else:
-            self.valInactive = valInactive.copy()
+            self.valInactive = np.ones(self.nC)
+            self.valInactive[self.indInactive] = valInactive.copy()
+            
         self.valInactive[self.indActive] = 0
 
         inds = np.nonzero(self.indActive)[0]
@@ -533,83 +535,6 @@ 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):
     """

SimPEG/Mesh/CurvilinearMesh.py

@@ -2,6 +2,7 @@ from SimPEG import Utils, np
 from BaseMesh import BaseRectangularMesh
 from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
+from View import CurvView
 
 # Some helper functions.
 length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
@@ -10,7 +11,7 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
 normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
 
 
-class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
+class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView):
     """
     CurvilinearMesh is a mesh class that deals with curvilinear meshes.
 
@@ -330,102 +331,6 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts):
 
 
 
-    #############################################
-    #            Plotting Functions             #
-    #############################################
-
-    def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True,  showIt=False):
-        """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
-
-
-        .. plot::
-            :include-source:
-
-            from SimPEG import Mesh, Utils
-            X, Y = Utils.exampleLrmGrid([3,3],'rotate')
-            M = Mesh.CurvilinearMesh([X, Y])
-            M.plotGrid(showIt=True)
-
-        """
-        import matplotlib.pyplot as plt
-        import matplotlib
-        from mpl_toolkits.mplot3d import Axes3D
-        mkvc = Utils.mkvc
-
-        axOpts = {'projection':'3d'} if self.dim == 3 else {}
-        if ax is None: ax = plt.subplot(111, **axOpts)
-
-        NN = self.r(self.gridN, 'N', 'N', 'M')
-        if self.dim == 2:
-
-            if lines:
-                X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
-                Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
-
-                X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
-                Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
-
-                X = np.r_[X1, X2]
-                Y = np.r_[Y1, Y2]
-
-                ax.plot(X, Y, 'b-')
-            if centers:
-                ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
-
-            # Nx = self.r(self.normals, 'F', 'Fx', 'V')
-            # Ny = self.r(self.normals, 'F', 'Fy', 'V')
-            # Tx = self.r(self.tangents, 'E', 'Ex', 'V')
-            # Ty = self.r(self.tangents, 'E', 'Ey', 'V')
-
-            # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
-
-            # nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
-            # ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
-            # ax.plot(nX, nY, 'r-')
-
-            # nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
-            # #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
-            # ax.plot(nX, nY, 'g-')
-
-            # tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
-            # tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
-            # ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
-            # ax.plot(tX, tY, 'r-')
-
-            # nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
-            # nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
-            # #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
-            # ax.plot(nX, nY, 'g-')
-
-        elif self.dim == 3:
-            X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
-            Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
-            Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
-
-            X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten()
-            Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten()
-            Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten()
-
-            X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten()
-            Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten()
-            Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten()
-
-            X = np.r_[X1, X2, X3]
-            Y = np.r_[Y1, Y2, Y3]
-            Z = np.r_[Z1, Z2, Z3]
-
-            ax.plot(X, Y, 'b', zs=Z)
-            ax.set_zlabel('x3')
-
-        ax.grid(True)
-        ax.set_xlabel('x1')
-        ax.set_ylabel('x2')
-
-        if showIt: plt.show()
-
-
 if __name__ == '__main__':
     nc = 5
     h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])

SimPEG/Mesh/CylMesh.py

@@ -330,7 +330,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
                 raise NotImplementedError('wrapping in the averaging is not yet implemented')
         return self._aveF2CCV
 
-    def getInterpolationMatCartMesh(self, Mrect, locType='CC'):
+    def getInterpolationMatCartMesh(self, Mrect, locType='CC', locTypeTo=None):
         """
             Takes a cartesian mesh and returns a projection to translate onto the cartesian grid.
         """
@@ -338,19 +338,22 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
         assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes"
 
 
+        if locTypeTo is None:
+            locTypeTo = locType
+
         if locType == 'F':
             # do this three times for each component
-            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx')
-            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy')
-            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz')
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Fx', locTypeTo=locTypeTo+'x')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy', locTypeTo=locTypeTo+'y')
+            Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz', locTypeTo=locTypeTo+'z')
             return sp.vstack((X,Y,Z))
         if locType == 'E':
-            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex')
-            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey')
-            Z = spzeros(Mrect.nEz, self.nE)
+            X = self.getInterpolationMatCartMesh(Mrect, locType='Ex', locTypeTo=locTypeTo+'x')
+            Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey', locTypeTo=locTypeTo+'y')
+            Z = spzeros(getattr(Mrect, 'n' + locTypeTo + 'z'), self.nE)
             return sp.vstack((X,Y,Z))
 
-        grid = getattr(Mrect, 'grid' + locType)
+        grid = getattr(Mrect, 'grid' + locTypeTo)
         # This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice
         theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2
         theta[theta < 0] += np.pi*2.0
@@ -366,7 +369,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView):
                         'Ex': Mrect.tangents[:Mrect.nEx,:],
                         'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:],
                         'Ez': Mrect.tangents[-Mrect.nEz:,:],
-                    }[locType]
+                    }[locTypeTo]
             if 'F' in locType:
                 normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)]
                 proj = ( normals * dotMe ).sum(axis=1)

SimPEG/Mesh/MeshIO.py

@@ -24,7 +24,6 @@ class TensorMeshIO(object):
                     re = int(sp[0])*(' ' + sp[1])
                     line = line.replace(st,re.strip())
             return np.array(line.split(),dtype=float)
-
         # Read the file as line strings, remove lines with comment = !
         msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!')
 

SimPEG/Mesh/View.py

@@ -552,7 +552,8 @@ class CurvView(object):
     def __init__(self):
         pass
 
-    def plotGrid(self, length=0.05, showIt=False):
+
+    def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True,  showIt=False):
         """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
 
 
@@ -560,60 +561,63 @@ class CurvView(object):
             :include-source:
 
             from SimPEG import Mesh, Utils
-            X, Y = Utils.exampleCurvGird([3,3],'rotate')
+            X, Y = Utils.exampleLrmGrid([3,3],'rotate')
             M = Mesh.CurvilinearMesh([X, Y])
             M.plotGrid(showIt=True)
 
         """
+        import matplotlib.pyplot as plt
+        import matplotlib
+        from mpl_toolkits.mplot3d import Axes3D
+
+        axOpts = {'projection':'3d'} if self.dim == 3 else {}
+        if ax is None: ax = plt.subplot(111, **axOpts)
+
         NN = self.r(self.gridN, 'N', 'N', 'M')
         if self.dim == 2:
-            fig = plt.figure(2)
-            fig.clf()
-            ax = plt.subplot(111)
-            X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
-            Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
 
-            X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
-            Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
+            if lines:
+                X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten()
+                Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten()
 
-            X = np.r_[X1, X2]
-            Y = np.r_[Y1, Y2]
+                X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten()
+                Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten()
 
-            plt.plot(X, Y)
+                X = np.r_[X1, X2]
+                Y = np.r_[Y1, Y2]
 
-            plt.hold(True)
-            Nx = self.r(self.normals, 'F', 'Fx', 'V')
-            Ny = self.r(self.normals, 'F', 'Fy', 'V')
-            Tx = self.r(self.tangents, 'E', 'Ex', 'V')
-            Ty = self.r(self.tangents, 'E', 'Ey', 'V')
+                ax.plot(X, Y, 'b-')
+            if centers:
+                ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro')
 
-            plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
+            # Nx = self.r(self.normals, 'F', 'Fx', 'V')
+            # Ny = self.r(self.normals, 'F', 'Fy', 'V')
+            # Tx = self.r(self.tangents, 'E', 'Ex', 'V')
+            # Ty = self.r(self.tangents, 'E', 'Ey', 'V')
 
-            nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
-            nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
-            plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
-            plt.plot(nX, nY, 'r-')
+            # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo')
 
-            nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
-            nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
-            #plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
-            plt.plot(nX, nY, 'g-')
+            # nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten()
+            # ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs')
+            # ax.plot(nX, nY, 'r-')
 
-            tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
-            tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
-            plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
-            plt.plot(tX, tY, 'r-')
+            # nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten()
+            # #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs')
+            # ax.plot(nX, nY, 'g-')
 
-            nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
-            nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
-            #plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
-            plt.plot(nX, nY, 'g-')
-            plt.axis('equal')
+            # tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten()
+            # tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten()
+            # ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^')
+            # ax.plot(tX, tY, 'r-')
+
+            # nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten()
+            # nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten()
+            # #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^')
+            # ax.plot(nX, nY, 'g-')
 
         elif self.dim == 3:
-            fig = plt.figure(3)
-            fig.clf()
-            ax = fig.add_subplot(111, projection='3d')
             X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten()
             Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten()
             Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten()
@@ -630,16 +634,50 @@ class CurvView(object):
             Y = np.r_[Y1, Y2, Y3]
             Z = np.r_[Z1, Z2, Z3]
 
-            plt.plot(X, Y, 'b', zs=Z)
+            ax.plot(X, Y, 'b', zs=Z)
             ax.set_zlabel('x3')
 
         ax.grid(True)
-        ax.hold(False)
         ax.set_xlabel('x1')
         ax.set_ylabel('x2')
 
         if showIt: plt.show()
 
+    def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None):
+        if self.dim == 3: raise NotImplementedError('This is not yet done!')
+
+        import matplotlib.pyplot as plt
+        import matplotlib
+        from mpl_toolkits.mplot3d import Axes3D
+        import matplotlib.colors as colors
+        import matplotlib.cm as cmx
+
+        if ax is None: ax = plt.subplot(111)
+        jet = cm = plt.get_cmap('jet')
+        cNorm  = colors.Normalize(
+            vmin=I.min() if clim is None else clim[0],
+            vmax=I.max() if clim is None else clim[1])
+
+        scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet)
+        # ax.set_xlim((self.x0[0], self.h[0].sum()))
+        # ax.set_ylim((self.x0[1], self.h[1].sum()))
+
+        Nx = self.r(self.gridN[:,0],'N','N','M')
+        Ny = self.r(self.gridN[:,1],'N','N','M')
+        cell = self.r(I,'CC','CC','M')
+
+        for ii in range(self.nCx):
+            for jj in range(self.nCy):
+                I = [ii,ii+1,ii+1,ii]
+                J = [jj,jj,jj+1,jj+1]
+                ax.add_patch(plt.Polygon(np.c_[Nx[I,J],Ny[I,J]], facecolor=scalarMap.to_rgba(cell[ii,jj]), edgecolor='k' if grid else 'none'))
+
+        scalarMap._A = []  # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots
+        ax.set_xlabel('x')
+        ax.set_ylabel('y')
+        if showIt: plt.show()
+        return [scalarMap]
+
 
 if __name__ == '__main__':
     from SimPEG import *

SimPEG/Optimization.py

@@ -1008,4 +1008,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
             indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
             delx[indx] = 0.
 
-            return delx
+            return delx

SimPEG/PropMaps.py

@@ -74,7 +74,7 @@ class Property(object):
                 if linkedMap is None:
                     return None
                 linkMap = linkMapClass(None) * linkedMap
-                m = getattr(self, '%s'%linkName)
+                m = getattr(self, '%sModel'%linkName)
                 return linkMap.deriv( m )
 
             m = getattr(self, '%sModel'%prop.name)
@@ -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

@@ -1,4 +1,6 @@
-import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
+import Utils, Maps, Mesh
+import numpy as np
+import scipy.sparse as sp
 
 class RegularizationMesh(object):
     """
@@ -311,6 +313,9 @@ class BaseRegularization(object):
             tmp = indActive
             indActive = np.zeros(mesh.nC, dtype=bool)
             indActive[tmp] = True
+        if indActive is not None and mapping is None:
+            mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
+
         self.regmesh = RegularizationMesh(mesh,indActive)
         self.mapping = mapping or self.mapPair(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
@@ -400,7 +405,238 @@ class BaseRegularization(object):
 
         return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
 
-class Tikhonov(BaseRegularization):
+class Simple(BaseRegularization):
+    """
+    Simple regularization that does not include length scales in the derivatives.
+    """
+
+    mrefInSmooth = False  #: include mref in the smoothness?
+    alpha_s      = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
+    alpha_x      = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'],     "Weight for the first derivative in the x direction")
+    alpha_y      = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'],     "Weight for the first derivative in the y direction")
+    alpha_z      = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'],     "Weight for the first derivative in the z direction")
+    cell_weights = 1.
+
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+
+        if isinstance(self.cell_weights,float):
+            self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
+
+    @property
+    def Wsmall(self):
+        """Regularization matrix Wsmall"""
+        if getattr(self,'_Wsmall', None) is None:
+            self._Wsmall = Utils.sdiag((self.alpha_s*self.cell_weights)**0.5)
+        return self._Wsmall
+
+    @property
+    def Wx(self):
+        """Regularization matrix Wx"""
+        if getattr(self, '_Wx', None) is None:
+            self._Wx = Utils.sdiag((self.alpha_x * (self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.regmesh.cellDiffxStencil
+        return self._Wx
+
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
+        if getattr(self, '_Wy', None) is None:
+            self._Wy = Utils.sdiag((self.alpha_y * (self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.regmesh.cellDiffyStencil
+        return self._Wy
+
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
+        if getattr(self, '_Wz', None) is None:
+            self._Wz = Utils.sdiag((self.alpha_z * (self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.regmesh.cellDiffzStencil
+        return self._Wz
+
+#    @property
+#    def Wsmooth(self):
+#        """Full smoothness regularization matrix W"""
+#        print 'wtf why are we using Wsmooth'
+#        raise NotImplementedError
+#        if getattr(self, '_Wsmooth', None) is None:
+#            wlist = (self.Wx,)
+#            if self.regmesh.dim > 1:
+#                wlist += (self.Wy,)
+#            if self.regmesh.dim > 2:
+#                wlist += (self.Wz,)
+#            self._Wsmooth = sp.vstack(wlist)
+#        return self._Wsmooth
+#
+#    @property
+#    def W(self):
+#        """Full regularization matrix W"""
+#        print 'wtf why are we using W'
+#        if getattr(self, '_W', None) is None:
+#            wlist = (self.Wsmall, self.Wx)
+#            if self.regmesh.dim > 1:
+#                wlist += (self.Wy,)
+#            if self.regmesh.dim > 2:
+#                wlist += (self.Wz,)
+#            self._W = sp.vstack(wlist)
+#        return self._W
+
+
+    @Utils.timeIt
+    def _evalSmall(self, m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmallDeriv(self, m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
+        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
+
+    @Utils.timeIt
+    def _evalSmall2Deriv(self, m, v = None):
+        rDeriv = self.Wsmall * ( self.mapping.deriv(m - self.mref) )
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothx(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wx * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wx * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothy(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wy * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wy * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothz(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wz * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wz * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmooth(self, m):
+        phiSmooth = self._evalSmoothx(m)
+        if self.regmesh.dim > 1:
+            phiSmooth += self._evalSmoothy(m)
+        if self.regmesh.dim > 2:
+            phiSmooth += self._evalSmoothz(m)
+        return phiSmooth
+
+    @Utils.timeIt
+    def _evalSmoothxDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wx * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wx * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wx * ( self.mapping * m )
+            return r.T * ( self.Wx * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothx2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wx * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wx * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothyDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wy * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wy * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wy * ( self.mapping * m )
+            return r.T * ( self.Wy * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothy2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wy * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wy * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothzDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wz * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wz * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wz * ( self.mapping * m )
+            return r.T * ( self.Wz * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothz2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wz * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wz * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothDeriv(self, m):
+        deriv = self._evalSmoothxDeriv(m)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyDeriv(m)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzDeriv(m)
+        return deriv
+
+    @Utils.timeIt
+    def _evalSmooth2Deriv(self, m, v=None):
+        deriv = self._evalSmoothx2Deriv(m, v)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothy2Deriv(m, v)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothz2Deriv(m, v)
+        return deriv
+
+
+    @Utils.timeIt
+    def eval(self, m):
+        return self._evalSmall(m) + self._evalSmooth(m)
+
+    @Utils.timeIt
+    def evalDeriv(self, m):
+        """
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+
+    @Utils.timeIt
+    def eval2Deriv(self, m, v=None):
+        return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v)
+
+
+
+class Tikhonov(Simple):
     """
     L2 Tikhonov regularization with both smallness and smoothness (first order
     derivative) contributions.
@@ -490,56 +726,131 @@ class Tikhonov(BaseRegularization):
             self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
         return self._Wzz
 
+
     @property
-    def Wsmooth(self):
+    def Wsmooth2(self):
         """Full smoothness regularization matrix W"""
         if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx, self.Wxx)
+            wlist = (self.Wxx)
             if self.regmesh.dim > 1:
-                wlist += (self.Wy, self.Wyy)
+                wlist += (self.Wyy)
             if self.regmesh.dim > 2:
-                wlist += (self.Wz, self.Wzz)
+                wlist += (self.Wzz)
             self._Wsmooth = sp.vstack(wlist)
         return self._Wsmooth
 
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
     @Utils.timeIt
-    def _evalSmall(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
+    def _evalSmoothxx(self, m):
         if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
+            r = self.Wxx * ( self.mapping * (m - self.mref) )
         elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * (m) )
+            r = self.Wxx * ( self.mapping * (m) )
         return 0.5 * r.dot(r)
 
+    @Utils.timeIt
+    def _evalSmoothyy(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wyy * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wyy * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothzz(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wzz * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wzz * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmooth2(self, m):
+        phiSmooth2 = self._evalSmoothxx(m)
+        if self.regmesh.dim > 1:
+            phiSmooth2 += self._evalSmoothyy(m)
+        if self.regmesh.dim > 2:
+            phiSmooth2 += self._evalSmoothzz(m)
+        return phiSmooth2
+
+    @Utils.timeIt
+    def _evalSmoothxxDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wxx * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wxx * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wxx * ( self.mapping * m )
+            return r.T * ( self.Wxx * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothyyDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wyy * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wyy * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wyy * ( self.mapping * m )
+            return r.T * ( self.Wyy * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothzzDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wzz * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wzz * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wzz * ( self.mapping * m )
+            return r.T * ( self.Wzz * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothxx2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wxx * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wxx * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothyy2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wyy * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wyy * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothzz2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wzz * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wzz * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothDeriv2(self, m):
+        deriv = self._evalSmoothxxDeriv(m)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyyDeriv(m)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzzDeriv(m)
+        return deriv
+
+    @Utils.timeIt
+    def _evalSmooth2Deriv2(self, m, v=None):
+        deriv = self._evalSmoothxx2Deriv(m, v)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyy2Deriv(m, v)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzz2Deriv(m, v)
+        return deriv
+
+
     @Utils.timeIt
     def eval(self, m):
-        return self._evalSmall(m) + self._evalSmooth(m)
-
-    @Utils.timeIt
-    def _evalSmallDeriv(self,m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
-
-    @Utils.timeIt
-    def _evalSmoothDeriv(self,m):
-        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * ( m - self.mref ) )
-            return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m )
-            return r.T * ( self.Wsmooth * self.mapping.deriv(m) )
+        return self._evalSmall(m) + self._evalSmooth(m) + self._evalSmooth2(m)
 
     @Utils.timeIt
     def evalDeriv(self, m):
@@ -557,185 +868,135 @@ class Tikhonov(BaseRegularization):
             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
-        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m)
+
+    def eval2Deriv(self, m, v=None):
+        """
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + self._evalSmooth2Deriv2(m, v)
 
 
-class Simple(Tikhonov):
+
+class Sparse(Simple):
     """
-    Simple regularization that does not include length scales in the derivatives.
+        The regularization is:
+    
+        .. math::
+    
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R W(m-m_\\text{ref})}
+    
+        where the IRLS weight
+    
+        .. math::
+    
+            R = \eta TO FINISH LATER!!!
+    
+        So the derivative is straight forward:
+    
+        .. math::
+    
+            R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})}
+    
+        The IRLS weights are recomputed after each beta solves.
+        It is strongly recommended to do a few Gauss-Newton iterations
+        before updating.
     """
-
-    mrefInSmooth = False  #: SMOOTH and SMOOTH_MOD_DIF options
-    alpha_s     = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
-    alpha_x     = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
-    alpha_y     = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
-    alpha_z     = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
-    wght        = 1.
+        
+    # set default values
+    eps_p = 1e-1        # Threshold value for the model norm
+    eps_q = 1e-1        # Threshold value for the model gradient norm
+    curModel = None     # Requires model to compute the weights
+    l2model = None
+    gamma = 1.          # Model norm scaling to smooth out convergence
+    norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz)
+    cell_weights = 1.        # Consider overwriting with sensitivity weights
 
     def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
 
-        if isinstance(self.wght,float):
-            self.wght = np.ones(self.regmesh.nC) * self.wght
+        if isinstance(self.cell_weights,float):
+            self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
 
     @property
     def Wsmall(self):
         """Regularization matrix Wsmall"""
         if getattr(self,'_Wsmall', None) is None:
-            self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
+            if getattr(self, 'curModel', None) is None:
+                self.Rs = Utils.speye(self.regmesh.nC)
+
+            else:
+                f_m = self.mapping * (self.curModel - self.reg.mref)
+                self.rs = self.R(f_m , self.eps_p, self.norms[0])
+                self.Rs = Utils.sdiag( self.rs )
+
+            self._Wsmall = Utils.sdiag((self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs
+
         return self._Wsmall
 
     @property
     def Wx(self):
         """Regularization matrix Wx"""
-        if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil
+        if getattr(self,'_Wx', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel)
+                self.rx = self.R( f_m , self.eps_q, self.norms[1])
+                self.Rx = Utils.sdiag( self.rx )
+
+            self._Wx = Utils.sdiag(( self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
+
         return self._Wx
 
     @property
     def Wy(self):
         """Regularization matrix Wy"""
-        if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil
+        if getattr(self,'_Wy', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel)
+                self.ry = self.R( f_m , self.eps_q, self.norms[2])
+                self.Ry = Utils.sdiag( self.ry )
+
+            self._Wy = Utils.sdiag((self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
+
         return self._Wy
 
     @property
     def Wz(self):
         """Regularization matrix Wz"""
-        if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil
+        if getattr(self,'_Wz', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel)
+                self.rz = self.R( f_m , self.eps_q, self.norms[3])
+                self.Rz = Utils.sdiag( self.rz )
+
+            self._Wz = Utils.sdiag((self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
+
         return self._Wz
 
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx,)
-            if self.regmesh.dim > 1:
-                wlist += (self.Wy,)
-            if self.regmesh.dim > 2:
-                wlist += (self.Wz,)
-            self._Wsmooth = sp.vstack(wlist)
-        return self._Wsmooth
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
-    @Utils.timeIt
-    def _evalSmall(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m)
-        return 0.5 * r.dot(r)
-
-
-class Sparse(Simple):
-
-    # set default values
-    eps_p      = 1e-1
-    eps_q      = 1e-1
-    curModel = None # use a model to compute the weights
-    gamma    = 1.
-    norms    = [0., 2., 2., 2.]
-    wght     = 1.
-
-    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
-
-        if isinstance(self.wght,float):
-            self.wght = np.ones(self.regmesh.nC) * self.wght
-
-    @property
-    def Wsmall(self):
-        """Regularization matrix Wsmall"""
-        if getattr(self, 'curModel', None) is None:
-            self.Rs = Utils.speye(self.regmesh.nC)
-
-        else:
-            f_m = self.curModel - self.reg.mref
-            self.rs = self.R(f_m , self.eps_p, self.norms[0])
-            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
-            self.Rs = Utils.sdiag( self.rs )
-
-        return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
-
-
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffxStencil * self.curModel
-            self.rx = self.R( f_m , self.eps_q, self.norms[1])
-            self.Rx = Utils.sdiag( self.rx )
-
-        return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
-
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffyStencil * self.curModel
-            self.ry = self.R( f_m , self.eps_q, self.norms[2])
-            self.Ry = Utils.sdiag( self.ry )
-
-        return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
-
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffzStencil * self.curModel
-            self.rz = self.R( f_m , self.eps_q, self.norms[3])
-            self.Rz = Utils.sdiag( self.rz )
-
-        return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
-
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        #if getattr(self, '_Wsmooth', None) is None:
-        wlist = (self.Wx,)
-        if self.regmesh.dim > 1:
-            wlist += (self.Wy,)
-        if self.regmesh.dim > 2:
-            wlist += (self.Wz,)
-        #self._Wsmooth = sp.vstack(wlist)
-        return sp.vstack(wlist)
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        #if getattr(self, '_W', None) is None:
-        wlist = (self.Wsmall, self.Wsmooth)
-        #self._W = sp.vstack(wlist)
-        return sp.vstack(wlist)
-
     def R(self, f_m , eps, exponent):
 
-        eta = (eps**(1-exponent/2.))**0.5
-        r = eta / (f_m**2.+ eps**2.)**((1-exponent/2.)/2.)
+        # Eta scaling is important for mix-norms...do not mess with it
+        eta = (eps**(1.-exponent/2.))**0.5
+        r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)
 
         return r

SimPEG/Utils/__init__.py

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

SimPEG/Utils/io_utils.py

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

SimPEG/Utils/modelutils.py

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

docs/examples/DC_Forward_PseudoSection.rst

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

docs/examples/Mesh_Basic_ForwardDC.rst

@@ -0,0 +1,25 @@
+.. _examples_Mesh_Basic_ForwardDC:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+Mesh: Basic Forward 2D DC Resistivity
+=====================================
+
+2D DC forward modeling example with Tensor and Curvilinear Meshes
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.Mesh_Basic_ForwardDC.run()
+
+.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_ForwardDC.py
+    :language: python
+    :linenos:

docs/examples/Utils_surface2ind_topo.rst

@@ -1,4 +1,4 @@
-.. _examples_Forward_BasicDirectCurrent:
+.. _examples_Utils_surface2ind_topo:
 
 .. --------------------------------- ..
 ..                                   ..
@@ -8,14 +8,17 @@
 ..                                   ..
 .. --------------------------------- ..
 
-Forward BasicDirectCurrent
-==========================
+
+Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below
+a topographic surface.
+
+
 
 .. plot::
 
     from SimPEG import Examples
-    Examples.Forward_BasicDirectCurrent.run()
+    Examples.Utils_surface2ind_topo.run()
 
-.. literalinclude:: ../../SimPEG/Examples/Forward_BasicDirectCurrent.py
+.. literalinclude:: ../../SimPEG/Examples/Utils_surface2ind_topo.py
     :language: python
     :linenos:

setup.py

@@ -5,16 +5,17 @@ SimPEG is a python package for simulation and gradient based
 parameter estimation in the context of geophysical applications.
 """
 
-import numpy as np
-
 import os
 import sys
 import subprocess
 
 from distutils.core import setup
+from distutils.command.build_ext import build_ext
 from setuptools import find_packages
 from distutils.extension import Extension
 
+
+
 CLASSIFIERS = [
 'Development Status :: 4 - Beta',
 'Intended Audience :: Developers',
@@ -51,11 +52,16 @@ if args.count("build_ext") > 0 and args.count("--inplace") == 0:
 try:
     from Cython.Build import cythonize
     from Cython.Distutils import build_ext
-    cythonKwargs = dict(cmdclass={'build_ext': build_ext})
     USE_CYTHON = True
 except Exception, e:
     USE_CYTHON = False
-    cythonKwargs = dict()
+
+class NumpyBuild(build_ext):
+    def finalize_options(self):
+        build_ext.finalize_options(self)
+        __builtins__.__NUMPY_SETUP__ = False
+        import numpy
+        self.include_dirs.append(numpy.get_include())
 
 ext = '.pyx' if USE_CYTHON else '.c'
 
@@ -94,8 +100,8 @@ setup(
     classifiers=CLASSIFIERS,
     platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"],
     use_2to3 = False,
-    include_dirs=[np.get_include()],
+    cmdclass={'build_ext':NumpyBuild},
+    setup_requires=['numpy'],
     ext_modules = extensions,
     scripts=scripts,
-    **cythonKwargs
 )

tests/base/test_PropMaps.py

@@ -1,6 +1,7 @@
 import unittest
 from SimPEG import *
 from scipy.constants import mu_0
+from SimPEG import Tests
 
 
 class MyPropMap(Maps.PropMap):
@@ -187,6 +188,34 @@ 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_regularization.py

@@ -65,10 +65,8 @@ class RegularizationTests(unittest.TestCase):
                     elif mesh.dim == 3:
                         indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5)
 
-                    mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size)
-
                     for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers
-                        reg = r(mesh, mapping=mapping, indActive=indAct)
+                        reg = r(mesh, indActive=indAct)
                         m = np.random.rand(mesh.nC)[indAct]
                         reg.mref = np.ones_like(m)*np.mean(m)
 

tests/em/fdem/forward/test_FDEM_analytics.py

@@ -28,12 +28,12 @@ class FDEM_analyticTests(unittest.TestCase):
 
         x = np.linspace(-10,10,5)
         XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
-        rxList = EM.FDEM.Rx(XYZ, 'exi')
+        rxList = EM.FDEM.Rx.Point_e(XYZ, orientation='x', component='imag')
         Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
 
         survey = EM.FDEM.Survey([Src0])
 
-        prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
+        prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping)
         prb.pair(survey)
 
         try:
@@ -125,8 +125,8 @@ class FDEM_analyticTests(unittest.TestCase):
 
         mapping = [('sigma', Maps.IdentityMap(mesh)),('mu', Maps.IdentityMap(mesh))]
 
-        prbe = EM.FDEM.Problem_h(mesh, mapping=mapping)
-        prbm = EM.FDEM.Problem_e(mesh, mapping=mapping)
+        prbe = EM.FDEM.Problem3D_h(mesh, mapping=mapping)
+        prbm = EM.FDEM.Problem3D_e(mesh, mapping=mapping)
 
         prbe.pair(surveye) # pair problem and survey
         prbm.pair(surveym)

tests/em/fdem/forward/test_FDEM_forwardHB.py

@@ -12,7 +12,7 @@ testBH = True
 verbose = False
 
 TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.)
-#TODO: choose better testing parameters to lower this 
+#TODO: choose better testing parameters to lower this
 
 SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop']
 
@@ -125,4 +125,4 @@ class FDEM_CrossCheck(unittest.TestCase):
             self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB))
 
 if __name__ == '__main__':
-    unittest.main()
+    unittest.main()

tests/em/static/test_DC_2D_analytic.py

@@ -18,9 +18,9 @@ class DCProblemAnalyticTests(unittest.TestCase):
         A0loc = np.r_[-150, 0.]
         A1loc = np.r_[-130, 0.]
         rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]]
-        data_anal = EM.Analytics.DCAnalyticHalf(np.r_[A0loc, 0.], rxloc, sighalf, flag="halfspace")
+        data_anal = EM.Analytics.DCAnalyticHalf(np.r_[A0loc, 0.], rxloc, sighalf, earth_type="halfspace")
 
-        rx = DC.Rx.Dipole(M, N)
+        rx = DC.Rx.Dipole_ky(M, N)
         src0 = DC.Src.Pole([rx], A0loc)
         survey = DC.Survey_ky([src0])
 

tests/em/static/test_DC_analytic.py

@@ -19,8 +19,8 @@ class DCProblemAnalyticTests(unittest.TestCase):
         Bloc = np.r_[200., 0., 0.]
         M = Utils.ndgrid(x-25.,y, np.r_[0.])
         N = Utils.ndgrid(x+25.,y, np.r_[0.])
-        phiA = EM.Analytics.DCAnalyticHalf(Aloc, [M,N], 1e-2, flag="halfspace")
-        phiB = EM.Analytics.DCAnalyticHalf(Bloc, [M,N], 1e-2, flag="halfspace")
+        phiA = EM.Analytics.DCAnalyticHalf(Aloc, [M,N], 1e-2, earth_type="halfspace")
+        phiB = EM.Analytics.DCAnalyticHalf(Bloc, [M,N], 1e-2, earth_type="halfspace")
         data_anal = phiA-phiB
 
         rx = DC.Rx.Dipole(M, N)

tests/mesh/test_cylMesh.py

@@ -146,6 +146,20 @@ class TestCyl2DMesh(unittest.TestCase):
 
         assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
 
+    def test_getInterpMatCartMesh_Cells2Nodes(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        mc = np.arange(Mc.nC)
+        xr = np.linspace(0,0.4,50)
+        xc = np.linspace(0,0.4,50) + 0.2
+        Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'N')
+        Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC')
+        Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC', locTypeTo='N')
+
+        assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3
+
     def test_getInterpMatCartMesh_Faces(self):
 
         Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
@@ -177,6 +191,37 @@ class TestCyl2DMesh(unittest.TestCase):
         assert np.abs(mag[dist > 0.1].min() - 1) < TOL
 
 
+    def test_getInterpMatCartMesh_Faces2Edges(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pf2e = Mc.getInterpolationMatCartMesh(Mr, 'F', locTypeTo='E')
+        mf = np.ones(Mc.nF)
+
+        ecart = Pf2e * mf
+
+        excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey')
+        ezcc = Mr.r(ecart, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 1) < TOL
+        assert np.abs(float(excc[indY]) - 0) < TOL
+        assert np.abs(float(eycc[indX]) - 0) < TOL
+        assert np.abs(float(eycc[indY]) - 1) < TOL
+        assert np.abs((ezcc - 1).sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
     def test_getInterpMatCartMesh_Edges(self):
 
         Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
@@ -185,11 +230,42 @@ class TestCyl2DMesh(unittest.TestCase):
         Pe = Mc.getInterpolationMatCartMesh(Mr, 'E')
         me = np.ones(Mc.nE)
 
-        erect = Pe * me
+        ecart = Pe * me
 
-        excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex')
-        eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey')
-        ezcc = Mr.r(erect, 'E', 'Ez')
+        excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex')
+        eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey')
+        ezcc = Mr.aveEz2CC*Mr.r(ecart, 'E', 'Ez')
+
+        indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
+        indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])
+
+        TOL = 1e-2
+        assert np.abs(float(excc[indX]) - 0) < TOL
+        assert np.abs(float(excc[indY]) + 1) < TOL
+        assert np.abs(float(eycc[indX]) - 1) < TOL
+        assert np.abs(float(eycc[indY]) - 0) < TOL
+        assert np.abs(ezcc.sum()) < TOL
+
+        mag = (excc**2 + eycc**2)**0.5
+        dist = ((Mr.gridCC[:,0] + 0.2)**2  + (Mr.gridCC[:,1] + 0.2)**2)**0.5
+
+        assert np.abs(mag[dist > 0.1].max() - 1) < TOL
+        assert np.abs(mag[dist > 0.1].min() - 1) < TOL
+
+
+    def test_getInterpMatCartMesh_Edges2Faces(self):
+
+        Mr = Mesh.TensorMesh([100,100,2], x0='CC0')
+        Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0])
+
+        Pe2f = Mc.getInterpolationMatCartMesh(Mr, 'E', locTypeTo='F')
+        me = np.ones(Mc.nE)
+
+        frect = Pe2f * me
+
+        excc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx')
+        eycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy')
+        ezcc = Mr.r(frect, 'F', 'Fz')
 
         indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5])
         indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5])

tests/mt/test_Problem3D_againstAnalytic.py

@@ -242,9 +242,6 @@ class TestAnalytics(unittest.TestCase):
     def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3))
     def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False))
 
-    # Do a derivative test
-    def test_derivProj1(self):self.assertTrue(DerivProjfieldsTest(halfSpace(1e-2)))
-
     # Do a derivative test of Jvec
     # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
     # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))