minor fix

- make PTv out from the loop
- do not use Data class in the loop for Jv

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

SimPEG/Directives.py

@@ -144,6 +144,7 @@ class BetaSchedule(InversionDirective):
             if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
             self.invProb.beta /= self.coolingFactor
 
+
 class TargetMisfit(InversionDirective):
 
     chifact = 1.
@@ -242,12 +243,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         # Save the file as a npz
         np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
 
-
-# class UpdateReferenceModel(Parameter):
-
-#     mref0 = None
-
-#     def nextIter(self):
 #         mref = getattr(self, 'm_prev', None)
 #         if mref is None:
 #             if self.debug: print 'UpdateReferenceModel is using mref0'
@@ -258,56 +253,138 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 class Update_IRLS(InversionDirective):
 
     eps_min = None
+    eps_p = None
+    eps_q = 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
 
+    # Solving parameter for IRLS (mode:2)
+    IRLSiter   = 0
+    minGNiter = 5
+    maxIRLSiter = 10
+    iterStart = 0
+
+    # Beta schedule
+    coolingFactor = 2.
+    coolingRate = 1
+
+    mode = 1
+
+    @property
+    def target(self):
+        if getattr(self, '_target', None) is None:
+            self._target = self.survey.nD*0.5
+        return self._target
+    @target.setter
+    def target(self, val):
+        self._target = val
 
     def initialize(self):
 
-        # Scale the regularization for changes in norm
-        if getattr(self, 'phi_m_last', None) is not None:
-
-            self.reg.curModel = self.invProb.curModel
-            self.reg.gamma = 1.
-            phim_new = self.reg.eval(self.invProb.curModel)
-            self.gamma = self.phi_m_last / phim_new
-
-        self.reg.curModel = self.invProb.curModel
-        self.reg.gamma = self.gamma
-
-        if getattr(self, 'phi_d_last', None) is None:
-            self.phi_d_last = self.invProb.phi_d
+        if self.mode == 1:
+            self.reg.norms = [2., 2., 2., 2.]
 
     def endIter(self):
-        # Cool the threshold parameter if required
-        if getattr(self, 'factor', None) is not None:
-            eps = self.reg.eps / self.factor
 
-            if getattr(self, 'eps_min', None) is not None:
-                self.reg.eps = np.max([self.eps_min,eps])
+        # After reaching target misfit with l2-norm, switch to IRLS (mode:2)
+        if self.invProb.phi_d < self.target and self.mode == 1:
+            print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
+
+            self.mode = 2
+            print self.eps_p, self.eps_q, self.norms
+            self.reg.eps_p = self.eps_p
+            self.reg.eps_q = self.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
 
-         # Update the model used for the IRLS weights
-        self.reg.curModel = self.invProb.curModel
+                if getattr(self, 'eps_min', None) is not None:
+                    self.reg.eps = np.max([self.eps_min,eps])
+                else:
+                    self.reg.eps = eps
 
-        # Temporarely set gamma to 1. to get raw phi_m
-        self.reg.gamma = 1.
+            # Get phi_m at the end of current iteration
+            self.phi_m_last = self.invProb.phi_m_last
 
-        # Compute new model objective function value
-        phim_new = self.reg.eval(self.invProb.curModel)
+            # 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
 
-        # Update gamma to scale the regularization between IRLS iterations
-        self.reg.gamma = self.phi_m_last / phim_new
+             # Update the model used for the IRLS weights
+            self.reg.curModel = self.invProb.curModel
 
-        # Set the weighting matrix to None so that it is recomputed next time
-        # it is called in the inversion
-        self.reg._W = None
+            # Temporarely set gamma to 1. to get raw phi_m
+            self.reg.gamma = 1.
+
+            # Compute new model objective function value
+            phim_new = self.reg.eval(self.invProb.curModel)
+
+            # Update gamma to scale the regularization between IRLS iterations
+            self.reg.gamma = self.phi_m_last / phim_new
+
+            # Reset the regularization matrices again for new gamma
+            self.reg._Wsmall = None
+            self.reg._Wx = None
+            self.reg._Wy = None
+            self.reg._Wz = None
+
+            # Check if misfit is within the tolerance, otherwise scale beta
+            val = self.invProb.phi_d / (self.survey.nD*0.5)
+
+            if np.abs(1.-val) > self.beta_tol:
+                self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
 
 class Update_lin_PreCond(InversionDirective):
     """
@@ -360,19 +437,3 @@ class Update_Wj(InversionDirective):
             JtJdiag = JtJdiag / max(JtJdiag)
 
             self.reg.wght = JtJdiag
-
-class Scale_Beta(InversionDirective):
-    """
-        Instead of a linear cooling schedule, beta is allowed to change based
-        on the ratio between the target misfit and the current data misfit. The
-        update is done only if the misfit is outside some threshold bounds.
-    """
-    tol = 0.05
-
-    def endIter(self):
-
-        # Check if misfit is within the tolerance, otherwise adjust beta
-        val = self.invProb.phi_d / (self.survey.nD*0.5)
-
-        if np.abs(1.-val) > self.tol:
-            self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d

SimPEG/EM/Base.py

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

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -60,6 +60,20 @@ class Fields(SimPEG.Problem.Fields):
 
         return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
 
+    def _bSecondary(self, solution, srcList):
+        """
+        Total magnetic flux density is sum of primary and secondary
+
+        :param numpy.ndarray solution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: total magnetic flux density
+        """
+        if getattr(self, '_bSecondary', None) is None:
+            raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
+
+        return  self._bSecondary(solution, srcList)
+
     def _h(self, solution, srcList):
         """
         Total magnetic field is sum of primary and secondary
@@ -124,6 +138,21 @@ class Fields(SimPEG.Problem.Fields):
             return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
         return np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex)
 
+    def _bSecondaryDeriv(self, src, du_dm_v, v, adjoint = False):
+        """
+        Total derivative of b with respect to the inversion model. Returns :math:`d\mathbf{b}/d\mathbf{m}` for forward and (:math:`d\mathbf{b}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
+
+        :param Src src: sorce
+        :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
+        :param numpy.ndarray v: vector to take sensitivity product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: derivative times a vector (or tuple for adjoint)
+        """
+        # TODO: modify when primary field is dependent on m
+
+        return self._bDeriv(src, du_dm_v, v, adjoint = adjoint)
+
     def _hDeriv(self, src, du_dm_v, v, adjoint = False):
         """
         Total derivative of h with respect to the inversion model. Returns :math:`d\mathbf{h}/d\mathbf{m}` for forward and (:math:`d\mathbf{h}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
@@ -471,6 +500,8 @@ class Fields3D_b(Fields):
             return 'E'
         elif fieldType == 'b':
             return 'F'
+        elif fieldType == 'bSecondary':
+            return 'F'
         elif (fieldType == 'h') or (fieldType == 'j'):
             return'CCV'
         else:

SimPEG/EM/FDEM/RxFDEM.py

@@ -97,6 +97,19 @@ class Point_b(BaseRx):
         self.projField = 'b'
         super(Point_b, self).__init__(locs, orientation, component)
 
+class Point_bSecondary(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 = 'bSecondary'
+        super(Point_bSecondary, self).__init__(locs, orientation, component)
+
 
 class Point_h(BaseRx):
     """

SimPEG/EM/FDEM/SrcFDEM.py

@@ -555,7 +555,7 @@ class CircularLoop(BaseSrc):
             a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
-            srcfct = MagneticDipoleVectorPotential
+            srcfct = MagneticLoopVectorPotential
             ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
             ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
             az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)

SimPEG/EM/Static/DC/ProblemDC.py

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

SimPEG/EM/Static/IP/ProblemIP.py

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

SimPEG/EM/Static/Utils/StaticUtils.py

@@ -315,107 +315,3 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
 
     return SrcList
 
-
-def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
-    """
-        Write UBC GIF DCIP 2D or 3D observation file
-
-        Input:
-        :string     fileName -> including path where the file is written out
-        :DCsurvey   DC survey class object
-        :string     dtype ->  either '2D' | '3D'
-        :string     stype ->  either 'SURFACE' | 'GENERAL'
-
-        Output:
-        :param UBC2D-Data file
-        :return
-
-        Last edit: February 16th, 2016
-
-        @author: dominiquef
-
-    """
-    from SimPEG import mkvc
-
-    assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'"
-    assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
-
-    fid = open(fileName,'w')
-
-
-    if iptype!=0:
-        fid.write('IPTYPE=%i\n'%iptype)
-
-    else:
-        fid.write('! ' + stype + ' FORMAT\n')
-
-    count = 0
-
-    for ii in range(DCsurvey.nSrc):
-
-        tx = np.c_[DCsurvey.srcList[ii].loc]
-
-        rx = DCsurvey.srcList[ii].rxList[0].locs
-
-        nD = DCsurvey.srcList[ii].nD
-
-        M = rx[0]
-        N = rx[1]
-
-        # Adapt source-receiver location for dtype and stype
-        if dtype=='2D':
-
-            if stype == 'SIMPLE':
-
-                #fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
-                A = np.repeat(tx[0,0],M.shape[0],axis=0)
-                B = np.repeat(tx[0,1],M.shape[0],axis=0)
-                M = M[:,0]
-                N = N[:,0]
-
-                np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
-
-
-            else:
-
-                if stype == 'SURFACE':
-
-                    fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
-                    M = M[:,0]
-                    N = N[:,0]
-
-                if stype == 'GENERAL':
-
-                    # Flip sign for z-elevation to depth
-                    tx[2::2,:] = -tx[2::2,:]
-
-                    fid.writelines("%e " % ii for ii in mkvc(tx[::2,:]))
-                    M = M[:,0::2]
-                    N = N[:,0::2]
-
-                    # Flip sign for z-elevation to depth
-                    M[:,1::2] = -M[:,1::2]
-                    N[:,1::2] = -N[:,1::2]
-
-                fid.write('%i\n'% nD)
-                np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
-
-        if dtype=='3D':
-
-            if stype == 'SURFACE':
-
-                fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:]))
-                M = M[:,0:2]
-                N = N[:,0:2]
-
-            if stype == 'GENERAL':
-
-                fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
-
-            fid.write('%i\n'% nD)
-            np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
-            fid.write('\n')
-
-        count += nD
-
-    fid.close()

SimPEG/EM/TDEM/FieldsTDEM.py

@@ -0,0 +1,142 @@
+import numpy as np
+import scipy.sparse as sp
+import SimPEG
+from SimPEG import Utils
+from SimPEG.EM.Utils import omega
+from SimPEG.Utils import Zero, Identity
+
+class Fields(SimPEG.Problem.TimeFields):
+    """
+
+    Fancy Field Storage for a TDEM survey. Only one field type is stored for
+    each problem, the rest are computed. The fields obejct acts like an array and is indexed by
+
+    .. code-block:: python
+
+        f = problem.fields(m)
+        e = f[srcList,'e']
+        b = f[srcList,'b']
+
+    If accessing all sources for a given field, use the :code:`:`
+
+    .. code-block:: python
+
+        f = problem.fields(m)
+        e = f[:,'e']
+        b = f[:,'b']
+
+    The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies)
+    """
+
+    knownFields = {}
+    dtype = float
+
+    def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+        if adjoint is True:
+            return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint)
+        return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v)
+
+    def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+        if adjoint is True:
+            return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
+        return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
+
+class Fields_Derivs(Fields):
+    knownFields = {
+                    'bDeriv': 'F',
+                    'eDeriv': 'E',
+                    'hDeriv': 'E',
+                    'jDeriv': 'F'
+                  }
+
+
+class Fields_b(Fields):
+    """Fancy Field Storage for a TDEM survey."""
+    knownFields = {'bSolution': 'F'}
+    aliasFields = {
+                    'b': ['bSolution', 'F', '_b'],
+                    'e': ['bSolution', 'E', '_e'],
+                  }
+
+    def startup(self):
+        self.MeSigmaI      = self.survey.prob.MeSigmaI
+        self.MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
+        self.edgeCurl      = self.survey.prob.mesh.edgeCurl
+        self.MfMui         = self.survey.prob.MfMui
+
+    def _b(self, bSolution, srcList, tInd):
+        return bSolution
+
+    def _bDeriv_u(self, tInd, src, dun_dm_v, adjoint=False):
+        return Identity()*dun_dm_v
+
+    def _bDeriv_m(self, tInd, src, v, adjoint=False):
+        return Zero()
+
+    # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+    #     if adjoint is True:
+    #         return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
+    #     return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
+
+    def _e(self, bSolution, srcList, tInd):
+        e = self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) )
+        for i, src in enumerate(srcList):
+            _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
+            e[:,i] = e[:,i] - self.MeSigmaI * S_e
+        return e
+
+    def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False):
+        if adjoint is True:
+            return self.MfMui.T * ( self.edgeCurl * ( self.MeSigmaI.T * dun_dm_v ) )
+        return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * dun_dm_v ) )
+
+    def _eDeriv_m(self, tInd, src, v, adjoint = False):
+        _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
+        bSolution = self[[src],'bSolution',tInd]
+
+        _, S_eDeriv = src.evalDeriv(self.survey.prob.times[tInd], self, adjoint=adjoint)
+
+        if adjoint is True:
+            return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution ) ).T * v - S_eDeriv(self.MeSigmaI.T * v)
+
+        return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution)) * v - self.MeSigmaI * S_eDeriv(v)
+
+
+
+class Fields_e(Fields):
+    """Fancy Field Storage for a TDEM survey."""
+    knownFields = {'eSolution': 'E'}
+    aliasFields = {
+                    'e': ['eSolution', 'E', '_e'],
+                    'b': ['eSolution', 'F', '_b'],
+                  }
+
+    def startup(self):
+        self.MeSigmaI      = self.survey.prob.MeSigmaI
+        self.MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
+        self.edgeCurl      = self.survey.prob.mesh.edgeCurl
+        self.MfMui         = self.survey.prob.MfMui
+
+
+    def _e(self, eSolution, srcList, tInd):
+        return eSolution
+
+    def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False):
+        return dun_dm_v
+
+    def _eDeriv_m(self, tInd, src, v, adjoint = False):
+        return Zero()
+
+    def _b(self, eSolution, srcList, tInd):
+        raise NotImplementedError
+
+    def _bDeriv_u(self, tInd, src, dun_dm_v, adjoint=False):
+        raise NotImplementedError
+
+    def _bDeriv_m(self, tInd, src, v, adjoint=False):
+        raise NotImplementedError
+
+    # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+    #     if adjoint is True:
+    #         return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
+    #     return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)

SimPEG/EM/TDEM/SrcTDEM.py

@@ -0,0 +1,246 @@
+import SimPEG
+from SimPEG import np, Utils
+from SimPEG.Utils import Zero, Identity
+from scipy.constants import mu_0
+from SimPEG.EM.Utils import *
+
+####################################################
+# Sources
+####################################################
+
+class BaseWaveform(object):
+
+    def __init__(self, offTime=0., hasInitialFields=False):
+        self.offTime = offTime
+        self.hasInitialFields = hasInitialFields
+
+    def _assertMatchesPair(self, pair):
+        assert (isinstance(self, pair)
+            ), "Waveform object must be an instance of a %s BaseWaveform class."%(pair.__name__)
+
+    def eval(self, time):
+        raise NotImplementedError
+
+    def evalDeriv(self, time):
+        raise NotImplementedError # needed for E-formulation
+
+
+class StepOffWaveform(BaseWaveform):
+
+    def __init__(self, offTime=0.):
+        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
+
+    def eval(self, time):
+        return 0.
+
+
+class RawWaveform(BaseWaveform):
+
+    def __init__(self, offTime=0.):
+        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
+
+    def eval(self, time):
+        raise NotImplementedError('RawWaveform has not been implemented, you should write it!')
+
+
+class TriangularWaveform(BaseWaveform):
+
+    def __init__(self, offTime=0.):
+        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
+
+    def eval(self, time):
+        raise NotImplementedError('TriangularWaveform has not been implemented, you should write it!')
+
+
+
+class BaseSrc(SimPEG.Survey.BaseSrc):
+
+    # rxPair = Rx
+    integrate = True
+    waveformPair = BaseWaveform
+
+    @property
+    def waveform(self):
+        "A waveform instance is not None"
+        return getattr(self, '_waveform', None)
+    @waveform.setter
+    def waveform(self, val):
+        if self.waveform is None:
+            val._assertMatchesPair(self.waveformPair)
+            self._mapping = val
+        else:
+            self._mapping = self.PropMap(val)
+
+
+    def __init__(self, rxList, waveform = StepOffWaveform(), **kwargs):
+        self.waveform = waveform
+        SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs)
+
+
+    def bInitial(self, prob):
+        return Zero()
+
+    def bInitialDeriv(self, prob, v=None, adjoint=False):
+        return Zero()
+
+    def eInitial(self, prob):
+        return Zero()
+
+    def eInitialDeriv(self, prob, v=None, adjoint=False):
+        return Zero()
+
+    def eval(self, prob, time):
+        S_m = self.S_m(prob, time)
+        S_e = self.S_e(prob, time)
+        return S_m, S_e
+
+    def evalDeriv(self, prob, time, v=None, adjoint=False):
+        if v is not None:
+            return self.S_mDeriv(prob, time, v, adjoint), self.S_eDeriv(prob, time, v, adjoint)
+        else:
+            return lambda v: self.S_mDeriv(prob, time, v, adjoint), lambda v: self.S_eDeriv(prob, time, v, adjoint)
+
+    def S_m(self, prob, time):
+        return Zero()
+
+    def S_e(self, prob, time):
+        return Zero()
+
+    def S_mDeriv(self, prob, time, v=None, adjoint=False):
+        return Zero()
+
+    def S_eDeriv(self, prob, time, v=None, adjoint=False):
+        return Zero()
+
+
+class MagDipole(BaseSrc):
+
+    waveform = None
+    loc = None
+    orientation = 'Z'
+    moment = 1.
+    mu = mu_0
+
+    def __init__(self, rxList, **kwargs):
+        assert self.orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
+        self.integrate = False
+        BaseSrc.__init__(self, rxList, **kwargs)
+
+    def _bfromVectorPotential(self, prob):
+        if prob._eqLocs is 'FE':
+            gridX = prob.mesh.gridEx
+            gridY = prob.mesh.gridEy
+            gridZ = prob.mesh.gridEz
+            C = prob.mesh.edgeCurl
+
+        elif prob._eqLocs is 'EF':
+            gridX = prob.mesh.gridFx
+            gridY = prob.mesh.gridFy
+            gridZ = prob.mesh.gridFz
+            C = prob.mesh.edgeCurl.T
+
+
+        if prob.mesh._meshType is 'CYL':
+            if not prob.mesh.isSymmetric:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
+
+        else:
+            srcfct = MagneticDipoleVectorPotential
+            ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
+            ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
+            az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
+            a = np.concatenate((ax, ay, az))
+
+        return C*a
+
+
+    def bInitial(self, prob):
+
+        if self.waveform.hasInitialFields is False:
+            return Zero()
+
+        return self._bfromVectorPotential(prob)
+
+    def eInitial(self, prob):
+
+        if self.waveform.hasInitialFields is False:
+            return Zero()
+
+        b = self.bInitial(prob)
+        MeSigmaI = prob.MeSigmaI
+        MfMui = prob.MfMui
+        C = prob.mesh.edgeCurl
+
+        return MeSigmaI * (C.T * (MfMui * b))
+
+    def eInitialDeriv(self, prob, v=None, adjoint=False):
+
+        if self.waveform.hasInitialFields is False:
+            return Zero()
+
+        b = self.bInitial(prob)
+        MeSigmaIDeriv = prob.MeSigmaIDeriv
+        MfMui = prob.MfMui
+        C = prob.mesh.edgeCurl
+        S_e = self.S_e(prob, prob.t0)
+
+        # S_e doesn't depend on the model
+
+        if adjoint:
+            return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ).T * v
+
+        return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ) * v
+
+
+    def S_m(self, prob, time):
+        if self.waveform.hasInitialFields is False:
+            raise NotImplementedError
+        return Zero()
+
+    def S_e(self, prob, time):
+        if self.waveform.hasInitialFields is False:
+            raise NotImplementedError
+        return Zero()
+
+class CircularLoop(MagDipole):
+
+    waveform = None
+    loc = None
+    orientation = 'Z'
+    radius = None
+    mu = mu_0
+
+    def __init__(self, rxList, **kwargs):
+        assert self.orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
+        self.integrate = False
+        BaseSrc.__init__(self, rxList, **kwargs)
+
+    def _bfromVectorPotential(self, prob):
+        if prob._eqLocs is 'FE':
+            gridX = prob.mesh.gridEx
+            gridY = prob.mesh.gridEy
+            gridZ = prob.mesh.gridEz
+            C = prob.mesh.edgeCurl
+
+        elif prob._eqLocs is 'EF':
+            gridX = prob.mesh.gridFx
+            gridY = prob.mesh.gridFy
+            gridZ = prob.mesh.gridFz
+            C = prob.mesh.edgeCurl.T
+
+
+        if prob.mesh._meshType is 'CYL':
+            if not prob.mesh.isSymmetric:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', radius=self.radius, mu=self.mu)
+
+        else:
+            srcfct = MagneticLoopVectorPotential
+            ax = srcfct(self.loc, gridX, 'x', mu=self.mu, radius=self.radius)
+            ay = srcfct(self.loc, gridY, 'y', mu=self.mu, radius=self.radius)
+            az = srcfct(self.loc, gridZ, 'z', mu=self.mu, radius=self.radius)
+            a = np.concatenate((ax, ay, az))
+
+        return C*a
+

SimPEG/EM/TDEM/SurveyTDEM.py

@@ -1,10 +1,16 @@
-from SimPEG import Utils, Survey, np
-from SimPEG.Survey import BaseSurvey
+import SimPEG
+from SimPEG import np, Utils
+from SimPEG.Utils import Zero, Identity
+from scipy.constants import mu_0
 from SimPEG.EM.Utils import *
-from BaseTDEM import FieldsTDEM
+import SrcTDEM as Src
 
 
-class RxTDEM(Survey.BaseTimeRx):
+####################################################
+# Receivers
+####################################################
+
+class Rx(SimPEG.Survey.BaseTimeRx):
 
     knownRxTypes = {
                     'ex':['e', 'Ex', 'N'],
@@ -21,7 +27,7 @@ class RxTDEM(Survey.BaseTimeRx):
                    }
 
     def __init__(self, locs, times, rxType):
-        Survey.BaseTimeRx.__init__(self, locs, times, rxType)
+        SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType)
 
     @property
     def projField(self):
@@ -56,144 +62,60 @@ class RxTDEM(Survey.BaseTimeRx):
         u_part = Utils.mkvc(u[src, self.projField, :])
         return P*u_part
 
-    def evalDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
+    def evalDeriv(self, src, mesh, timeMesh, v, adjoint=False):
         P = self.getP(mesh, timeMesh)
-
         if not adjoint:
-            return P * Utils.mkvc(v[src, self.projField, :])
+            return P * v #Utils.mkvc(v[src, self.projField+'Deriv', :])
         elif adjoint:
-            return P.T * v[src, self]
+            # dP_dF_T = P.T * v #[src, self]
+            # newshape = (len(dP_dF_T)/timeMesh.nN, timeMesh.nN )
+            return P.T * v #np.reshape(dP_dF_T, newshape, order='F')
 
 
-class SrcTDEM(Survey.BaseSrc):
-    rxPair = RxTDEM
-    radius = None
+####################################################
+# Survey
+####################################################
 
-    def getInitialFields(self, mesh):
-        F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
-        return F0
-
-    def getJs(self, mesh, time):
-        return None
-
-
-class SrcTDEM_VMD_MVP(SrcTDEM):
-
-    def __init__(self,rxList,loc,waveformType="STEPOFF"):
-        self.loc = loc
-        self.waveformType = waveformType
-        SrcTDEM.__init__(self,rxList)
-
-    def getInitialFields(self, mesh):
-        """Vertical magnetic dipole, magnetic vector potential"""
-        if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"
-            if mesh._meshType is 'CYL':
-                if mesh.isSymmetric:
-                    MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
-                else:
-                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            elif mesh._meshType is 'TENSOR':
-                MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
-            else:
-                raise Exception('Unknown mesh for VMD')
-            return {"b": mesh.edgeCurl*MVP}
-        elif self.waveformType == "GENERAL":
-            print ">> General waveform: Zero initial condition"
-            return {"b": np.zeros(mesh.nF)}
-        else:
-            raise NotImplementedError("Only use STEPOFF or GENERAL")
-
-    def getMeS(self, mesh, MfMui):
-        if mesh._meshType is 'CYL':
-            if mesh.isSymmetric:
-                MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
-            else:
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-        elif mesh._meshType is 'TENSOR':
-            MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
-        else:
-            raise Exception('Unknown mesh for VMD')
-        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
-
-
-class SrcTDEM_CircularLoop_MVP(SrcTDEM):
-    def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
-        self.loc = loc
-        self.radius = radius
-        self.waveformType = waveformType
-        SrcTDEM.__init__(self,rxList)
-
-    def getInitialFields(self, mesh):
-        """Circular Loop, magnetic vector potential"""
-        if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"
-            if mesh._meshType is 'CYL':
-                if mesh.isSymmetric:
-                    MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
-                else:
-                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            elif mesh._meshType is 'TENSOR':
-                MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
-            else:
-                raise Exception('Unknown mesh for CircularLoop')
-            return {"b": mesh.edgeCurl*MVP}
-        elif self.waveformType == "GENERAL":
-            print ">> General waveform: Zero initial condition"
-            return {"b": np.zeros(mesh.nF)}
-        else:
-            raise NotImplementedError("Only use STEPOFF or GENERAL")
-
-    def getMeS(self, mesh, MfMui):
-        if mesh._meshType is 'CYL':
-            if mesh.isSymmetric:
-                MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
-            else:
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-        elif mesh._meshType is 'TENSOR':
-            MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
-        else:
-            raise Exception('Unknown mesh for CircularLoop')
-        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
-
-
-class SurveyTDEM(Survey.BaseSurvey):
+class Survey(SimPEG.Survey.BaseSurvey):
     """
-        docstring for SurveyTDEM
+    Time domain electromagnetic survey
     """
-    srcPair = SrcTDEM
+
+    srcPair = Src.BaseSrc
+    rxPair = Rx
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
         self.srcList = srcList
-        Survey.BaseSurvey.__init__(self, **kwargs)
+        SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
 
     def eval(self, u):
-        data = Survey.Data(self)
+        data = SimPEG.Survey.Data(self)
         for src in self.srcList:
             for rx in src.rxList:
                 data[src, rx] = rx.eval(src, self.mesh, self.prob.timeMesh, u)
         return data
 
     def evalDeriv(self, u, v=None, adjoint=False):
-        assert v is not None, 'v to multiply must be provided.'
+        raise Exception('Use Receivers to project fields deriv.')
+        # assert v is not None, 'v to multiply must be provided.'
 
-        if not adjoint:
-            data = Survey.Data(self)
-            for src in self.srcList:
-                for rx in src.rxList:
-                    data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
-            return data
-        else:
-            f = FieldsTDEM(self.mesh, self)
-            for src in self.srcList:
-                for rx in src.rxList:
-                    Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
-                    Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
-                    if rx.projField not in f: # first time we are projecting
-                        f[src, rx.projField, :] = Ptv
-                    else: # there are already fields, so let's add to them!
-                        f[src, rx.projField, :] += Ptv
-            return f
+        # if not adjoint:
+        #     data = SimPEG.Survey.Data(self)
+        #     for src in self.srcList:
+        #         for rx in src.rxList:
+        #             data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
+        #     return data
+        # else:
+        #     f = FieldsTDEM(self.mesh, self)
+        #     for src in self.srcList:
+        #         for rx in src.rxList:
+        #             Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
+        #             Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
+        #             if rx.projField not in f: # first time we are projecting
+        #                 f[src, rx.projField, :] = Ptv
+        #             else: # there are already fields, so let's add to them!
+        #                 f[src, rx.projField, :] += Ptv
+        #     return f
 
 

SimPEG/EM/TDEM/TDEM.py

@@ -0,0 +1,553 @@
+from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
+from SimPEG.EM.Base import BaseEMProblem
+from SimPEG.EM.TDEM.SurveyTDEM import Survey as SurveyTDEM
+from SimPEG.EM.TDEM.FieldsTDEM import *
+from scipy.constants import mu_0
+import time
+
+class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
+    """
+    We start with the first order form of Maxwell's equations
+    """
+    surveyPair = SurveyTDEM
+    fieldsPair = Fields
+
+    def __init__(self, mesh, mapping=None, **kwargs):
+        Problem.BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
+
+    def fields(self, m):
+        """
+        Solve the forward problem for the fields.
+
+        :param numpy.array m: inversion model (nP,)
+        :rtype numpy.array:
+        :return F: fields
+        """
+
+        tic = time.time()
+        self.curModel = m
+
+        F = self.fieldsPair(self.mesh, self.survey)
+
+        # set initial fields
+        F[:,self._fieldType+'Solution',0] = self.getInitialFields()
+
+        # timestep to solve forward
+        if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
+        Ainv = None
+        for tInd, dt in enumerate(self.timeSteps):
+            if Ainv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same
+                Ainv.clean()
+                Ainv = None
+
+            if Ainv is None:
+                A = self.getAdiag(tInd)
+                if self.verbose: print 'Factoring...   (dt = %e)'%dt
+                Ainv = self.Solver(A, **self.solverOpts)
+                if self.verbose: print 'Done'
+
+            rhs = self.getRHS(tInd+1) # this is on the nodes of the time mesh
+            Asubdiag = self.getAsubdiag(tInd)
+
+            if self.verbose: print ('    Solving...   (tInd = %i)')% (tInd+1)
+            sol = Ainv * (rhs - Asubdiag * F[:,self._fieldType+'Solution',tInd]) # taking a step
+
+            if self.verbose: print '    Done...'
+
+            if sol.ndim == 1:
+                sol.shape = (sol.size,1)
+            F[:,self._fieldType+'Solution',tInd+1] = sol
+        if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
+        Ainv.clean()
+        return F
+
+
+    def Jvec(self, m, v, f=None):
+        """
+        Jvec computes the sensitivity times a vector
+
+        .. math::
+            \mathbf{J} \mathbf{v} = \\frac{d\mathbf{P}}{d\mathbf{F}} \left( \\frac{d\mathbf{F}}{d\mathbf{u}} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{\partial\mathbf{F}}{\partial\mathbf{m}} \\right) \mathbf{v}
+
+        where
+
+        .. math::
+            \mathbf{A} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} = \\frac{d \mathbf{RHS}}{d \mathbf{m}}
+        """
+
+        if f is None:
+           f = self.fields(m)
+
+        ftype = self._fieldType + 'Solution' # the thing we solved for
+        self.curModel = m
+
+        # mat to store previous time-step's solution deriv times a vector for each source
+        # size: nu x nSrc
+
+        # this is a bit silly
+
+        # if self._fieldType is 'b' or self._fieldType is 'j':
+        #     ifields = np.zeros((self.mesh.nF, len(Srcs)))
+        # elif self._fieldType is 'e' or self._fieldType is 'h':
+        #     ifields = np.zeros((self.mesh.nE, len(Srcs)))
+
+        # for i, src in enumerate(self.survey.srcList):
+        dun_dm_v = np.hstack([Utils.mkvc(self.getInitialFieldsDeriv(src,v),2) for src in self.survey.srcList]) # can over-write this at each timestep
+        #
+        df_dm_v = Fields_Derivs(self.mesh, self.survey) # store the field derivs we need to project to calc full deriv
+
+        Adiaginv = None
+
+        for tInd, dt in zip(range(self.nT), self.timeSteps):
+            if Adiaginv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same
+                Adiaginv.clean()
+                Adiaginv = None
+
+            if Adiaginv is None:
+                A = self.getAdiag(tInd)
+                Adiaginv = self.Solver(A, **self.solverOpts)
+
+            Asubdiag = self.getAsubdiag(tInd)
+
+            for i, src in enumerate(self.survey.srcList):
+
+                # here, we are lagging by a timestep, so filling in as we go
+                for projField in set([rx.projField for rx in src.rxList]):
+                    # Seogi: df_duFun?
+                    df_dmFun = getattr(f, '_%sDeriv'%projField, None)
+                    # df_dm_v is dense, but we only need the times at (rx.P.T * ones > 0)
+                    # This should be called rx.footprint
+                    df_dm_v[src, '%sDeriv'%projField , tInd] = df_dmFun(tInd, src, dun_dm_v[:,i], v)
+
+                un_src = f[src,ftype,tInd+1]
+
+                dA_dm_v   = self.getAdiagDeriv(tInd, un_src, v) # cell centered on time mesh
+                dRHS_dm_v = self.getRHSDeriv(tInd+1, src, v) # on nodes of time mesh
+
+                dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, f[src,ftype,tInd], v)
+
+                JRHS = dRHS_dm_v - dAsubdiag_dm_v - dA_dm_v
+
+                # step in time and overwrite
+                if tInd != len(self.timeSteps+1):
+                    dun_dm_v[:,i] = Adiaginv * (JRHS - Asubdiag * dun_dm_v[:,i])
+
+        # Seogi: suspcious spot
+        # Jv = self.dataPair(self.survey)
+        Jv = []
+        for src in self.survey.srcList:
+            for rx in src.rxList:
+                # Looping over data class append memory as well!!
+                # Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(df_dm_v[src,'%sDeriv'%rx.projField,:]))
+                Jv.append(rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(df_dm_v[src,'%sDeriv'%rx.projField,:])))
+        Adiaginv.clean()
+        # del df_dm_v, dun_dm_v, Asubdiag
+        # return Utils.mkvc(Jv)
+        return np.hstack(Jv)
+
+    def Jtvec(self, m, v, f=None):
+
+        """
+        Jvec computes the adjoint of the sensitivity times a vector
+
+        .. math::
+            \mathbf{J}^\\top \mathbf{v} =  \left( \\frac{d\mathbf{u}}{d\mathbf{m}} ^ \\top \\frac{d\mathbf{F}}{d\mathbf{u}} ^ \\top  + \\frac{\partial\mathbf{F}}{\partial\mathbf{m}} ^ \\top \\right) \\frac{d\mathbf{P}}{d\mathbf{F}} ^ \\top \mathbf{v}
+
+        where
+
+        .. math::
+            \\frac{d\mathbf{u}}{d\mathbf{m}} ^\\top \mathbf{A}^\\top  + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} ^ \\top = \\frac{d \mathbf{RHS}}{d \mathbf{m}} ^ \\top
+        """
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+        ftype = self._fieldType + 'Solution' # the thing we solved for
+
+        # Ensure v is a data object.
+        if not isinstance(v, self.dataPair):
+            v = self.dataPair(self.survey, v)
+
+        df_duT_v = Fields_Derivs(self.mesh, self.survey)
+        ATinv_df_duT_v = np.zeros((len(self.survey.srcList), len(f[self.survey.srcList[0],ftype,0])), dtype=float) # same size as fields at a single timestep
+
+        JTv = np.zeros(m.shape, dtype=float)
+
+        # Loop over sources and receivers to create a fields object: PT_v, df_duT_v, df_dmT_v
+        PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources)
+        for src in self.survey.srcList:
+            # Looping over initializing field class is appending memory!
+            # PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources)
+            # initialize size
+            df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(f[src, self._fieldType, :])
+
+            for rx in src.rxList:
+                print ('_%sDeriv')%(rx.projField)
+                PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) # this is +=
+
+                # PT_v = np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F')
+                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+
+                for tInd in range(self.nT+1):
+                    cur = df_duTFun(tInd, src, None, Utils.mkvc(PT_v[src,'%sDeriv'%rx.projField,tInd]), adjoint=True)
+                    df_duT_v[src, '%sDeriv'%self._fieldType, tInd] = df_duT_v[src, '%sDeriv'%self._fieldType, tInd] + Utils.mkvc(cur[0],2)
+                    JTv = cur[1] + JTv
+
+        del PT_v # no longer need this
+
+        AdiagTinv = None
+
+        # Do the back-solve through time
+        for tIndP in reversed(range(self.nT + 1)):
+            tInd = tIndP - 1
+            if AdiagTinv is not None and (tInd <= self.nT and self.timeSteps[tInd] != self.timeSteps[tInd+1]): # if the previous timestep is the same --> no need to refactor the matrix
+                AdiagTinv.clean()
+                AdiagTinv = None
+
+            # refactor if we need to
+            if AdiagTinv is None and tInd > -1:
+                Adiag = self.getAdiag(tInd)
+                AdiagTinv = self.Solver(Adiag.T, **self.solverOpts)
+
+            dAsubdiag_dm_v = Zero()
+
+            if tInd < self.nT - 1:
+                Asubdiag = self.getAsubdiag(tInd+1)
+
+
+            for isrc, src in enumerate(self.survey.srcList):
+                # solve against df_duT_v
+                if tInd >= self.nT-1:
+                    # last timestep (first to be solved)
+                    ATinv_df_duT_v[isrc,:] = AdiagTinv * df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]
+                elif tInd > -1:
+                # else:
+                    ATinv_df_duT_v[isrc,:] = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+                else:
+                    # AdiagTinv = I
+                    ATinv_df_duT_v[isrc,:] = Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])
+                    # - Utils.mkvc(Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+                    # (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+
+                if tInd < self.nT - 1:
+                    dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], ATinv_df_duT_v[isrc,:], adjoint = True)
+
+                if tInd > -1:
+                    un_src = f[src,ftype,tInd+1]
+                    dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
+                    dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
+
+                    JTv = JTv + Utils.mkvc(- dAT_dm_v - dAsubdiag_dm_v + dRHST_dm_v)
+                else:
+                    # dA_dm_v = self.getInitialFieldsDeriv(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True)
+                    # print np.linalg.norm(self.getInitialFieldsDeriv(src, df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True))
+                    # print np.linalg.norm(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
+                    # vec = - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]) + Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
+                    # dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint = True)
+                    dRHST_dm_v = Utils.mkvc(self.getInitialFieldsDeriv(src, Utils.mkvc(ATinv_df_duT_v[isrc,:]) , adjoint=True))
+
+                    JTv = JTv + Utils.mkvc( -dAsubdiagT_dm_v + dRHST_dm_v) #
+
+
+
+                    # # dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
+                    # dRHST_dm_v0 = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
+                    # dRHST_dm_v1 = self.getInitialFieldsDeriv( Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
+                    # JTv = JTv + Utils.mkvc(dRHST_dm_v0 + dRHST_dm_v1)
+
+                    # print 'here'
+                    # inFields = self.getInitialFieldsDeriv(f[src,ftype,tInd+1], Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
+                    #  # - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint=True)
+                    # print inFields.shape
+                    # JTv = JTv + inFields
+                # dAsubdiag_dm_v = 0
+
+
+
+        # Missing the 0 step
+
+        # adding du_dm^T * dF_du^T * P^T vfor time 0 (no dRHS_dm_v at time 0)
+        # Asubdiag = self.getAsubdiag(0)
+        # for src in self.survey.srcList:
+        #     for projField in set(rx.projField):
+        #         v = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+        #         JTv = JTv - Utils.mkvc(self.getAdiagDeriv(0, f[src, ftype, tInd], v, adjoint = True))
+        #         # JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0] - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])), adjoint=True)
+
+        # del df_duT_v, ATinv_df_duT_v, A, Asubdiag
+        if AdiagTinv is not None:
+            AdiagTinv.clean()
+
+        return Utils.mkvc(JTv).astype(float)
+
+
+
+    def getSourceTerm(self, tInd):
+
+        Srcs = self.survey.srcList
+
+        if self._eqLocs is 'FE':
+            S_m = np.zeros((self.mesh.nF,len(Srcs)))
+            S_e = np.zeros((self.mesh.nE,len(Srcs)))
+        elif self._eqLocs is 'EF':
+            S_m = np.zeros((self.mesh.nE,len(Srcs)))
+            S_e = np.zeros((self.mesh.nF,len(Srcs)))
+
+        for i, src in enumerate(Srcs):
+            smi, sei = src.eval(self, self.times[tInd])
+            S_m[:,i] = S_m[:,i] + smi
+            S_e[:,i] = S_e[:,i] + sei
+
+        return S_m, S_e
+
+    def getInitialFields(self):
+
+        Srcs = self.survey.srcList
+
+        if self._fieldType is 'b' or self._fieldType is 'j':
+            ifields = np.zeros((self.mesh.nF, len(Srcs)))
+        elif self._fieldType is 'e' or self._fieldType is 'h':
+            ifields = np.zeros((self.mesh.nE, len(Srcs)))
+
+        for i,src in enumerate(Srcs):
+            ifields[:,i] = ifields[:,i] + getattr(src, '%sInitial'%self._fieldType, None)(self)
+
+        return ifields
+
+    def getInitialFieldsDeriv(self, src, v, adjoint=False):
+
+        if adjoint is False:
+            if self._fieldType is 'b' or self._fieldType is 'j':
+                ifieldsDeriv = np.zeros(self.mesh.nF)
+            elif self._fieldType is 'e' or self._fieldType is 'h':
+                ifieldsDeriv = np.zeros(self.mesh.nE)
+
+        elif adjoint is True:
+            ifieldsDeriv = np.zeros(self.mapping.nP)
+
+        ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv
+
+            # ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv
+            # ifieldsDeriv = self.getAdiagDeriv(None, u, v, adjoint)
+            # ifieldsDeriv = ifieldsDeriv.sum()
+
+        return ifieldsDeriv
+
+
+##########################################################################################
+################################ E-B Formulation #########################################
+##########################################################################################
+
+# ------------------------------- Problem_b -------------------------------------------- #
+
+class Problem_b(BaseTDEMProblem):
+    """
+    Starting from the quasi-static E-B formulation of Maxwell's equations (semi-discretized)
+
+    .. math::
+
+        \mathbf{C} \mathbf{e} + \\frac{\partial \mathbf{b}}{\partial t} = \mathbf{s_m} \\\\
+        \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}
+
+    where :math:`\mathbf{s_e}` is an integrated quantity, we eliminate :math:`\mathbf{e}` using
+
+    .. math::
+        \mathbf{e} = \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}
+
+    to obtain a second order semi-discretized system in :math:`\mathbf{b}`
+
+    .. math::
+        \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b}  + \\frac{\partial \mathbf{b}}{\partial t} = \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e} + \mathbf{s_m}
+
+    and moving everything except the time derivative to the rhs gives
+
+    .. math::
+        \\frac{\partial \mathbf{b}}{\partial t} = -\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} + \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e} + \mathbf{s_m}
+
+    For the time discretization, we use backward euler. To solve for the :math:`n+1`th time step, we have
+
+    .. math::
+        \\frac{\mathbf{b}^{n+1} - \mathbf{b}^{n}}{\mathbf{dt}} = -\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b}^{n+1} + \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}^{n+1} + \mathbf{s_m}^{n+1}
+
+    re-arranging to put :math:`\mathbf{b}^{n+1}` on the left hand side gives
+
+    .. math::
+        (\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}) \mathbf{b}^{n+1} = \mathbf{b}^{n} + \mathbf{dt}(\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}^{n+1} + \mathbf{s_m}^{n+1})
+
+    :param Mesh mesh: mesh
+    :param Mapping mapping: mapping
+    """
+
+    _fieldType = 'b'
+    _eqLocs    = 'FE'
+    fieldsPair = Fields_b
+    surveyPair = SurveyTDEM
+
+    def __init__(self, mesh, mapping=None, **kwargs):
+        BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs)
+
+    def getAdiag(self, tInd):
+        """
+        System matrix at a given time index
+
+        .. math::
+            (\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f})
+
+        """
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+        C = self.mesh.edgeCurl
+        MeSigmaI = self.MeSigmaI
+        MfMui = self.MfMui
+        I = Utils.speye(self.mesh.nF)
+
+        A = 1./dt * I + ( C * ( MeSigmaI * (C.T * MfMui ) ) )
+
+        if self._makeASymmetric is True:
+            return MfMui.T * A
+        return A
+
+    def getAdiagDeriv(self, tInd, u, v, adjoint=False):
+        C = self.mesh.edgeCurl
+        MeSigmaIDeriv = lambda x: self.MeSigmaIDeriv(x)
+        MfMui = self.MfMui
+
+        if adjoint:
+            if self._makeASymmetric is True:
+                v = MfMui * v
+            return  MeSigmaIDeriv(C.T * ( MfMui * u )).T * ( C.T * v )
+
+        ADeriv = ( C * ( MeSigmaIDeriv(C.T * ( MfMui * u )) * v ) )
+        if self._makeASymmetric is True:
+            return MfMui.T * ADeriv
+        return ADeriv
+
+
+    def getAsubdiag(self, tInd):
+
+        dt = self.timeSteps[tInd]
+        MfMui = self.MfMui
+        Asubdiag = - 1./dt * sp.eye(self.mesh.nF)
+
+        if self._makeASymmetric is True:
+            return MfMui.T * Asubdiag
+
+        return Asubdiag
+
+    def getAsubdiagDeriv(self, tInd, u, v, adjoint=False):
+        return Zero() * v
+
+
+
+    def getRHS(self, tInd):
+        C = self.mesh.edgeCurl
+        MeSigmaI = self.MeSigmaI
+        MfMui = self.MfMui
+
+        S_m, S_e = self.getSourceTerm(tInd)
+
+        rhs =  (C * (MeSigmaI * S_e) + S_m)
+        if self._makeASymmetric is True:
+            return MfMui.T * rhs
+        return rhs
+
+    def getRHSDeriv(self, tInd, src, v, adjoint=False):
+
+        C = self.mesh.edgeCurl
+        MeSigmaI = self.MeSigmaI
+        MeSigmaIDeriv = lambda u: self.MeSigmaIDeriv(u)
+        MfMui = self.MfMui
+
+        _, S_e = src.eval(tInd, self)
+        S_mDeriv, S_eDeriv = src.evalDeriv(self.times[tInd], self, adjoint=adjoint)
+
+        if adjoint:
+            if self._makeASymmetric is True:
+                v = self.MfMui * v
+            if isinstance(S_e, Utils.Zero):
+                MeSigmaIDerivT_v = Utils.Zero()
+            else:
+                MeSigmaIDerivT_v = MeSigmaIDeriv(S_e).T * v
+            RHSDeriv = MeSigmaIDerivT_v + S_eDeriv( MeSigmaI.T *  ( C.T * v ) ) + S_mDeriv(v)
+            return RHSDeriv
+
+        if isinstance(S_e, Utils.Zero):
+            MeSigmaIDeriv_v = Utils.Zero()
+        else:
+            MeSigmaIDeriv_v = MeSigmaIDeriv(S_e) * v
+
+        RHSDeriv = (C * (MeSigmaIDeriv_v + MeSigmaI * S_eDeriv(v) + S_mDeriv(v)))
+
+        if self._makeASymmetric is True:
+            return self.MfMui.T * RHSDeriv
+        return RHSDeriv
+
+
+# ------------------------------- Problem_e -------------------------------------------- #
+
+class Problem_e(BaseTDEMProblem):
+
+    _fieldType = 'e'
+    _eqLocs    = 'FE'
+    fieldsPair = Fields_e
+    surveyPair = SurveyTDEM
+
+    def __init__(self, mesh, mapping=None, **kwargs):
+        BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs)
+
+    def getAdiag(self, tInd):
+        """
+        System matrix at a given time index
+
+        """
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+        C = self.mesh.edgeCurl
+        MfMui = self.MfMui
+        MeSigma = self.MeSigma
+
+        return C.T * ( MfMui * C ) + 1./dt * MeSigma
+
+
+    def getAdiagDeriv(self, tInd, u, v, adjoint=False):
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+        C = self.mesh.edgeCurl
+        MfMui = self.MfMui
+        MeSigmaDeriv = self.MeSigmaDeriv(u)
+
+        if adjoint:
+            return 1./dt * MeSigmaDeriv.T * v
+
+        return 1./dt * MeSigmaDeriv * v
+
+
+    def getAsubdiag(self, tInd):
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+
+        return - 1./dt * self.MeSigma
+
+    def getAsubdiagDeriv(self, tInd, u, v, adjoint=False):
+        dt = self.timeSteps[tInd]
+
+        if adjoint:
+            return - 1./dt * self.MeSigmaDeriv(u).T * v
+
+        return - 1./dt * self.MeSigmaDeriv(u) * v
+
+    def getRHS(self, tInd):
+        return Zero()
+
+    def getRHSDeriv(self, tInd, src, v, adjoint=False):
+        return Zero()
+
+
+
+
+

SimPEG/EM/TDEM/__init__.py

@@ -1,3 +1,3 @@
-from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
-from BaseTDEM import BaseTDEMProblem, FieldsTDEM
-from TDEM_b import ProblemTDEM_b
+from TDEM import BaseTDEMProblem, Problem_b, Problem_e
+from FieldsTDEM import Fields, Fields_b
+from SurveyTDEM import Survey, Src, Rx

SimPEG/EM/TDEM_old/SurveyTDEM.py

@@ -0,0 +1,199 @@
+from SimPEG import Utils, Survey, np
+from SimPEG.Survey import BaseSurvey
+from SimPEG.EM.Utils import *
+from BaseTDEM import FieldsTDEM
+import SrcTDEM as Src
+
+class RxTDEM(Survey.BaseTimeRx):
+
+    knownRxTypes = {
+                    'ex':['e', 'Ex', 'N'],
+                    'ey':['e', 'Ey', 'N'],
+                    'ez':['e', 'Ez', 'N'],
+
+                    'bx':['b', 'Fx', 'N'],
+                    'by':['b', 'Fy', 'N'],
+                    'bz':['b', 'Fz', 'N'],
+
+                    'dbxdt':['b', 'Fx', 'CC'],
+                    'dbydt':['b', 'Fy', 'CC'],
+                    'dbzdt':['b', 'Fz', 'CC'],
+                   }
+
+    def __init__(self, locs, times, rxType):
+        Survey.BaseTimeRx.__init__(self, locs, times, rxType)
+
+    @property
+    def projField(self):
+        """Field Type projection (e.g. e b ...)"""
+        return self.knownRxTypes[self.rxType][0]
+
+    @property
+    def projGLoc(self):
+        """Grid Location projection (e.g. Ex Fy ...)"""
+        return self.knownRxTypes[self.rxType][1]
+
+    @property
+    def projTLoc(self):
+        """Time Location projection (e.g. CC N)"""
+        return self.knownRxTypes[self.rxType][2]
+
+    def getTimeP(self, timeMesh):
+        """
+            Returns the time projection matrix.
+
+            .. note::
+
+                This is not stored in memory, but is created on demand.
+        """
+        if self.rxType in ['dbxdt','dbydt','dbzdt']:
+            return timeMesh.getInterpolationMat(self.times, self.projTLoc)*timeMesh.faceDiv
+        else:
+            return timeMesh.getInterpolationMat(self.times, self.projTLoc)
+
+    def eval(self, src, mesh, timeMesh, u):
+        P = self.getP(mesh, timeMesh)
+        u_part = Utils.mkvc(u[src, self.projField, :])
+        return P*u_part
+
+    def evalDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
+        P = self.getP(mesh, timeMesh)
+
+        if not adjoint:
+            return P * Utils.mkvc(v[src, self.projField, :])
+        elif adjoint:
+            return P.T * v[src, self]
+
+
+class SrcTDEM(Survey.BaseSrc):
+    rxPair = RxTDEM
+    radius = None
+
+    def getInitialFields(self, mesh):
+        F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
+        return F0
+
+    def getJs(self, mesh, time):
+        return None
+
+
+class SrcTDEM_VMD_MVP(SrcTDEM):
+
+    def __init__(self,rxList,loc,waveformType="STEPOFF"):
+        self.loc = loc
+        self.waveformType = waveformType
+        SrcTDEM.__init__(self,rxList)
+
+    def getInitialFields(self, mesh):
+        """Vertical magnetic dipole, magnetic vector potential"""
+        if self.waveformType == "STEPOFF":
+            print ">> Step waveform: Non-zero initial condition"
+            if mesh._meshType is 'CYL':
+                if mesh.isSymmetric:
+                    MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
+                else:
+                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            elif mesh._meshType is 'TENSOR':
+                MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
+            else:
+                raise Exception('Unknown mesh for VMD')
+            return {"b": mesh.edgeCurl*MVP}
+        elif self.waveformType == "GENERAL":
+            print ">> General waveform: Zero initial condition"
+            return {"b": np.zeros(mesh.nF)}
+        else:
+            raise NotImplementedError("Only use STEPOFF or GENERAL")
+
+    def getMeS(self, mesh, MfMui):
+        if mesh._meshType is 'CYL':
+            if mesh.isSymmetric:
+                MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
+            else:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+        elif mesh._meshType is 'TENSOR':
+            MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
+        else:
+            raise Exception('Unknown mesh for VMD')
+        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
+
+
+class SrcTDEM_CircularLoop_MVP(SrcTDEM):
+    def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
+        self.loc = loc
+        self.radius = radius
+        self.waveformType = waveformType
+        SrcTDEM.__init__(self,rxList)
+
+    def getInitialFields(self, mesh):
+        """Circular Loop, magnetic vector potential"""
+        if self.waveformType == "STEPOFF":
+            print ">> Step waveform: Non-zero initial condition"
+            if mesh._meshType is 'CYL':
+                if mesh.isSymmetric:
+                    MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
+                else:
+                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            elif mesh._meshType is 'TENSOR':
+                MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
+            else:
+                raise Exception('Unknown mesh for CircularLoop')
+            return {"b": mesh.edgeCurl*MVP}
+        elif self.waveformType == "GENERAL":
+            print ">> General waveform: Zero initial condition"
+            return {"b": np.zeros(mesh.nF)}
+        else:
+            raise NotImplementedError("Only use STEPOFF or GENERAL")
+
+    def getMeS(self, mesh, MfMui):
+        if mesh._meshType is 'CYL':
+            if mesh.isSymmetric:
+                MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
+            else:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+        elif mesh._meshType is 'TENSOR':
+            MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
+        else:
+            raise Exception('Unknown mesh for CircularLoop')
+        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
+
+
+class SurveyTDEM(Survey.BaseSurvey):
+    """
+        docstring for SurveyTDEM
+    """
+    srcPair = SrcTDEM
+
+    def __init__(self, srcList, **kwargs):
+        # Sort these by frequency
+        self.srcList = srcList
+        Survey.BaseSurvey.__init__(self, **kwargs)
+
+    def projectFields(self, u):
+        data = Survey.Data(self)
+        for src in self.srcList:
+            for rx in src.rxList:
+                data[src, rx] = rx.projectFields(src, self.mesh, self.prob.timeMesh, u)
+        return data
+
+    def projectFieldsDeriv(self, u, v=None, adjoint=False):
+        assert v is not None, 'v to multiply must be provided.'
+
+        if not adjoint:
+            data = Survey.Data(self)
+            for src in self.srcList:
+                for rx in src.rxList:
+                    data[src, rx] = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v)
+            return data
+        else:
+            f = FieldsTDEM(self.mesh, self)
+            for src in self.srcList:
+                for rx in src.rxList:
+                    Ptv = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
+                    Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
+                    if rx.projField not in f: # first time we are projecting
+                        f[src, rx.projField, :] = Ptv
+                    else: # there are already fields, so let's add to them!
+                        f[src, rx.projField, :] += Ptv
+            return f
+
+

SimPEG/EM/TDEM_old/__init__.py

@@ -0,0 +1,3 @@
+from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
+from BaseTDEM import BaseTDEMProblem, FieldsTDEM
+from TDEM_b import ProblemTDEM_b

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_TDEM_1D_Inversion.py

@@ -42,10 +42,10 @@ def run(plotIt=True):
 
 
     rxOffset=1e-3
-    rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
-    src = EM.TDEM.SrcTDEM_VMD_MVP([rx], np.array([0., 0., 80]))
-    survey = EM.TDEM.SurveyTDEM([src])
-    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
+    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 80]))
+    survey = EM.TDEM.Survey([src])
+    prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
 
     prb.Solver = SolverLU
     prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)]
@@ -53,9 +53,9 @@ def run(plotIt=True):
 
     # create observed data
     std = 0.05
-    
+
     survey.dobs = survey.makeSyntheticData(mtrue,std)
-    survey.std = std 
+    survey.std = std
     survey.eps = 1e-5*np.linalg.norm(survey.dobs)
 
     if plotIt:

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.
@@ -50,57 +52,47 @@ def run(N=200, plotIt=True):
     wr = np.sum(prob.G**2.,axis=0)**0.5
     wr = ( wr/np.max(wr) )
 
-    reg = Regularization.Simple(mesh)
-    reg.wght = wr
-
+#    reg = Regularization.Simple(mesh)
+#    reg.mref = mref
+#    reg.cell_weights = wr
+#
     dmis = DataMisfit.l2_DataMisfit(survey)
     dmis.Wd = 1./wd
-
-    opt = Optimization.ProjectedGNCG(maxIter=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()
-
+#
+#    opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4)
+#    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
+#    invProb.curModel = m0
+#
+#    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
+#    target = Directives.TargetMisfit()
+#
     betaest = Directives.BetaEstimate_ByEig()
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
-
-
-    mrec = inv.run(m0)
-    ml2 = mrec
-    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
-
-    # Switch regularization to sparse
-    phim = invProb.phi_m_last
-    phid =  invProb.phi_d
+#    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
+#
+#
+#    mrec = inv.run(m0)
+#    ml2 = mrec
+#    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
+#
+#    # Switch regularization to sparse
+#    phim = invProb.phi_m_last
+#    phid =  invProb.phi_d
 
     reg = Regularization.Sparse(mesh)
+    reg.mref = mref
+    reg.cell_weights = wr
 
-#==============================================================================
-#     fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-#     dmdx = reg.mesh.cellDiffxStencil * mrec
-#     plt.plot(np.sort(dmdx))
-#==============================================================================
-
-    #reg.recModel = mrec
-    reg.wght = np.ones(mesh.nC)
     reg.mref = np.zeros(mesh.nC)
-    reg.eps_p = 5e-2
-    reg.eps_q = 1e-2
-    reg.norms   = [0., 0., 2., 2.]
-    reg.wght = wr
+    eps_p = 5e-2
+    eps_q = 5e-2
+    norms   = [0., 0., 2., 2.]
 
-    opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.)
-    beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
-    #betaest = Directives.BetaEstimate_ByEig()
-    target = Directives.TargetMisfit()
-    IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
+    opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
+    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
+    update_Jacobi = Directives.Update_lin_PreCond()
+    IRLS = Directives.Update_IRLS( norms=norms,  eps_p=eps_p, eps_q=eps_q)
 
-    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 +109,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/Examples/simpegEMpaper_example_1.py

@@ -0,0 +1,157 @@
+from SimPEG import np, Mesh, Maps, Utils, DataMisfit, Regularization, Optimization, Inversion, InvProblem, Directives
+from SimPEG import SolverLU
+from SimPEG.EM import FDEM, TDEM, mu_0
+import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.rcParams['font.size'] = 14
+
+def run(plotIt=True):
+    # Set up cylindrically symmeric mesh
+    cs, ncx, ncz, npad = 10., 15, 25, 13  # padded cyl mesh
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hz], '00C')
+
+    # Conductivity model
+    layerz = np.r_[-200., -100.]
+    layer = (mesh.vectorCCz>=layerz[0]) & (mesh.vectorCCz<=layerz[1])
+    active = mesh.vectorCCz<0.
+    sig_half  = 1e-2  # Half-space conductivity
+    sig_air   = 1e-8  # Air conductivity
+    sig_layer = 5e-2  # Layer conductivity
+    sigma = np.ones(mesh.nCz)*sig_air
+    sigma[active] = sig_half
+    sigma[layer] = sig_layer
+
+    # Mapping
+    actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
+    mtrue = np.log(sigma[active])
+
+    # FDEM problem & survey
+    rxlocs = Utils.ndgrid([np.r_[50.], np.r_[0], np.r_[0.]])
+    bzi = FDEM.Rx.Point_bSecondary(rxlocs, 'z', 'real')
+    bzr = FDEM.Rx.Point_bSecondary(rxlocs, 'z', 'imag')
+
+    freqs = np.logspace(2, 3, 5)
+    srcLoc = np.array([0., 0., 0.])
+
+    print 'min skin depth = ', 500./np.sqrt(freqs.max() * sig_half), 'max skin depth = ', 500./np.sqrt(freqs.min() * sig_half)
+    print 'max x ', mesh.vectorCCx.max(), 'min z ', mesh.vectorCCz.min(), 'max z ', mesh.vectorCCz.max()
+
+    srcList = []
+    [srcList.append(FDEM.Src.MagDipole([bzr, bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
+
+    surveyFD = FDEM.Survey(srcList)
+    prbFD = FDEM.Problem3D_b(mesh, mapping=mapping)
+    prbFD.pair(surveyFD)
+    std = 0.03
+    surveyFD.makeSyntheticData(mtrue, std)
+    surveyFD.eps = np.linalg.norm(surveyFD.dtrue)*1e-5
+
+    # FDEM inversion
+    np.random.seed(1)
+    dmisfit = DataMisfit.l2_DataMisfit(surveyFD)
+    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
+    reg = Regularization.Simple(regMesh)
+    opt = Optimization.InexactGaussNewton(maxIterCG=10, maxIter=4)
+    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
+    # Inversion Directives
+    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=3)
+    # betaest = Directives.BetaEstimate_ByEig(beta0_ratio=10.)
+    invProb.beta = 1.
+    target = Directives.TargetMisfit()
+
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta,target])
+    m0 = np.log(np.ones(mtrue.size)*sig_half)
+    reg.alpha_s = 5e-1
+    reg.alpha_x = 1.
+    prbFD.counter = opt.counter = Utils.Counter()
+    opt.LSshorten = 0.5
+    opt.tolG = 1e-10
+    opt.eps = 1e-10
+    opt.remember('xc')
+    moptFD = inv.run(m0)
+
+    # TDEM problem
+    times = np.logspace(-4, np.log10(2e-3), 10)
+    print 'min diffusion distance ', 1.28*np.sqrt(times.min()/(sig_half*mu_0)), 'max diffusion distance ', 1.28*np.sqrt(times.max()/(sig_half*mu_0))
+    rx = TDEM.Rx(rxlocs, times, 'bz')
+    src = TDEM.Src.MagDipole([rx], waveform=TDEM.Src.StepOffWaveform(), loc=srcLoc) # same src location as FDEM problem
+
+    surveyTD = TDEM.Survey([src])
+    prbTD = TDEM.Problem_b(mesh, mapping=mapping)
+    prbTD.timeSteps = [(5e-5, 10),(1e-4, 10),(5e-4, 10)]
+    prbTD.pair(surveyTD)
+    prbTD.Solver = SolverLU
+
+    std = 0.03
+    surveyTD.makeSyntheticData(mtrue, std)
+    surveyTD.std = std
+    surveyTD.eps = np.linalg.norm(surveyTD.dtrue)*1e-5
+
+    # TDEM inversion
+    dmisfit = DataMisfit.l2_DataMisfit(surveyTD)
+    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
+    reg = Regularization.Simple(regMesh)
+    opt = Optimization.InexactGaussNewton(maxIterCG=10, maxIter=4)
+    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
+
+    # Inversion Directives
+    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=3)
+    invProb.beta = 1.
+    # betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1.)
+    target = Directives.TargetMisfit()
+
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta, target])
+    m0 = np.log(np.ones(mtrue.size)*sig_half)
+    reg.alpha_s = 5e-1
+    reg.alpha_x = 1.
+    prbTD.counter = opt.counter = Utils.Counter()
+    opt.LSshorten = 0.5
+    opt.remember('xc')
+    moptTD = inv.run(m0)
+
+    if plotIt:
+        fig, ax = plt.subplots(1,1, figsize = (4, 6))
+        plt.semilogx(sigma[active], mesh.vectorCCz[active], 'k-', lw=2)
+        plt.semilogx(np.exp(moptFD), mesh.vectorCCz[active], 'ko', ms=3)
+        plt.semilogx(np.exp(moptTD), mesh.vectorCCz[active], 'k*')
+        ax.set_ylim(-1000, 0)
+        ax.set_xlim(5e-3, 1e-1)
+
+        ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
+        ax.set_ylabel('Depth (m)', fontsize = 14)
+        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
+        plt.legend(['True', 'Pred (FD)', 'Pred (TD)'], fontsize=13, loc=4)
+        plt.show()
+
+        fig = plt.figure(figsize = (10*1.3, 5*1.3))
+        ax2 = plt.subplot(122)
+        ax2.plot(times, surveyTD.dobs, 'k-', lw=2)
+        ax2.plot(times, surveyTD.dpred(moptTD), 'ko', ms=4)
+        ax2.set_xscale('log')
+        ax2.set_yscale('log')
+        ax2.set_xlim(times.min(), times.max())
+        ax1 = plt.subplot(121)
+        ax1.plot(freqs, -surveyFD.dobs[::2], 'k-', lw=2)
+        ax1.plot(freqs, -surveyFD.dobs[1::2], 'k--', lw=2)
+        dpredFD = surveyFD.dpred(moptTD)
+        ax1.plot(freqs, -dpredFD[::2], 'ko', ms=4)
+        ax1.plot(freqs, -dpredFD[1::2], 'k+', markeredgewidth=2., ms=10)
+        ax1.set_xscale('log')
+        ax1.set_yscale('log')
+        ax2.set_xlabel('Time (s)', fontsize = 14)
+        ax1.set_xlabel('Frequency (Hz)', fontsize = 14)
+        ax1.set_ylabel('Vertical magnetic field (T)', fontsize = 14)
+        ax2.grid(True,which='minor')
+        ax1.grid(True,which='minor')
+        ax2.set_title("(b) TD observed vs. predicted", fontsize = 14)
+        ax1.set_title("(a) FD observed vs. predicted", fontsize = 14)
+        ax2.legend(("Obs", "Pred"), fontsize = 12)
+        ax1.legend(("Obs", "Pred (real)", "Pred (imag)"), fontsize = 12, loc=3)
+        ax1.set_xlim(freqs.max(), freqs.min())
+        plt.show()
+
+if __name__ == '__main__':
+    run()

SimPEG/Maps.py

@@ -533,83 +533,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/MeshIO.py

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

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

SimPEG/Survey.py

@@ -1,4 +1,5 @@
 import Utils, numpy as np, scipy.sparse as sp, uuid
+import gc
 
 class BaseRx(object):
     """SimPEG Receiver Object"""

SimPEG/Utils/__init__.py

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

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/em/api_TDEM.rst

@@ -347,10 +347,10 @@ and
 
 
 
-TDEM - B formulation
-====================
+TDEM Problem
+============
 
-.. automodule:: SimPEG.EM.TDEM.TDEM_b
+.. automodule:: SimPEG.EM.TDEM.TDEM
     :show-inheritance:
     :members:
     :undoc-members:
@@ -359,7 +359,7 @@ TDEM - B formulation
 Field Storage
 =============
 
-.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.FieldsTDEM
+.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.Fields
     :show-inheritance:
     :members:
     :undoc-members:
@@ -369,19 +369,19 @@ Field Storage
 TDEM Survey Classes
 ===================
 
-.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.SurveyTDEM
+.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.Survey
     :show-inheritance:
     :members:
     :undoc-members:
     :inherited-members:
 
 
-Base Classes
-============
+.. Base Classes
+.. ============
 
-.. automodule:: SimPEG.EM.TDEM.BaseTDEM
-    :show-inheritance:
-    :members:
-    :undoc-members:
-    :inherited-members:
+.. .. automodule:: SimPEG.EM.TDEM.BaseTDEM
+..     :show-inheritance:
+..     :members:
+..     :undoc-members:
+..     :inherited-members:
 

docs/examples/DC_Forward_PseudoSection.rst

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

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:

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/em/tdem/test_TDEM_b_DerivAdjoint.py

@@ -3,310 +3,218 @@ from SimPEG import *
 from SimPEG import EM
 
 plotIt = False
-tol = 1e-6
 
-class TDEM_bDerivTests(unittest.TestCase):
+testDeriv   = True
+testAdjoint = True
 
-    def setUp(self):
+TOL = 1e-5
 
-        cs = 5.
-        ncx = 20
-        ncy = 6
-        npad = 20
-        hx = [(cs,ncx), (cs,npad,1.3)]
-        hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
-        mesh = Mesh.CylMesh([hx,1,hy], '00C')
+def setUp(prbtype='b', rxcomp='bz'):
+    cs = 5.
+    ncx = 20
+    ncy = 15
+    npad = 20
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hy], '00C')
+#
+    active = mesh.vectorCCz<0.
+    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
 
-        active = mesh.vectorCCz<0.
-        activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
-        mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
+    rxOffset = 10.
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,]
+    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 0.]))
 
-        rxOffset = 40.
-        rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
-        src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
+    survey = EM.TDEM.Survey([src])
 
-        survey = EM.TDEM.SurveyTDEM([src])
+    if prbtype == 'b':
+        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+    elif prbtype == 'e':
+        prb = EM.TDEM.Problem_e(mesh, mapping=mapping)
 
-        self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
-        # self.prb.timeSteps = [1e-5]
-        self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-        # self.prb.timeSteps = [(1e-05, 100)]
+    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+    # prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)]
 
-        try:
-            from pymatsolver import MumpsSolver
-            self.prb.Solver = MumpsSolver
-        except ImportError, e:
-            self.prb.Solver = SolverLU
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+        prb.Solver = SolverLU
 
-        self.sigma = np.ones(mesh.nCz)*1e-8
-        self.sigma[mesh.vectorCCz<0] = 1e-1
-        self.sigma = np.log(self.sigma[active])
+    m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP)
 
-        self.prb.pair(survey)
-        self.mesh = mesh
+    prb.pair(survey)
+    mesh = mesh
 
-    def test_AhVec(self):
-        """
-            Test that fields and AhVec produce consistent results
-        """
+    return prb, m, mesh
 
-        prb = self.prb
-        sigma = self.sigma
 
-        u = prb.fields(sigma)
-        Ahu = prb._AhVec(sigma, u)
+class TDEM_DerivTests(unittest.TestCase):
 
-        V1 = Ahu[:,'b',1]
-        V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0]
-        self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6)
+# ====== TEST A ========== #
 
-        V1 = Ahu[:,'e',1]
-        return np.linalg.norm(V1) < 1.e-6
+    def AderivTest(self, prbtype):
+        prb, m0, mesh = setUp(prbtype)
+        tInd = 2
+        if prbtype == 'b':
+            nu = mesh.nF
+        elif prbtype == 'e':
+            nu = mesh.nE
+        v = np.random.rand(nu)
 
-        for i in range(2,prb.nT):
+        def AderivFun(m):
+            prb.curModel = m
+            A = prb.getAdiag(tInd)
+            Av = A*v
+            prb.curModel = m0
+            ADeriv_dm = lambda dm: prb.getAdiagDeriv(tInd, v, dm)
 
-            dt = prb.timeSteps[i]
+            return Av, ADeriv_dm
 
-            V1 = Ahu[:,'b',i]
-            V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1]
-            # print np.linalg.norm(V1), np.linalg.norm(V2)
-            self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
+        print '\n Testing ADeriv %s'%(prbtype)
+        Tests.checkDerivative(AderivFun, m0, plotIt=False, num=4, eps=1e-20)
 
-            V1 = Ahu[:,'e',i]
-            V2 = prb.MeSigma*u[:,'e',i]
-            # print np.linalg.norm(V1), np.linalg.norm(V2)
-            return np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6
-
-    def test_AhVecVSMat_OneTS(self):
-
-        prb = self.prb
-        prb.timeSteps = [1e-05]
-        sigma = self.sigma
-        prb.curModel = sigma
-
-        dt = prb.timeSteps[0]
-        a11 = 1/dt*prb.MfMui*sp.identity(prb.mesh.nF)
-        a12 = prb.MfMui*prb.mesh.edgeCurl
-        a21 = prb.mesh.edgeCurl.T*prb.MfMui
-        a22 = -prb.MeSigma
-        A = sp.bmat([[a11,a12],[a21,a22]])
-
-        f = prb.fields(sigma)
-        u1 = A*f.tovec()
-        u2 = prb._AhVec(sigma,f).tovec()
-
-        self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
-
-    def test_solveAhVSMat_OneTS(self):
-        prb = self.prb
-
-        prb.timeSteps = [1e-05]
-
-        sigma = self.sigma
-        prb.curModel = sigma
-
-        dt = prb.timeSteps[0]
-        a11 = 1.0/dt*prb.MfMui*sp.identity(prb.mesh.nF)
-        a12 = prb.MfMui*prb.mesh.edgeCurl
-        a21 = prb.mesh.edgeCurl.T*prb.MfMui
-        a22 = -prb.MeSigma
-        A = sp.bmat([[a11,a12],[a21,a22]])
-
-        f = prb.fields(sigma)
-        f[:,:,0] = {'b':0}
-        f[:,'b',1] = 0
-
-        self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec()))
-
-        u1 = prb.solveAh(sigma,f).tovec().flatten()
-        u2 = sp.linalg.spsolve(A.tocsr(),f.tovec())
-
-        self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
-
-    def test_solveAhVsAhVec(self):
-
-        prb = self.prb
-        mesh = self.prb.mesh
-        sigma = self.sigma
-        self.prb.curModel = sigma
-
-        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        f[:,'b',:] = 0.0
-        for i in range(prb.nT):
-            f[:,'e', i] = np.random.rand(mesh.nE, 1)
-
-        Ahf = prb._AhVec(sigma, f)
-        f_test = prb.solveAh(sigma, Ahf)
-
-        u1 = f.tovec()
-        u2 = f_test.tovec()
-        self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
-
-    def test_DerivG(self):
-        """
-            Test the derivative of c with respect to sigma
-        """
-
-        # Random model and perturbation
-        sigma = np.random.rand(self.prb.mapping.nP)
-
-        f = self.prb.fields(sigma)
-        dm = 1000*np.random.rand(self.prb.mapping.nP)
-        h = 0.01
-
-        derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
-        print '\ntest_DerivG'
-        passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-        return passed
-
-    def test_Deriv_dUdM(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        dm = 10*np.random.rand(prb.mapping.nP)
-        f = prb.fields(sigma)
-
-        derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
-        print '\n'
-        print 'test_Deriv_dUdM'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-
-    def test_Deriv_J(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
-        d_sig = 10*np.random.rand(prb.mapping.nP)
-
-
-        derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
-        print '\n'
-        print 'test_Deriv_J'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
-
-    def test_projectAdjoint(self):
-        prb = self.prb
-        survey = prb.survey
-        mesh = self.mesh
-
-        # Generate random fields and data
-        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(prb.nT):
-            f[:,'b',i] = np.random.rand(mesh.nF, 1)
-            f[:,'e',i] = np.random.rand(mesh.nE, 1)
-        d_vec = np.random.rand(survey.nD)
-        d = Survey.Data(survey,v=d_vec)
-
-        # Check that d.T*Q*f = f.T*Q.T*d
-        V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec())
-        V2 = f.tovec().dot(survey.evalDeriv(None, v=d, adjoint=True).tovec())
-
-        self.assertTrue((V1-V2)/np.abs(V1) < tol)
-
-    def test_adjointAhVsAht(self):
-        prb = self.prb
-        mesh = self.mesh
-        sigma = self.sigma
-
-        f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            f1[:,'b',i] = np.random.rand(mesh.nF, 1)
-            f1[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-        f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            f2[:,'b',i] = np.random.rand(mesh.nF, 1)
-            f2[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-        V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec())
-        V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec())
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
-
-    # def test_solveAhtVsAhtVec(self):
-    #     prb = self.prb
-    #     mesh = self.mesh
-    #     sigma = np.random.rand(prb.mapping.nP)
-
-    #     f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey)
-    #     for i in range(1,prb.nT+1):
-    #         f1[:,'b',i] = np.random.rand(mesh.nF, 1)
-    #         f1[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-    #     f2 = prb.solveAht(sigma, f1)
-    #     f3 = prb._AhtVec(sigma, f2)
-
-    #     if True:
-    #         import matplotlib.pyplot as plt
-    #         plt.plot(f3.tovec(),'b')
-    #         plt.plot(f1.tovec(),'r')
-    #         plt.show()
-    #     V1 = np.linalg.norm(f3.tovec()-f1.tovec())
-    #     V2 = np.linalg.norm(f1.tovec())
-    #     print 'AhtVsAhtVec', V1, V2, f1.tovec()
-    #     print 'I am gunna fail this one: boo. :('
-    #     self.assertLess(V1/V2, 1e-6)
-
-    # def test_adjointsolveAhVssolveAht(self):
-    #     prb = self.prb
-    #     mesh = self.mesh
-    #     sigma = self.sigma
-
-    #     f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-    #     for i in range(1,prb.nT+1):
-    #         f1[:,'b',i] = np.random.rand(mesh.nF, 1)
-    #         f1[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-    #     f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-    #     for i in range(1,prb.nT+1):
-    #         f2[:,'b',i] = np.random.rand(mesh.nF, 1)
-    #         f2[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-    #     V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
-    #     V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
-    #     print V1, V2
-    #     self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
-
-    def test_adjointGvecVsGtvec(self):
-        mesh = self.mesh
-        prb = self.prb
+    def A_adjointTest(self,prbtype):
+        prb, m0, mesh = setUp(prbtype)
+        tInd = 2
 
+        print '\n Testing A_adjoint'
         m = np.random.rand(prb.mapping.nP)
-        sigma = np.random.rand(prb.mapping.nP)
+        if prbtype == 'b':
+            nu = prb.mesh.nF
+        elif prbtype == 'e':
+            nu = prb.mesh.nE
 
-        u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            u[:,'b',i] = np.random.rand(mesh.nF, 1)
-            u[:,'e',i] = np.random.rand(mesh.nE, 1)
+        v = np.random.rand(nu)
+        u = np.random.rand(nu)
+        prb.curModel = m0
 
-        v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            v[:,'b',i] = np.random.rand(mesh.nF, 1)
-            v[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-        V1 = m.dot(prb.Gtvec(sigma, v, u))
-        V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
-
-    def test_adjointJvecVsJtvec(self):
-        mesh = self.mesh
-        prb = self.prb
-        sigma = self.sigma
-
-        m = np.random.rand(prb.mapping.nP)
-        d = np.random.rand(prb.survey.nD)
-
-        V1 = d.dot(prb.Jvec(sigma, m))
-        V2 = m.dot(prb.Jtvec(sigma, d))
-        passed = np.abs(V1-V2)/np.abs(V1) < tol
-        print 'AdjointTest', V1, V2, passed
+        tInd = 2 # not actually used
+        V1 = v.dot(prb.getAdiagDeriv(tInd, u, m))
+        V2 = m.dot(prb.getAdiagDeriv(tInd, u, v, adjoint=True))
+        passed = np.abs(V1-V2) < TOL * (np.abs(V1) + np.abs(V2))/2.
+        print 'AdjointTest %s'%(prbtype), V1, V2, passed
         self.assertTrue(passed)
 
+    def test_Aderiv_b(self):
+        self.AderivTest('b')
+    def test_Aderiv_e(self):
+        self.AderivTest('e')
+
+    def test_Aadjoint_b(self):
+        self.A_adjointTest('b')
+    def test_Aadjoint_e(self):
+        self.A_adjointTest('e')
+
+# ====== TEST Fields Deriv Pieces ========== #
+
+    def test_eDeriv_m_adjoint(self):
+        prb, m0, mesh = setUp()
+        tInd = 0
+
+        v = np.random.rand(mesh.nF)
+
+        print '\n Testing eDeriv_m Adjoint'
+
+        prb, m0, mesh = setUp()
+        f = prb.fields(m0)
+
+        m = np.random.rand(prb.mapping.nP)
+        e = np.random.randn(prb.mesh.nE)
+        V1 = e.dot(f._eDeriv_m(1, prb.survey.srcList[0], m))
+        V2 = m.dot(f._eDeriv_m(1, prb.survey.srcList[0], e, adjoint=True))
+        tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+        passed = np.abs(V1-V2) < tol
+
+        print '    ', V1, V2, np.abs(V1-V2), tol, passed
+        self.assertTrue(passed)
+
+    def test_eDeriv_u_adjoint(self):
+        print '\n Testing eDeriv_u Adjoint'
+
+        prb, m0, mesh = setUp()
+        f = prb.fields(m0)
+
+        b = np.random.rand(prb.mesh.nF)
+        e = np.random.randn(prb.mesh.nE)
+        V1 = e.dot(f._eDeriv_u(1, prb.survey.srcList[0], b))
+        V2 = b.dot(f._eDeriv_u(1, prb.survey.srcList[0], e, adjoint=True))
+        tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+        passed = np.abs(V1-V2) < tol
+
+        print '    ', V1, V2, np.abs(V1-V2), tol, passed
+        self.assertTrue(passed)
+
+
+# ====== TEST Jvec ========== #
+
+    if testDeriv:
+
+        def JvecTest(self, prbtype, rxcomp):
+            prb, m, mesh = setUp(prbtype, rxcomp)
+
+            derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)]
+            print '\n'
+            print 'test_Jvec_%s_%s' %(prbtype, rxcomp)
+            Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20)
+
+        def test_Jvec_b_bx(self):
+            self.JvecTest('b','bx')
+
+        def test_Jvec_b_bz(self):
+            self.JvecTest('b','bz')
+
+        def test_Jvec_b_dbxdt(self):
+            self.JvecTest('b','dbxdt')
+
+        def test_Jvec_b_dbzdt(self):
+            self.JvecTest('b','dbzdt')
+
+        def test_Jvec_b_ey(self):
+            self.JvecTest('b','ey')
+
+        def test_Jvec_e_ey(self):
+            self.JvecTest('e','ey')
+
+
+# ====== TEST Jtvec ========== #
+
+    if testAdjoint:
+
+        def JvecVsJtvecTest(self, prbtype='b', rxcomp='bz'):
+
+            print '\nAdjoint Testing Jvec, Jtvec %s' %(rxcomp)
+
+            prb, m0, mesh = setUp(prbtype, rxcomp)
+            m = np.random.rand(prb.mapping.nP)
+            d = np.random.randn(prb.survey.nD)
+            V1 = d.dot(prb.Jvec(m0, m))
+            V2 = m.dot(prb.Jtvec(m0, d))
+            tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+            passed = np.abs(V1-V2) < tol
+
+            print '    ', V1, V2, np.abs(V1-V2), tol, passed
+            self.assertTrue(passed)
+
+        def test_Jvec_adjoint_b_bx(self):
+            self.JvecVsJtvecTest('b', 'bx')
+
+        def test_Jvec_adjoint_b_bz(self):
+            self.JvecVsJtvecTest('b', 'bz')
+
+        def test_Jvec_adjoint_b_dbxdt(self):
+            self.JvecVsJtvecTest('b', 'bx')
+
+        def test_Jvec_adjoint_b_dbzdt(self):
+            self.JvecVsJtvecTest('b', 'bz')
+
+        def test_Jvec_adjoint_b_ey(self):
+            self.JvecVsJtvecTest('b', 'ey')
+
+        # This is not working because Problem_e has not done
+        # def test_Jvec_adjoint_e_ey(self):
+        #     self.JvecVsJtvecTest('e', 'ey')
 
 
 

tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py

@@ -3,10 +3,12 @@ from SimPEG import *
 from SimPEG import EM
 
 plotIt = False
+testDeriv = True
+testAdjoint = True
 
-class TDEM_bDerivTests(unittest.TestCase):
+TOL = 1e-5
 
-    def setUp(self):
+def setUp(self, rxcomp='bz'):
 
         cs = 5.
         ncx = 20
@@ -21,131 +23,78 @@ class TDEM_bDerivTests(unittest.TestCase):
         mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
 
         rxOffset = 40.
-        rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
-        src = EM.TDEM.SrcTDEM_VMD_MVP( [rx], loc=np.array([0., 0., 0.]))
-        rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz')
-        src2 = EM.TDEM.SrcTDEM_VMD_MVP( [rx2], loc=np.array([0., 0., 0.]))
+        rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), rxcomp)
+        src = EM.TDEM.Src.MagDipole( [rx], loc=np.array([0., 0., 0.]))
+        rx2 = EM.TDEM.Rx(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), rxcomp)
+        src2 = EM.TDEM.Src.MagDipole( [rx2], loc=np.array([0., 0., 0.]))
 
-        survey = EM.TDEM.SurveyTDEM([src,src2])
+        survey = EM.TDEM.Survey([src,src2])
 
-        self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
-        # self.prb.timeSteps = [1e-5]
-        self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-        # self.prb.timeSteps = [(1e-05, 100)]
+        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+        # prb.timeSteps = [1e-5]
+        prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+        # prb.timeSteps = [(1e-05, 100)]
 
         try:
             from pymatsolver import MumpsSolver
-            self.prb.Solver = MumpsSolver
+            prb.Solver = MumpsSolver
         except ImportError, e:
-            self.prb.Solver  = SolverLU
+            prb.Solver  = SolverLU
 
-        self.sigma = np.ones(mesh.nCz)*1e-8
-        self.sigma[mesh.vectorCCz<0] = 1e-1
-        self.sigma = np.log(self.sigma[active])
+        m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP)
 
-        self.prb.pair(survey)
-        self.mesh = mesh
+        prb.pair(survey)
 
-    def test_DerivG(self):
-        """
-            Test the derivative of c with respect to sigma
-        """
+        return mesh, prb, m
 
-        # Random model and perturbation
-        sigma = np.random.rand(self.prb.mapping.nP)
-
-        f = self.prb.fields(sigma)
-        dm = 1000*np.random.rand(self.prb.mapping.nP)
-        h = 0.01
-
-        derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
-        print '\ntest_DerivG'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-
-    def test_Deriv_dUdM(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        dm = 10*np.random.rand(prb.mapping.nP)
-        f = prb.fields(sigma)
-
-        derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
-        print '\n'
-        print 'test_Deriv_dUdM'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-
-    def test_Deriv_J(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
-        d_sig = 10*np.random.rand(prb.mapping.nP)
+class TDEM_bDerivTests(unittest.TestCase):
 
 
-        derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
-        print '\n'
-        print 'test_Deriv_J'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
+    if testDeriv:
+        def Deriv_J(self, rxcomp='bz'):
 
-    def test_projectAdjoint(self):
-        prb = self.prb
-        survey = prb.survey
-        nSrc = survey.nSrc
-        mesh = self.mesh
+            mesh, prb, m0 = setUp(rxcomp)
 
-        # Generate random fields and data
-        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(prb.nT):
-            f[:,'b',i] = np.random.rand(mesh.nF, nSrc)
-            f[:,'e',i] = np.random.rand(mesh.nE, nSrc)
-        d_vec = np.random.rand(survey.nD)
-        d = Survey.Data(survey,v=d_vec)
+            prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
 
-        # Check that d.T*Q*f = f.T*Q.T*d
-        V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec())
-        V2 = np.sum((f.tovec())*(survey.evalDeriv(None, v=d, adjoint=True).tovec()))
+            derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m0, mx)]
+            print '\n'
+            print 'test_Deriv_J %s'%rxcomp
+            Tests.checkDerivative(derChk, m0, plotIt=False, num=3, eps=1e-20)
 
-        self.assertTrue((V1-V2)/np.abs(V1) < 1e-6)
+        def test_Jvec_bx(self):
+            self.Deriv_J('bx')
 
-    def test_adjointGvecVsGtvec(self):
-        mesh = self.mesh
-        prb = self.prb
+        def test_Jvec_bz(self):
+            self.Deriv_J('bz')
 
-        m = np.random.rand(prb.mapping.nP)
-        sigma = np.random.rand(prb.mapping.nP)
+        def test_Jvec_ey(self):
+            self.Deriv_J('ey')
 
-        u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            u[:,'b',i] = np.random.rand(mesh.nF, 2)
-            u[:,'e',i] = np.random.rand(mesh.nE, 2)
+    if testAdjoint:
+        def adjointJvecVsJtvec(self, rxcomp='bz'):
+            print ' \n Testing Adjoint %s' %rxcomp
+            mesh, prb, m0 = setUp(rxcomp)
 
-        v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            v[:,'b',i] = np.random.rand(mesh.nF, 2)
-            v[:,'e',i] = np.random.rand(mesh.nE, 2)
+            m = np.random.rand(prb.mapping.nP)
+            d = np.random.rand(prb.survey.nD)
 
-        V1 = m.dot(prb.Gtvec(sigma, v, u))
-        V2 = np.sum(v.tovec()*prb.Gvec(sigma, m, u).tovec())
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) <1e-6)
+            V1 = d.dot(prb.Jvec(m0, m))
+            V2 = m.dot(prb.Jtvec(m0, d))
 
-    def test_adjointJvecVsJtvec(self):
-        mesh = self.mesh
-        prb = self.prb
-        sigma = self.sigma
+            tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+            passed = np.abs(V1-V2) < tol
+            print '    ', V1, V2, np.abs(V1-V2), tol, passed
+            self.assertTrue(passed)
 
-        m = np.random.rand(prb.mapping.nP)
-        d = np.random.rand(prb.survey.nD)
+        def test_JvecVsJtvec_bx(self):
+            self.adjointJvecVsJtvec('bx')
 
-        V1 = d.dot(prb.Jvec(sigma, m))
-        V2 = m.dot(prb.Jtvec(sigma, d))
-        print 'AdjointTest', V1, V2
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < 1e-6)
+        def test_JvecVsJtvec_bz(self):
+            self.adjointJvecVsJtvec('bz')
+
+        def test_JvecVsJtvec_ey(self):
+            self.adjointJvecVsJtvec('ey')
 
 
 

tests/em/tdem/test_TDEM_combos.py

@@ -1,94 +0,0 @@
-import unittest
-from SimPEG import *
-from SimPEG import EM
-
-plotIt = False
-
-def getProb(meshType='CYL',rxTypes='bx,bz',nSrc=1):
-    cs = 5.
-    ncx = 20
-    ncy = 6
-    npad = 20
-    hx = [(cs,ncx), (cs,npad,1.3)]
-    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
-    mesh = Mesh.CylMesh([hx,1,hy], '00C')
-
-    active = mesh.vectorCCz<0.
-    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
-    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
-
-    rxOffset = 40.
-
-    srcs = []
-    for ii in range(nSrc):
-        rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')]
-        srcs += [EM.TDEM.SrcTDEM_VMD_MVP(rxs,np.array([0., 0., 0.]))]
-
-    survey = EM.TDEM.SurveyTDEM(srcs)
-
-    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
-    # prb.timeSteps = [1e-5]
-    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-    # prb.timeSteps = [(1e-05, 100)]
-
-    try:
-        from pymatsolver import MumpsSolver
-        prb.Solver = MumpsSolver
-    except ImportError, e:
-        prb.Solver  = SolverLU
-
-    sigma = np.ones(mesh.nCz)*1e-8
-    sigma[mesh.vectorCCz<0] = 1e-1
-    sigma = np.log(sigma[active])
-
-    prb.pair(survey)
-    return prb, mesh, sigma
-
-def dotestJvec(prb, mesh, sigma):
-    prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-    # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
-    d_sig = 10*np.random.rand(prb.mapping.nP)
-    derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
-    return Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20)
-
-def dotestAdjoint(prb, mesh, sigma):
-    m = np.random.rand(prb.mapping.nP)
-    d = np.random.rand(prb.survey.nD)
-
-    V1 = d.dot(prb.Jvec(sigma, m))
-    V2 = m.dot(prb.Jtvec(sigma, d))
-    print 'AdjointTest', V1, V2
-    return np.abs(V1-V2)/np.abs(V1), 1e-6
-
-class TDEM_bDerivTests(unittest.TestCase):
-
-    def test_Jvec_bx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx')))
-    def test_Adjoint_bx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx')))
-
-    def test_Jvec_bxbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz')))
-    def test_Adjoint_bxbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz')))
-
-    def test_Jvec_bxbz_2src(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nSrc=2)))
-    def test_Adjoint_bxbz_2src(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nSrc=2)))
-
-    def test_Jvec_bxbzbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz,bz')))
-    def test_Adjoint_bxbzbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz,bz')))
-
-    def test_Jvec_dbxdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbxdt')))
-    def test_Adjoint_dbxdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbxdt')))
-
-    def test_Jvec_dbzdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbzdt')))
-    def test_Adjoint_dbzdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbzdt')))
-
-    def test_Jvec_dbxdtbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbxdt,bz')))
-    def test_Adjoint_dbxdtbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbxdt,bz')))
-
-    def test_Jvec_ey(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey')))
-    def test_Adjoint_ey(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey')))
-
-    def test_Jvec_eybzdbxdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey,bz,dbxdt')))
-    def test_Adjoint_eybzdbxdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey,bz,dbxdt')))
-
-
-if __name__ == '__main__':
-    unittest.main()

tests/em/tdem/test_TDEM_crosscheck.py

@@ -0,0 +1,76 @@
+import unittest
+from SimPEG import *
+from SimPEG import EM
+
+TOL = 1e-5
+FLR = 1e-20
+
+np.random.seed(seed=25) # set a seed so that the same conductivity model is used for all runs
+
+def setUp(prbtype = 'b', rxcomp='bz'):
+    cs = 5.
+    ncx = 20
+    ncy = 15
+    npad = 20
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hy], '00C')
+#
+    active = mesh.vectorCCz<0.
+    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
+
+    rxOffset = 10.
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,]
+    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 0.]))
+
+    survey = EM.TDEM.Survey([src])
+
+    if prbtype == 'b':
+        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+    elif prbtype == 'e':
+        prb = EM.TDEM.Problem_e(mesh, mapping=mapping)
+
+    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+    # prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)]
+
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+        prb.Solver = SolverLU
+
+    m = np.log(1e-1)*np.ones(prb.mapping.nP) #+ 1e-2*np.random.randn(prb.mapping.nP)
+
+    prb.pair(survey)
+    mesh = mesh
+
+    return prb, m, mesh
+
+def CrossCheck(prbtype1='b', prbtype2='e', rxcomp='bz'):
+
+    prb1,m1,mesh1 = setUp(prbtype1, rxcomp)
+    prb2,m2,mesh2 = setUp(prbtype2, rxcomp)
+
+    assert (m1 == m2).all(), 'Models for two formulations are different'
+
+    d1 = prb1.survey.dpred(m1)
+    d2 = prb2.survey.dpred(m2)
+
+
+    check = np.linalg.norm(d1 - d2)
+    tol   = 0.5 * (np.linalg.norm(d1) + np.linalg.norm(d2)) * TOL
+    passed = check < tol
+
+    print 'Checking %s, %s for %s data'%(prbtype1, prbtype2, rxcomp)
+    print '  ', np.linalg.norm(d1), np.linalg.norm(d2), np.linalg.norm(check), tol, passed
+
+    assert passed
+
+class TDEM_cross_check_EB(unittest.TestCase):
+    def test_EB_ey(self):
+        CrossCheck('b','e','ey')
+
+if __name__ == '__main__':
+    unittest.main()
+

tests/em/tdem/test_TDEM_forward_Analytic.py

@@ -29,12 +29,12 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
     actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
     mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
 
-    rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
-    src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
+    src = EM.TDEM.Src.MagDipole([rx], waveform= EM.TDEM.Src.StepOffWaveform(), loc=np.array([0., 0., 0.]))
     # src = EM.TDEM.SrcTDEM([rx], loc=np.array([0., 0., 0.]))
 
-    survey = EM.TDEM.SurveyTDEM([src])
-    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
+    survey = EM.TDEM.Survey([src])
+    prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
     prb.Solver = MumpsSolver
 
     prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40), (0.0001, 40), (0.0005, 40)]
@@ -50,6 +50,8 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
 
     ind = np.logical_and(rx.times > bounds[0],rx.times < bounds[1])
     log10diff = np.linalg.norm(np.log10(np.abs(bz_calc[ind])) - np.log10(np.abs(bz_ana[ind])))/np.linalg.norm(np.log10(np.abs(bz_ana[ind])))
+
+    print ' |bz_ana| = ',np.linalg.norm(bz_ana), ' |bz_num| = ', np.linalg.norm(bz_calc), ' |bz_ana - bz_num| =', np.linalg.norm(bz_ana-bz_calc)
     print 'Difference: ', log10diff
 
     if showIt == True:
@@ -61,6 +63,12 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
     return log10diff
 
 
+class TDEM_SimpleSrcTests(unittest.TestCase):
+    def test_source(self):
+        waveform = EM.TDEM.Src.StepOffWaveform()
+        assert waveform.eval(0.) == 0.
+
+
 class TDEM_bTests(unittest.TestCase):
 
     def test_analytic_p2_CYL_50m(self):