Minor changes...test runs fine from here

Fixed the boolean error in Directives.

SimPEG/DCIP/DCIPUtils.py

@@ -169,7 +169,7 @@ def readUBC_DC2DModel(fileName):
 
     return model
 
-def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
+def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None):
     """
         Read list of 2D tx-rx location and plot a speudo-section of apparent
         resistivity.
@@ -192,9 +192,6 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
     from scipy.interpolate import griddata
     import pylab as plt
 
-    # Set depth to 0 for now
-    z0 = 0.
-
     # Pre-allocate
     midx = []
     midz = []
@@ -259,38 +256,53 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
         midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
         midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ])
 
-    ax = axs
-
     # Grid points
     grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
     grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear')
-
+    
+    # Scale the color scheme
     if clim == None:
         vmin, vmax = rho.min(), rho.max()
     else:
         vmin, vmax = clim[0], clim[1]
 
+    # Plot data
     grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho)
-    ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, clim=(vmin, vmax))
-    cbar = plt.colorbar(format="$10^{%.1f}$",fraction=0.04,orientation="horizontal")
-
-    cmin,cmax = cbar.get_clim()
-    ticks = np.linspace(cmin,cmax,3)
-    cbar.set_ticks(ticks)
-    cbar.ax.tick_params(labelsize=10)
+    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 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)
+    if contour is not None:
+        plt.contour(grid_x,grid_z,grid_rho,levels = contour,colors = 'r', vmin = vmin, vmax = vmax)
+        
+    # Add scatter points
+    axs.scatter(midx,midz,s=10,c=rho.T, vmin = vmin, vmax = vmax)
+    
+    if colorbar:
+        
+        if dtype == 'volt':
+            cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal")
 
-    # Plot apparent resistivity
-    ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
+        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)
 
-    #ax.set_xticklabels([])
-    #ax.set_yticklabels([])
+
+
+    if not axlabel:
+        axs.set_xticklabels([])
+        axs.set_yticklabels([])
 
     plt.gca().set_aspect('equal', adjustable='box')
 
@@ -448,15 +460,15 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
     survey = DC.SurveyDC(SrcList)
     return survey, Tx, Rx
 
-def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
+def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
     """
         Write UBC GIF DCIP 2D or 3D observation file
 
         Input:
-        :string fileName -> including path where the file is written out
-        :DCsurvey -> DC survey class object
-        :string dtype ->  either '2D' | '3D'
-        :string  stype ->  either 'SURFACE' | 'GENERAL'
+        :string     fileName -> including path where the file is written out
+        :DCsurvey   DC survey class object
+        :string     dtype ->  either '2D' | '3D'
+        :string     stype ->  either 'SURFACE' | 'GENERAL'
 
         Output:
         :param UBC2D-Data file
@@ -471,10 +483,16 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
 
     assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'"
     assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
-
+    
     fid = open(fileName,'w')
-    fid.write('! ' + stype + ' FORMAT\n')
-
+    
+    
+    if iptype!=0:
+        fid.write('IPTYPE=%i\n'%iptype)
+        
+    else:
+        fid.write('! ' + stype + ' FORMAT\n')
+        
     count = 0
 
     for ii in range(DCsurvey.nSrc):
@@ -498,7 +516,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
                 B = np.repeat(tx[0,1],M.shape[0],axis=0)
                 M = M[:,0]
                 N = N[:,0]
-
+                
                 np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
 
 
@@ -506,18 +524,25 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
 
                 if stype == 'SURFACE':
 
-                    fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
+                    fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
                     M = M[:,0]
                     N = N[:,0]
 
                 if stype == 'GENERAL':
 
+                    # Flip sign for z-elevation to depth
+                    tx[2::2,:] = -tx[2::2,:]
+                    
                     fid.writelines("%e " % ii for ii in mkvc(tx[::2,:]))
                     M = M[:,0::2]
                     N = N[:,0::2]
 
+                    # Flip sign for z-elevation to depth
+                    M[:,1::2] = -M[:,1::2]
+                    N[:,1::2] = -N[:,1::2]
+                    
                 fid.write('%i\n'% nD)
-                np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
+                np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n')
 
         if dtype=='3D':
 
@@ -529,11 +554,12 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
 
             if stype == 'GENERAL':
 
-                fid.writelines("%e " % ii for ii in mkvc(tx))
+                fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
 
             fid.write('%i\n'% nD)
             np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
-
+            fid.write('\n')
+            
         count += nD
 
     fid.close()
@@ -640,51 +666,59 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'):
     DCsurvey2D.std = np.asarray(DCsurvey.std)
 
     return DCsurvey2D
-
-def readUBC_DC3Dobs(fileName):
+    
+def readUBC_DC3Dobs(fileName, dtype = 'DC'):
     """
-        Read UBC GIF DCIP 3D observation file and generate survey
+        Read UBC GIF IP 3D observation file and generate survey
 
         Input:
         :param fileName, path to the UBC GIF 3D obs file
 
         Output:
-        :param DCIPsurvey
+        :param IPsurvey
         :return
 
-        Created on Mon April 6th, 2015
-
         @author: dominiquef
 
     """
-
+    zflag = True # Flag for z value provided
+    
     # Load file
-    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
-
+    if dtype == 'IP':
+        obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='IPTYPE')
+        
+    elif dtype == 'DC':
+        obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
+        
+    else:
+        print "dtype must be 'DC'(default) | 'IP'"
+        
     # Pre-allocate
     srcLists = []
     Rx = []
     d = []
     wd = []
-    zflag = True # Flag for z value provided
+    
 
     # Countdown for number of obs/tx
     count = 0
     for ii in range(obsfile.shape[0]):
 
+        # Skip if blank line
         if not obsfile[ii]:
             continue
 
-        # First line is transmitter with number of receivers
+        # First line or end of a transmitter block, read transmitter info
         if count==0:
-
-            temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
+            # Read the line
+            temp = (np.fromstring(obsfile[ii], dtype=float, sep=' ').T)
             count = int(temp[-1])
 
             # Check if z value is provided, if False -> nan
             if len(temp)==5:
-                tx = np.r_[temp[0:2],np.nan,temp[0:2],np.nan]
-                zflag = False
+                tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan]
+                
+                zflag = False # Pass on the flag to the receiver loc
 
             else:
                 tx = temp[:-1]
@@ -692,8 +726,16 @@ def readUBC_DC3Dobs(fileName):
             rx = []
             continue
 
-        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
+        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string
 
+        # Filter out negative IP
+#        if temp[-2] < 0:       
+#            count = count -1
+#            print "Negative!"
+#            
+#        else:
+            
+        # If the Z-location is provided, otherwise put nan
         if zflag:
 
             rx.append(temp[:-2])
@@ -703,7 +745,7 @@ def readUBC_DC3Dobs(fileName):
                 wd.append(temp[-1])
 
         else:
-            rx.append(np.r_[temp[0:2],np.nan,temp[0:2],np.nan] )
+            rx.append(np.r_[temp[0:2],np.nan,temp[2:4],np.nan] )
             # Check if there is data with the location
             if len(temp)==6:
                 d.append(temp[-2])
@@ -711,7 +753,7 @@ def readUBC_DC3Dobs(fileName):
 
         count = count -1
 
-        # Reach the end of transmitter block
+        # Reach the end of transmitter block, append the src, rx and continue
         if count == 0:
             rx = np.asarray(rx)
             Rx = DC.RxDipole(rx[:,:3],rx[:,3:])

SimPEG/Directives.py

@@ -222,7 +222,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
             mref = 0
         mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
         phi_mx = 0.5 * mx.dot(mx)
-        if self.prob.mesh.dim==2:
+        if self.prob.mesh.dim >= 2:
             my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
             phi_my = 0.5 * my.dot(my)
         else:
@@ -237,41 +237,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
         # Save the file as a npz
         np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
 
-#==============================================================================
-# class SaveOutputDictEveryIteration(_SaveEveryIteration):
-#     """SaveOutputDictEveryIteration
-#     A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
-#     """
-#
-#     def initialize(self):
-#         print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName
-#
-#     def endIter(self):
-#         # Save the data.
-#         ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
-#         phi_ms = 0.5*ms.dot(ms)
-#         if self.reg.mrefInSmooth == True:
-#             mref = self.reg.mref
-#         else:
-#             mref = 0
-#         mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
-#         phi_mx = 0.5 * mx.dot(mx)
-#         if self.prob.mesh.dim==2:
-#             my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
-#             phi_my = 0.5 * my.dot(my)
-#         else:
-#             phi_my = 'NaN'
-#         if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType:
-#             mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
-#             phi_mz = 0.5 * mz.dot(mz)
-#         else:
-#             phi_mz = 'NaN'
-#
-#
-#         # Save the file as a npz
-#         np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
-#
-#==============================================================================
 
 # class UpdateReferenceModel(Parameter):
 
@@ -293,6 +258,7 @@ class Update_IRLS(InversionDirective):
     phi_m_last = None
     phi_d_last = None
 
+
     def initialize(self):
 
         # Scale the regularization for changes in norm
@@ -305,12 +271,13 @@ class Update_IRLS(InversionDirective):
 
         self.reg.curModel = self.invProb.curModel
         self.reg.gamma = self.gamma
-
+        
+        self.reg._W = None
         if getattr(self, 'phi_d_last', None) is None:
             self.phi_d_last = self.invProb.phi_d
 
     def endIter(self):
-        # Cool the threshold parameter
+        # Cool the threshold parameter if required
         if getattr(self, 'factor', None) is not None:
             eps = self.reg.eps / self.factor
 
@@ -321,32 +288,49 @@ class Update_IRLS(InversionDirective):
 
         # Get phi_m at the end of current iteration
         self.phi_m_last = self.invProb.phi_m_last
-
+        
+        self.reg._W = None
          # Update the model used for the IRLS weights
         self.reg.curModel = self.invProb.curModel
 
-        # Temporarely set gamma to 1.
+        # Temporarely set gamma to 1. to get raw phi_m
         self.reg.gamma = 1.
 
-        # Compute change in model objective function and update scaling
+        # 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
-
-        self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
+        
+        # Set the weighting matrix to None so that it is recomputed next time
+        # it is called in the inversion
+        self.reg._W = None
 
 class Update_lin_PreCond(InversionDirective):
-
-
+    """
+    Create a Jacobi preconditioner for the linear problem
+    """
+    onlyOnStart=False
+    
+    def initialize(self):
+    
+        if getattr(self.opt, 'approxHinv', None) is None:
+            # Update the pre-conditioner
+            diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
+            self.opt.approxHinv = PC
+            
     def endIter(self):
         # Cool the threshold parameter
-
+        if self.onlyOnStart==True:
+            return
+            
         if getattr(self.opt, 'approxHinv', None) is not None:
             # Update the pre-conditioner
-            diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
-            PC     = Utils.sdiag(diagA**-1.)
+            diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
             self.opt.approxHinv = PC
-            print 'Updated pre-cond'
+
 
 class Update_Wj(InversionDirective):
     """
@@ -373,3 +357,19 @@ 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/Examples/Inversion_IRLS.py

@@ -86,12 +86,12 @@ def run(N=200, plotIt=True):
     #reg.recModel = mrec
     reg.wght = np.ones(mesh.nC)
     reg.mref = np.zeros(mesh.nC)
-    reg.eps_p = 2e-3
-    reg.eps_q = 2e-3
+    reg.eps_p = 5e-2
+    reg.eps_q = 1e-2
     reg.norms   = [0., 0., 2., 2.]
     reg.wght = wr
 
-    opt = Optimization.ProjectedGNCG(maxIter=5 ,lower=-2.,upper=2., maxIterCG= 100, tolCG = 1e-3)
+    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()

SimPEG/MT/Utils/ediFilesUtils.py

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

SimPEG/Mesh/MeshIO.py

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

SimPEG/Regularization.py

@@ -663,9 +663,8 @@ class Sparse(Simple):
             self.Rs = Utils.speye(self.regmesh.nC)
 
         else:
-            f_m = self.curModel - self.reg.mref
+            f_m = self.mapping * (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
@@ -679,7 +678,7 @@ class Sparse(Simple):
             self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
 
         else:
-            f_m = self.regmesh.cellDiffxStencil * self.curModel
+            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 )
 
@@ -693,7 +692,7 @@ class Sparse(Simple):
             self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
 
         else:
-            f_m = self.regmesh.cellDiffyStencil * self.curModel
+            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 )
 
@@ -707,7 +706,7 @@ class Sparse(Simple):
             self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
 
         else:
-            f_m = self.regmesh.cellDiffzStencil * self.curModel
+            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 )
 
@@ -728,14 +727,15 @@ class Sparse(Simple):
     @property
     def W(self):
         """Full regularization matrix W"""
-        #if getattr(self, '_W', None) is None:
-        wlist = (self.Wsmall, self.Wsmooth)
-        #self._W = sp.vstack(wlist)
-        return sp.vstack(wlist)
+        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 = (eps**(1-exponent/2.))**0.5
-        r = eta / (f_m**2.+ eps**2.)**((1-exponent/2.)/2.)
+        eta = (eps**(1.-exponent/2.))**0.5
+        r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)
 
         return r

SimPEG/Utils/io_utils.py

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

tests/mt/test_Problem3D_againstAnalytic.py

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