no message

README.rst

@@ -25,10 +25,6 @@ SimPEG
     :target: https://coveralls.io/r/simpeg/simpeg?branch=master
     :alt: Coverage status
 
-.. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square
-    :alt: gitter chat room at https://gitter.im/simpeg/simpeg
-    :target: https://gitter.im/simpeg/simpeg
-    
 Simulation and Parameter Estimation in Geophysics  -  A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
 
 The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:

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, cblabel=True, axlabel = True, colorbar = True, contour = None):
+def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None):
     """
         Read list of 2D tx-rx location and plot a speudo-section of apparent
         resistivity.
@@ -192,6 +192,9 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl
     from scipy.interpolate import griddata
     import pylab as plt
 
+    # Set depth to 0 for now
+    z0 = 0.
+
     # Pre-allocate
     midx = []
     midz = []
@@ -254,55 +257,40 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cbl
                 break
 
         midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
-        midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ])
+        midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + z0 ])
+
+    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, vmin = vmin, vmax = vmax)
-    plt.gca().tick_params(axis='both', which='major', labelsize=8)
+    ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, clim=(vmin, vmax))
+    cbar = plt.colorbar(format="$10^{%.1f}$",fraction=0.04,orientation="horizontal")
+
+    cmin,cmax = cbar.get_clim()
+    ticks = np.linspace(cmin,cmax,3)
+    cbar.set_ticks(ticks)
+    cbar.ax.tick_params(labelsize=10)
     
-    if 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")
+    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)
 
-        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)
+    # Plot apparent resistivity
+    ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax))
 
-
-
-    if not axlabel:
-        axs.set_xticklabels([])
-        axs.set_yticklabels([])
+    #ax.set_xticklabels([])
+    #ax.set_yticklabels([])
 
     plt.gca().set_aspect('equal', adjustable='box')
 
@@ -460,15 +448,15 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
     survey = DC.SurveyDC(SrcList)
     return survey, Tx, Rx
 
-def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
+def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
     """
         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
@@ -483,16 +471,10 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
 
     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')
-        
+    fid.write('! ' + stype + ' FORMAT\n')
+
     count = 0
 
     for ii in range(DCsurvey.nSrc):
@@ -516,7 +498,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 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')
 
 
@@ -524,25 +506,18 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
 
                 if stype == 'SURFACE':
 
-                    fid.writelines("%f " % ii for ii in mkvc(tx[0,:]))
+                    fid.writelines("%e " % 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')
+                np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
 
         if dtype=='3D':
 
@@ -554,12 +529,11 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0):
 
             if stype == 'GENERAL':
 
-                fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:]))
+                fid.writelines("%e " % ii for ii in mkvc(tx))
 
             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()
@@ -666,59 +640,51 @@ 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):
     """
-        Read UBC GIF IP 3D observation file and generate survey
+        Read UBC GIF DCIP 3D observation file and generate survey
 
         Input:
         :param fileName, path to the UBC GIF 3D obs file
 
         Output:
-        :param IPsurvey
+        :param DCIPsurvey
         :return
 
+        Created on Mon April 6th, 2015
+
         @author: dominiquef
 
     """
-    zflag = True # Flag for z value provided
-    
+
     # Load file
-    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'"
-        
+    obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
+
     # 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 or end of a transmitter block, read transmitter info
+        # First line is transmitter with number of receivers
         if count==0:
-            # Read the line
-            temp = (np.fromstring(obsfile[ii], dtype=float, sep=' ').T)
+
+            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[2:4],np.nan]
-                
-                zflag = False # Pass on the flag to the receiver loc
+                tx = np.r_[temp[0:2],np.nan,temp[0:2],np.nan]
+                zflag = False
 
             else:
                 tx = temp[:-1]
@@ -726,16 +692,8 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'):
             rx = []
             continue
 
-        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string
+        temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
 
-        # 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])
@@ -745,7 +703,7 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'):
                 wd.append(temp[-1])
 
         else:
-            rx.append(np.r_[temp[0:2],np.nan,temp[2:4],np.nan] )
+            rx.append(np.r_[temp[0:2],np.nan,temp[0:2],np.nan] )
             # Check if there is data with the location
             if len(temp)==6:
                 d.append(temp[-2])
@@ -753,7 +711,7 @@ def readUBC_DC3Dobs(fileName, dtype = 'DC'):
 
         count = count -1
 
-        # Reach the end of transmitter block, append the src, rx and continue
+        # Reach the end of transmitter block
         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,6 +237,39 @@ 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):
 
@@ -250,33 +283,30 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 #         self.m_prev = self.invProb.m_current
 #         return mref
 
-class Update_IRLS(InversionDirective):
+class update_IRLS(InversionDirective):
 
     eps_min = None
     factor = None
     gamma = None
     phi_m_last = None
     phi_d_last = None
-
-
+    
     def initialize(self):
 
         # Scale the regularization for changes in norm
         if getattr(self, 'phi_m_last', None) is not None:
-
-            self.reg.curModel = self.invProb.curModel
             self.reg.gamma = 1.
             phim_new = self.reg.eval(self.invProb.curModel)
             self.gamma = self.phi_m_last / phim_new
 
         self.reg.curModel = self.invProb.curModel
         self.reg.gamma = self.gamma
-
+        
         if getattr(self, 'phi_d_last', None) is None:
             self.phi_d_last = self.invProb.phi_d
-
+            
     def endIter(self):
-        # Cool the threshold parameter if required
+        # Cool the threshold parameter
         if getattr(self, 'factor', None) is not None:
             eps = self.reg.eps / self.factor
 
@@ -287,74 +317,28 @@ class Update_IRLS(InversionDirective):
 
         # Get phi_m at the end of current iteration
         self.phi_m_last = self.invProb.phi_m_last
-
+        
          # Update the model used for the IRLS weights
         self.reg.curModel = self.invProb.curModel
 
-        # Temporarely set gamma to 1. to get raw phi_m
+        # 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.)
+        self.opt.approxHinv = PC
+         
+        # Temporarely set gamma to 1.
         self.reg.gamma = 1.
-
-        # Compute new model objective function value
+        
+        # Compute change in model objective function and update scaling
         phim_new = self.reg.eval(self.invProb.curModel)
-
-        # Update gamma to scale the regularization between IRLS iterations
         self.reg.gamma = self.phi_m_last / phim_new
-
-        # Set the weighting matrix to None so that it is recomputed next time
-        # it is called in the inversion
-        self.reg._W = None
-
-class Update_lin_PreCond(InversionDirective):
-
-
-    def endIter(self):
-        # Cool the threshold parameter
-
-        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.)
-            self.opt.approxHinv = PC
-
-
-class Update_Wj(InversionDirective):
-    """
-        Create approx-sensitivity base weighting using the probing method
-    """
-    k = None # Number of probing cycles
-    itr = None # Iteration number to update Wj, or always update if None
-
-    def endIter(self):
-
-        if self.itr is None or self.itr == self.opt.iter:
-
-            m = self.invProb.curModel
-            if self.k is None:
-                self.k = int(self.survey.nD/10)
-
-            def JtJv(v):
-
-                Jv = self.prob.Jvec(m, v)
-
-                return self.prob.Jtvec(m,Jv)
-
-            JtJdiag = Utils.diagEst(JtJv,len(m),k=self.k)
-            JtJdiag = JtJdiag / max(JtJdiag)
-
-            self.reg.wght = JtJdiag
-
-class Scale_Beta(InversionDirective):
-    """
-        Instead of a linear cooling schedule, beta is allowed to change based
-        on the ratio between the target misfit and the current data misfit. The
-        update is done only if the misfit is outside some threshold bounds.
-    """
-    tol = 0.05
-
-    def endIter(self):
-
-        # Check if misfit is within the tolerance, otherwise adjust beta
-        val = self.invProb.phi_d / (self.survey.nD*0.5)
-
-        if np.abs(1.-val) > self.tol:
-            self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
+        
+        # TO DO: Re-scale beta if too much change in misfit
+        self.invProb.beta = self.invProb.beta * self.phi_d_last / self.invProb.phi_d
+        
+#==============================================================================
+#          import pylab as plt
+#          plt.figure()
+#          ax = plt.subplot(221)
+#          self.prob.mesh.plotSlice(self.invProb.curModel, ax = ax, normal = 'Z', ind=-5, clim = (0, 0.005))
+#==============================================================================

SimPEG/Examples/DC_Forward_PseudoSection.py

@@ -7,17 +7,17 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p
         DC Forward Simulation
         =====================
 
-        Forward model two conductive spheres in a half-space and plot a
+        Forward model two conductive spheres in a half-space and plot a 
         pseudo-section. Assumes an infinite line source and measures along the
         center of the spheres.
-
+        
         INPUT:
         loc     = Location of spheres [[x1,y1,z1],[x2,y2,z2]]
         radi    = Radius of spheres [r1,r2]
         param   = Conductivity of background and two spheres [m0,m1,m2]
         stype   = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
         dtype   = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
-        Created by @fourndo
+        Created by @fourndo on Mon Feb 01 19:28:06 2016
 
     """
 

SimPEG/Examples/Inversion_IRLS.py

@@ -1,132 +0,0 @@
-from SimPEG import *
-
-
-def run(N=200, plotIt=True):
-    """
-        Inversion: Linear Problem
-        =========================
-
-        Here we go over the basics of creating a linear problem and inversion.
-
-    """
-
-
-    np.random.seed(1)
-
-    std_noise = 1e-2
-
-    mesh = Mesh.TensorMesh([N])
-
-    m0 = np.ones(mesh.nC) * 1e-4
-    nk = 10
-    jk = np.linspace(1.,nk,nk)
-    p = -2.
-    q = 1.
-
-    g = lambda k: np.exp(p*jk[k]*mesh.vectorCCx)*np.cos(np.pi*q*jk[k]*mesh.vectorCCx)
-
-    G = np.empty((nk, mesh.nC))
-
-    for i in range(nk):
-        G[i,:] = g(i)
-
-    mtrue = np.zeros(mesh.nC)
-    mtrue[mesh.vectorCCx > 0.3] = 1.
-    mtrue[mesh.vectorCCx > 0.45] = -0.5
-    mtrue[mesh.vectorCCx > 0.6] = 0
-
-
-    prob = Problem.LinearProblem(mesh, G)
-    survey = Survey.LinearSurvey()
-    survey.pair(prob)
-    survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk)
-    #survey.makeSyntheticData(mtrue, std=std_noise)
-
-    wd = np.ones(nk) * std_noise
-
-    #print survey.std[0]
-    #M = prob.mesh
-    # Distance weighting
-    wr = np.sum(prob.G**2.,axis=0)**0.5
-    wr = ( wr/np.max(wr) )
-
-    reg = Regularization.Simple(mesh)
-    reg.wght = wr
-
-    dmis = DataMisfit.l2_DataMisfit(survey)
-    dmis.Wd = 1./wd
-
-    opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
-    invProb.curModel = m0
-
-    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
-    target = Directives.TargetMisfit()
-
-    betaest = Directives.BetaEstimate_ByEig()
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
-
-
-    mrec = inv.run(m0)
-    ml2 = mrec
-    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
-
-    # Switch regularization to sparse
-    phim = invProb.phi_m_last
-    phid =  invProb.phi_d
-
-    reg = Regularization.Sparse(mesh)
-
-#==============================================================================
-#     fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-#     dmdx = reg.mesh.cellDiffxStencil * mrec
-#     plt.plot(np.sort(dmdx))
-#==============================================================================
-
-    #reg.recModel = mrec
-    reg.wght = np.ones(mesh.nC)
-    reg.mref = np.zeros(mesh.nC)
-    reg.eps_p = 5e-2
-    reg.eps_q = 1e-2
-    reg.norms   = [0., 0., 2., 2.]
-    reg.wght = wr
-
-    opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.)
-    beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
-    #betaest = Directives.BetaEstimate_ByEig()
-    target = Directives.TargetMisfit()
-    IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
-
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
-
-    m0 = mrec
-
-    # Run inversion
-    mrec = inv.run(m0)
-
-    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
-
-
-    if plotIt:
-        import matplotlib.pyplot as plt
-
-        fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-        for i in range(prob.G.shape[0]):
-            axes[0].plot(prob.G[i,:])
-        axes[0].set_title('Columns of matrix G')
-
-        axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
-        axes[1].plot(mesh.vectorCCx, ml2, 'r-')
-        #axes[1].legend(('True Model', 'Recovered Model'))
-        axes[1].set_ylim(-1.0,1.25)
-
-        axes[1].plot(mesh.vectorCCx, mrec, 'k-',lw = 2)
-        axes[1].legend(('True Model', 'Smooth l2-l2',
-        'Sparse lp:' + str(reg.norms[0]) + ', lqx:' + str(reg.norms[1]) ), fontsize = 12)
-        plt.show()
-
-    return prob, survey, mesh, mrec
-
-if __name__ == '__main__':
-    run()

SimPEG/Examples/__init__.py

@@ -9,7 +9,6 @@ import EM_Schenkel_Morrison_Casing
 import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
 import Forward_BasicDirectCurrent
-import Inversion_IRLS
 import Inversion_Linear
 import Mesh_Basic_PlotImage
 import Mesh_Basic_Types
@@ -21,7 +20,7 @@ import Mesh_Tensor_Creation
 import MT_1D_ForwardAndInversion
 import MT_3D_Foward
 
-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_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", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Inversion.py

@@ -33,9 +33,7 @@ class BaseInversion(object):
         self._directiveList = value
         self._directiveList.inversion = self
 
-    def __init__(self, invProb, directiveList=None, **kwargs):
-        if directiveList is None:
-            directiveList = []
+    def __init__(self, invProb, directiveList=[], **kwargs):
         self.directiveList = directiveList
         Utils.setKwargs(self, **kwargs)
 

SimPEG/MT/Utils/ediFilesUtils.py

@@ -7,16 +7,17 @@ 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)
 
@@ -25,8 +26,8 @@ class EDIimporter:
     comps = None
 
     # Hidden properties
-    _outEPSG = None # Project info
-    _2out = None    # The projection operator
+    _outEPSG = None
+    _2out = None
 
 
     def __init__(self, EDIfilesList, compList=None, outEPSG=None):
@@ -112,12 +113,6 @@ class EDIimporter:
     # nOutData=length(obj.data);
     # obj.data(nOutData+1:nOutData+length(TEMP.data),:) = TEMP.data;
     def _transfromPoints(self,longD,latD):
-        # Import the coordinate projections
-        try:
-            import osr
-        except ImportError as e:
-            print 'Could not import osr, missing the gdal package\nCan not project coordinates'
-            raise e
         # Coordinates convertor
         if self._2out is None:
             src = osr.SpatialReference()

SimPEG/Maps.py

@@ -759,29 +759,15 @@ class PolyMap(IdentityMap):
 
             m = [\sigma_1, \sigma_2, c]
 
-        Can take in an actInd vector to account for topography.
-
     """
-    def __init__(self, mesh, order, logSigma=True, normal='X', actInd = None):
+    def __init__(self, mesh, order, logSigma=True, normal='X'):
         IdentityMap.__init__(self, mesh)
         self.logSigma = logSigma
         self.order = order
         self.normal = normal
-        self.actInd = actInd
-
-        if getattr(self, 'actInd', None) is None:
-            self.actInd = range(self.mesh.nC)
-            self.nC = self.mesh.nC
-
-        else:
-            self.nC = len(self.actInd)
 
     slope = 1e4
 
-    @property
-    def shape(self):
-        return (self.nC, self.nP)
-
     @property
     def nP(self):
         if np.isscalar(self.order):
@@ -799,8 +785,8 @@ class PolyMap(IdentityMap):
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
         #2D
         if self.mesh.dim == 2:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
+            X = self.mesh.gridCC[:,0]
+            Y = self.mesh.gridCC[:,1]
             if self.normal =='X':
                 f = polynomial.polyval(Y, c) - X
             elif self.normal =='Y':
@@ -809,9 +795,9 @@ class PolyMap(IdentityMap):
                 raise(Exception("Input for normal = X or Y or Z"))
         #3D
         elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
-            Z = self.mesh.gridCC[self.actInd,2]
+            X = self.mesh.gridCC[:,0]
+            Y = self.mesh.gridCC[:,1]
+            Z = self.mesh.gridCC[:,2]
             if self.normal =='X':
                 f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
             elif self.normal =='Y':
@@ -820,7 +806,6 @@ class PolyMap(IdentityMap):
                 f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
-
         else:
             raise(Exception("Only supports 2D"))
 
@@ -834,8 +819,8 @@ class PolyMap(IdentityMap):
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
         #2D
         if self.mesh.dim == 2:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
+            X = self.mesh.gridCC[:,0]
+            Y = self.mesh.gridCC[:,1]
 
             if self.normal =='X':
                 f = polynomial.polyval(Y, c) - X
@@ -847,9 +832,9 @@ class PolyMap(IdentityMap):
                 raise(Exception("Input for normal = X or Y or Z"))
         #3D
         elif self.mesh.dim == 3:
-            X = self.mesh.gridCC[self.actInd,0]
-            Y = self.mesh.gridCC[self.actInd,1]
-            Z = self.mesh.gridCC[self.actInd,2]
+            X = self.mesh.gridCC[:,0]
+            Y = self.mesh.gridCC[:,1]
+            Z = self.mesh.gridCC[:,2]
 
             if self.normal =='X':
                 f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X

SimPEG/Mesh/DiffOperators.py

@@ -584,7 +584,67 @@ class DiffOperators(object):
 
         return Pbc, Pin, Pout
 
+    def getBCProjWF_simple(self, discretization='CC'):
+        """
 
+        The weak form boundary condition projection matrices
+        when mixed boundary condition is used
+
+
+        """
+
+        if discretization is not 'CC':
+            raise NotImplementedError('Boundary conditions only implemented for CC discretization.')
+
+        def projBC(n):
+            ij = ([0,n], [0,1])
+            vals = [0,0]
+            vals[0] = 1
+            vals[1] = 1
+            return sp.csr_matrix((vals, ij), shape=(n+1,2))
+
+        def projDirichlet(n, bc):
+            bc = checkBC(bc)
+            ij = ([0,n], [0,1])
+            vals = [0,0]
+            if(bc[0] == 'dirichlet'):
+                vals[0] = -1
+            if(bc[1] == 'dirichlet'):
+                vals[1] = 1
+            return sp.csr_matrix((vals, ij), shape=(n+1,2))
+
+        BC = [['dirichlet','dirichlet'],['dirichlet','dirichlet'],['dirichlet','dirichlet']]
+        n = self.vnC
+        indF = self.faceBoundaryInd
+        if(self.dim == 1):
+            Pbc = projDirichlet(n[0], BC[0])
+            B = projBC(n[0])
+            indF = indF[0] | indF[1]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+        elif(self.dim == 2):
+            Pbc1 = sp.kron(speye(n[1]), projDirichlet(n[0], BC[0]))
+            Pbc2 = sp.kron(projDirichlet(n[1], BC[1]), speye(n[0]))
+            Pbc = sp.block_diag((Pbc1, Pbc2), format="csr")
+            B1 = sp.kron(speye(n[1]), projBC(n[0]))
+            B2 = sp.kron(projBC(n[1]), speye(n[0]))
+            B = sp.block_diag((B1, B2), format="csr")
+            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3])]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+        elif(self.dim == 3):
+            Pbc1 = kron3(speye(n[2]), speye(n[1]), projDirichlet(n[0], BC[0]))
+            Pbc2 = kron3(speye(n[2]), projDirichlet(n[1], BC[1]), speye(n[0]))
+            Pbc3 = kron3(projDirichlet(n[2], BC[2]), speye(n[1]), speye(n[0]))
+            Pbc = sp.block_diag((Pbc1, Pbc2, Pbc3), format="csr")
+            B1 = kron3(speye(n[2]), speye(n[1]), projBC(n[0]))
+            B2 = kron3(speye(n[2]), projBC(n[1]), speye(n[0]))
+            B3 = kron3(projBC(n[2]), speye(n[1]), speye(n[0]))
+            B = sp.block_diag((B1, B2, B3), format="csr")
+            indF = np.r_[(indF[0] | indF[1]), (indF[2] | indF[3]), (indF[4] | indF[5])]
+            Pbc = Pbc*sdiag(self.area[indF])
+
+        return Pbc, B.T
     # --------------- Averaging ---------------------
 
     @property

SimPEG/Mesh/MeshIO.py

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

SimPEG/Mesh/TreeMesh.py

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

SimPEG/Mesh/View.py

@@ -42,9 +42,9 @@ class TensorView(object):
 
     def plotImage(self, v, vType='CC', grid=False, view='real',
               ax=None, clim=None, showIt=False,
-              pcolorOpts=None,
-              streamOpts=None,
-              gridOpts=None,
+              pcolorOpts={},
+              streamOpts={'color':'k'},
+              gridOpts={'color':'k'},
               numbering=True, annotationColor='w'
               ):
         """
@@ -84,12 +84,6 @@ class TensorView(object):
             M.plotImage(v, annotationColor='k', showIt=True)
 
         """
-        if pcolorOpts is None:
-            pcolorOpts = {}
-        if streamOpts is None:
-            streamOpts = {'color':'k'}
-        if gridOpts is None:
-            gridOpts = {'color':'k'}
 
         if ax is None:
             fig = plt.figure()
@@ -180,9 +174,9 @@ class TensorView(object):
     def plotSlice(self, v, vType='CC',
                   normal='Z', ind=None, grid=False, view='real',
                   ax=None, clim=None, showIt=False,
-                  pcolorOpts=None,
-                  streamOpts=None,
-                  gridOpts=None
+                  pcolorOpts={},
+                  streamOpts={'color':'k'},
+                  gridOpts={'color':'k', 'alpha':0.5}
                   ):
 
         """
@@ -203,12 +197,6 @@ class TensorView(object):
             M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8})
 
         """
-        if pcolorOpts is None:
-            pcolorOpts = {}
-        if streamOpts is None:
-            streamOpts = {'color':'k'}
-        if gridOpts is None:
-            gridOpts = {'color':'k', 'alpha':0.5}
         if type(vType) in [list, tuple]:
             assert ax is None, "cannot specify an axis to plot on with this function."
             fig, axs = plt.subplots(1,len(vType))
@@ -301,17 +289,11 @@ class TensorView(object):
 
     def _plotImage2D(self, v, vType='CC', grid=False, view='real',
               ax=None, clim=None, showIt=False,
-              pcolorOpts=None,
-              streamOpts=None,
-              gridOpts=None
+              pcolorOpts={},
+              streamOpts={'color':'k'},
+              gridOpts={'color':'k'}
               ):
 
-        if pcolorOpts is None:
-            pcolorOpts = {}
-        if streamOpts is None:
-            streamOpts = {'color':'k'}
-        if gridOpts is None:
-            gridOpts = {'color':'k'}
         vTypeOptsCC = ['N','CC','Fx','Fy','Ex','Ey']
         vTypeOptsV = ['CCv','F','E']
         vTypeOpts = vTypeOptsCC + vTypeOptsV

SimPEG/Optimization.py

@@ -1003,7 +1003,6 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
                 # perturb inactive set off of bounds so that they are included in the step
                 delx = delx + self.stepOffBoundsFact * (rhs_a * dm_i / dm_a)
 
-
             # Only keep gradients going in the right direction on the active set
             indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
             delx[indx] = 0.

SimPEG/Regularization.py

@@ -344,6 +344,7 @@ class BaseRegularization(object):
     def W(self):
         """Full regularization weighting matrix W."""
         return sp.identity(self.regmesh.nC)
+        # self.regmesh._Pac.T * sp.identity(self.regmesh.nC) * self.regmesh._Pac # or do we want sp.identity(self.mesh.nC) or even just Utils.Identity() ?
 
     @Utils.timeIt
     def eval(self, m):
@@ -374,12 +375,11 @@ class BaseRegularization(object):
     @Utils.timeIt
     def eval2Deriv(self, m, v=None):
         """
-        Second derivative
 
-        :param numpy.array m: geophysical model
-        :param numpy.array v: vector to multiply
-        :rtype: scipy.sparse.csr_matrix or numpy.ndarray
-        :return: WtW or WtW*v
+            :param numpy.array m: geophysical model
+            :param numpy.array v: vector to multiply
+            :rtype: scipy.sparse.csr_matrix or numpy.ndarray
+            :return: WtW or WtW*v
 
         The regularization is:
 
@@ -402,48 +402,25 @@ class BaseRegularization(object):
 
 class Tikhonov(BaseRegularization):
     """
-    L2 Tikhonov regularization with both smallness and smoothness (first order
-    derivative) contributions.
-
-    .. math::
-        \phi_m(\mathbf{m}) = \\alpha_s \| W_s (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
-        + \\alpha_x \| W_x \\frac{\partial}{\partial x} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
-        + \\alpha_y \| W_y \\frac{\partial}{\partial y} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
-        + \\alpha_z \| W_z \\frac{\partial}{\partial z} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
-
-    Note if the key word argument `mrefInSmooth` is False, then mref is not
-    included in the smoothness contribution.
-
-    :param Mesh mesh: SimPEG mesh
-    :param Maps mapping: regularization mapping, takes the model from model space to the thing you want to regularize
-    :param numpy.ndarray indActive: active cell indices for reducing the size of differential operators in the definition of a regularization mesh
-    :param bool mrefInSmooth: (default = False) put mref in the smoothness component?
-    :param float alpha_s: (default 1e-6) smallness weight
-    :param float alpha_x: (default 1) smoothness weight for first derivative in the x-direction
-    :param float alpha_y: (default 1) smoothness weight for first derivative in the y-direction
-    :param float alpha_z: (default 1) smoothness weight for first derivative in the z-direction
-    :param float alpha_xx: (default 1) smoothness weight for second derivative in the x-direction
-    :param float alpha_yy: (default 1) smoothness weight for second derivative in the y-direction
-    :param float alpha_zz: (default 1) smoothness weight for second derivative in the z-direction
     """
-    mrefInSmooth = False  # put mref in the smoothness contribution
-    alpha_s      = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Wsmall'], "Smallness weight")
-    alpha_x      = Utils.dependentProperty('_alpha_x',  1.0, ['_W', '_Wx'],     "Weight for the first derivative in the x direction")
-    alpha_y      = Utils.dependentProperty('_alpha_y',  1.0, ['_W', '_Wy'],     "Weight for the first derivative in the y direction")
-    alpha_z      = Utils.dependentProperty('_alpha_z',  1.0, ['_W', '_Wz'],     "Weight for the first derivative in the z direction")
-    alpha_xx     = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'],    "Weight for the second derivative in the x direction")
-    alpha_yy     = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'],    "Weight for the second derivative in the y direction")
-    alpha_zz     = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'],    "Weight for the second derivative in the z direction")
+    mrefInSmooth = True  #: SMOOTH and SMOOTH_MOD_DIF options
+    alpha_s  = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight")
+    alpha_x  = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
+    alpha_y  = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
+    alpha_z  = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
+    alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
+    alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
+    alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
 
     def __init__(self, mesh, mapping=None, indActive = None, **kwargs):
         BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
 
     @property
-    def Wsmall(self):
-        """Regularization matrix Wsmall"""
-        if getattr(self,'_Wsmall', None) is None:
-            self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5)
-        return self._Wsmall
+    def Ws(self):
+        """Regularization matrix Ws"""
+        if getattr(self,'_Ws', None) is None:
+            self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5)
+        return self._Ws
 
     @property
     def Wx(self):
@@ -506,44 +483,25 @@ class Tikhonov(BaseRegularization):
     def W(self):
         """Full regularization matrix W"""
         if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
+            wlist = (self.Ws, 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)
-
     @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) )
+            r1 = self.Wsmooth * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            return 0.5*(r1.dot(r1)+r2.dot(r2))
         elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m )
-            return r.T * ( self.Wsmooth * self.mapping.deriv(m) )
+            r = self.W * ( self.mapping * (m - self.mref) )
+            return 0.5*r.dot(r)
+
 
     @Utils.timeIt
     def evalDeriv(self, m):
         """
+
         The regularization is:
 
         .. math::
@@ -557,33 +515,45 @@ class Tikhonov(BaseRegularization):
             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
-        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+        if self.mrefInSmooth == True:
+            mD1 = self.mapping.deriv(m)
+            mD2 = self.mapping.deriv(m - self.mref)
+            r1 = self.Wsmooth * ( self.mapping * (m))
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            out1 = mD1.T * ( self.Wsmooth.T * r1 )
+            out2 = mD2.T * ( self.Ws.T * r2 )
+            out = out1+out2
+        elif self.mrefInSmooth == False:
+            mD = self.mapping.deriv(m - self.mref)
+            r = self.W * ( self.mapping * (m - self.mref) )
+            out = mD.T * ( self.W.T * r )
+        return out
 
 
-class Simple(Tikhonov):
+class Simple(BaseRegularization):
     """
-    Simple regularization that does not include length scales in the derivatives.
+        Only for tensor mesh
     """
 
-    mrefInSmooth = False  #: SMOOTH and SMOOTH_MOD_DIF options
-    alpha_s     = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
+    mrefInSmooth = True  #: SMOOTH and SMOOTH_MOD_DIF options
+    alpha_s     = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight")
     alpha_x     = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
     alpha_y     = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
     alpha_z     = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
     wght        = 1.
-
+    
     def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
         BaseRegularization.__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,'_Wsmall', None) is None:
-            self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
-        return self._Wsmall
+    def Ws(self):
+        """Regularization matrix Ws"""
+        if getattr(self,'_Ws', None) is None:
+                self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
+        return self._Ws
 
     @property
     def Wx(self):
@@ -622,65 +592,99 @@ class Simple(Tikhonov):
     def W(self):
         """Full regularization matrix W"""
         if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
+            wlist = (self.Ws, 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 eval(self, m):
         if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
+            r1 = self.Wsmooth * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            return 0.5*(r1.dot(r1)+r2.dot(r2))
         elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m)
-        return 0.5 * r.dot(r)
+            r = self.W * ( self.mapping * (m - self.mref) )
+            return 0.5*r.dot(r)
+        return phim
+
+
+
+    @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})}
+
+        """
+        if self.mrefInSmooth == True:
+            mD1 = self.mapping.deriv(m)
+            mD2 = self.mapping.deriv(m - self.mref)
+            r1 = self.Wsmooth * ( self.mapping * (m))
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            out1 = mD1.T * ( self.Wsmooth.T * r1 )
+            out2 = mD2.T * ( self.Ws.T * r2 )
+            out = out1+out2
+        elif self.mrefInSmooth == False:
+            mD = self.mapping.deriv(m - self.mref)
+            r = self.W * ( self.mapping * (m - self.mref) )
+            out = mD.T * ( self.W.T * r )
+        return out
 
 
 class Sparse(Simple):
 
     # set default values
-    eps_p      = 1e-1
-    eps_q      = 1e-1
+    eps      = 1e-1
     curModel = None # use a model to compute the weights
     gamma    = 1.
-    norms    = [0., 2., 2., 2.]
+    p        = 0.
+    qx       = 2.
+    qy       = 2.
+    qz       = 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"""
+    def Ws(self):
+        """Regularization matrix Ws"""
         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])
+            self.rs = self.R(f_m , self.p)
             #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 = self.R( f_m , self.qx)
             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
@@ -694,9 +698,9 @@ class Sparse(Simple):
 
         else:
             f_m = self.regmesh.cellDiffyStencil * self.curModel
-            self.ry = self.R( f_m , self.eps_q, self.norms[2])
+            self.ry = self.R( f_m , self.qy)
             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
@@ -708,7 +712,7 @@ class Sparse(Simple):
 
         else:
             f_m = self.regmesh.cellDiffzStencil * self.curModel
-            self.rz = self.R( f_m , self.eps_q, self.norms[3])
+            self.rz = self.R( f_m , self.qz)
             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
@@ -728,14 +732,14 @@ 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 self._W
+        #if getattr(self, '_W', None) is None:
+        wlist = (self.Ws, self.Wsmooth)
+        #self._W = sp.vstack(wlist)
+        return sp.vstack(wlist)
+        
+    def R(self, f_m , exponent):
 
-    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 = (self.eps**(1-exponent/2.))**0.5
+        r = eta / (f_m**2.+self.eps**2.)**((1-exponent/2.)/2.)
 
         return r

SimPEG/Utils/ModelBuilder.py

@@ -88,14 +88,12 @@ def getIndicesBlock(p0,p1,ccMesh):
     # Return a tuple
     return ind
 
-def defineBlock(ccMesh,p0,p1,vals=None):
+def defineBlock(ccMesh,p0,p1,vals=[0,1]):
     """
         Build a block with the conductivity specified by condVal.  Returns an array.
         vals[0]  conductivity of the block
         vals[1]  conductivity of the ground
     """
-    if vals is None:
-        vals = [0,1]
     sigma = np.zeros(ccMesh.shape[0]) + vals[1]
     ind   = getIndicesBlock(p0,p1,ccMesh)
 
@@ -103,11 +101,7 @@ def defineBlock(ccMesh,p0,p1,vals=None):
 
     return mkvc(sigma)
 
-def defineElipse(ccMesh, center=None, anisotropy=None, slope=10., theta=0.):
-    if center is None:
-        center = [0,0,0]
-    if anisotropy is None:
-        anisotropy = [1,1,1]
+def defineElipse(ccMesh, center=[0,0,0], anisotropy=[1,1,1], slope=10., theta=0.):
     G = ccMesh.copy()
     dim = ccMesh.shape[1]
     for i in range(dim):
@@ -162,7 +156,7 @@ def getIndicesSphere(center,radius,ccMesh):
     # Return a tuple
     return ind
 
-def defineTwoLayers(ccMesh,depth,vals=None):
+def defineTwoLayers(ccMesh,depth,vals=[0,1]):
     """
     Define a two layered model.  Depth of the first layer must be specified.
     CondVals vector with the conductivity values of the layers.  Eg:
@@ -173,8 +167,6 @@ def defineTwoLayers(ccMesh,depth,vals=None):
         0                          depth                                 zf
              1st layer                       2nd layer
     """
-    if vals is None:
-        vals = [0,1]
     sigma = np.zeros(ccMesh.shape[0]) + vals[1]
 
     dim = np.size(ccMesh[0,:])
@@ -260,7 +252,7 @@ def layeredModel(ccMesh, layerTops, layerValues):
 
 
 
-def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=None):
+def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
     """
         Create a random model by convolving a kernel with a
         uniformly distributed model.
@@ -284,8 +276,6 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=None):
 
 
     """
-    if bounds is None:
-        bounds = [0,1]
 
     if seed is None:
         seed = np.random.randint(1e3)

SimPEG/Utils/codeutils.py

@@ -55,10 +55,8 @@ def hook(obj, method, name=None, overwrite=False, silent=False):
         print 'Method '+name+' was not overwritten.'
 
 
-def setKwargs(obj, ignore=None,  **kwargs):
+def setKwargs(obj, ignore=[],  **kwargs):
     """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
-    if ignore is None:
-        ignore = []
     for attr in kwargs:
         if attr in ignore:
             continue

SimPEG/Utils/io_utils.py

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

docs/examples/DC_Forward_PseudoSection.rst

@@ -12,16 +12,8 @@
 DC Forward Simulation
 =====================
 
-Forward model two conductive spheres in a half-space and plot a
-pseudo-section. Assumes an infinite line source and measures along the
-center of the spheres.
+Forward model conductive spheres in a half-space and plot a pseudo-section
 
-INPUT:
-loc     = Location of spheres [[x1,y1,z1],[x2,y2,z2]]
-radi    = Radius of spheres [r1,r2]
-param   = Conductivity of background and two spheres [m0,m1,m2]
-stype   = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole)
-dtype   = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential)
 Created by @fourndo on Mon Feb 01 19:28:06 2016
 
 

docs/examples/Inversion_IRLS.rst

@@ -1,26 +0,0 @@
-.. _examples_Inversion_IRLS:
-
-.. --------------------------------- ..
-..                                   ..
-..    THIS FILE IS AUTO GENEREATED   ..
-..                                   ..
-..    SimPEG/Examples/__init__.py    ..
-..                                   ..
-.. --------------------------------- ..
-
-
-Inversion: Linear Problem
-=========================
-
-Here we go over the basics of creating a linear problem and inversion.
-
-
-
-.. plot::
-
-    from SimPEG import Examples
-    Examples.Inversion_IRLS.run()
-
-.. literalinclude:: ../../SimPEG/Examples/Inversion_IRLS.py
-    :language: python
-    :linenos:

tests/em/tdem/test_TDEM_forward_Analytic.py

@@ -10,9 +10,7 @@ except ImportError, e:
     MumpsSolver = SolverLU
 
 
-def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None, showIt=False):
-    if bounds is None:
-        bounds = [1e-5,1e-3]
+def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,1e-3], showIt=False):
     if meshType == 'CYL':
         cs, ncx, ncz, npad = 5., 30, 10, 15
         hx = [(cs,ncx), (cs,npad,1.3)]

tests/mesh/test_Mixed_boundaryPoisson.py

@@ -0,0 +1,396 @@
+import numpy as np
+import scipy.sparse as sp
+import unittest
+import matplotlib.pyplot as plt
+from SimPEG import *
+
+MESHTYPES = ['uniformTensorMesh']
+
+def getxBCyBC(mesh, alpha, beta, gamma):
+# def getxBCyBC(mesh, alpha, beta, gamma):
+    """
+
+    """
+    if mesh.dim == 1: #1D
+        if (len(alpha) != 2 or len(beta) != 2 or len(gamma) != 2):
+            raise Exception("Lenght of list, alpha should be 2")
+        fCCxm,fCCxp = mesh.cellBoundaryInd
+        nBC = fCCxm.sum()+fCCxp.sum()
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+
+        xBC = np.r_[xBC_xm, xBC_xp]
+        yBC = np.r_[yBC_xm, yBC_xp]
+
+    elif mesh.dim == 2: #2D
+        if (len(alpha) != 4 or len(beta) != 4 or len(gamma) != 4):
+            raise Exception("Lenght of list, alpha should be 4")
+
+        fCCxm,fCCxp,fCCym,fCCyp = mesh.cellBoundaryInd
+        fxm,fxp,fym,fyp = mesh.faceBoundaryInd
+        nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum()
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+        h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp]
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
+        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
+
+        h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
+        h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
+        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
+        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
+        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+        xBC_ym = 0.5*a_ym
+        xBC_yp = 0.5*a_yp/b_yp
+        yBC_ym = 0.5*(1.-b_ym)
+        yBC_yp = 0.5*(1.-1./b_yp)
+
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
+
+        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
+        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
+        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
+        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
+
+        xBC = np.r_[xBC_x, xBC_y]
+        yBC = np.r_[yBC_x, yBC_y]
+
+    elif mesh.dim == 3: #3D
+        if (len(alpha) != 6 or len(beta) != 6 or len(gamma) != 6):
+            raise Exception("Lenght of list, alpha should be 6")
+        fCCxm,fCCxp,fCCym,fCCyp,fCCzm,fCCzp = mesh.cellBoundaryInd
+        fxm,fxp,fym,fyp,fzm,fzp = mesh.faceBoundaryInd
+        nBC = fCCxm.sum()+fCCxp.sum()+fCCxm.sum()+fCCxp.sum()
+        h_xm, h_xp = mesh.gridCC[fCCxm], mesh.gridCC[fCCxp]
+        h_ym, h_yp = mesh.gridCC[fCCym], mesh.gridCC[fCCyp]
+        h_zm, h_zp = mesh.gridCC[fCCzm], mesh.gridCC[fCCzp]
+
+        alpha_xm, beta_xm, gamma_xm = alpha[0], beta[0], gamma[0]
+        alpha_xp, beta_xp, gamma_xp = alpha[1], beta[1], gamma[1]
+        alpha_ym, beta_ym, gamma_ym = alpha[2], beta[2], gamma[2]
+        alpha_yp, beta_yp, gamma_yp = alpha[3], beta[3], gamma[3]
+        alpha_zm, beta_zm, gamma_zm = alpha[2], beta[2], gamma[2]
+        alpha_zp, beta_zp, gamma_zp = alpha[3], beta[3], gamma[3]
+
+        h_xm, h_xp = mesh.gridCC[fCCxm,0], mesh.gridCC[fCCxp,0]
+        h_ym, h_yp = mesh.gridCC[fCCym,1], mesh.gridCC[fCCyp,1]
+        h_zm, h_zp = mesh.gridCC[fCCzm,2], mesh.gridCC[fCCzp,2]
+
+        a_xm = gamma_xm/(0.5*alpha_xm-beta_xm/h_xm)
+        b_xm = (0.5*alpha_xm+beta_xm/h_xm)/(0.5*alpha_xm-beta_xm/h_xm)
+        a_xp = gamma_xp/(0.5*alpha_xp-beta_xp/h_xp)
+        b_xp = (0.5*alpha_xp+beta_xp/h_xp)/(0.5*alpha_xp-beta_xp/h_xp)
+
+        a_ym = gamma_ym/(0.5*alpha_ym-beta_ym/h_ym)
+        b_ym = (0.5*alpha_ym+beta_ym/h_ym)/(0.5*alpha_ym-beta_ym/h_ym)
+        a_yp = gamma_yp/(0.5*alpha_yp-beta_yp/h_yp)
+        b_yp = (0.5*alpha_yp+beta_yp/h_yp)/(0.5*alpha_yp-beta_yp/h_yp)
+
+        a_zm = gamma_zm/(0.5*alpha_zm-beta_zm/h_zm)
+        b_zm = (0.5*alpha_zm+beta_zm/h_zm)/(0.5*alpha_zm-beta_zm/h_zm)
+        a_zp = gamma_zp/(0.5*alpha_zp-beta_zp/h_zp)
+        b_zp = (0.5*alpha_zp+beta_zp/h_zp)/(0.5*alpha_zp-beta_zp/h_zp)
+
+        xBC_xm = 0.5*a_xm
+        xBC_xp = 0.5*a_xp/b_xp
+        yBC_xm = 0.5*(1.-b_xm)
+        yBC_xp = 0.5*(1.-1./b_xp)
+        xBC_ym = 0.5*a_ym
+        xBC_yp = 0.5*a_yp/b_yp
+        yBC_ym = 0.5*(1.-b_ym)
+        yBC_yp = 0.5*(1.-1./b_yp)
+        xBC_zm = 0.5*a_zm
+        xBC_zp = 0.5*a_zp/b_zp
+        yBC_zm = 0.5*(1.-b_zm)
+        yBC_zp = 0.5*(1.-1./b_zp)
+
+        sortindsfx = np.argsort(np.r_[np.arange(mesh.nFx)[fxm], np.arange(mesh.nFx)[fxp]])
+        sortindsfy = np.argsort(np.r_[np.arange(mesh.nFy)[fym], np.arange(mesh.nFy)[fyp]])
+        sortindsfz = np.argsort(np.r_[np.arange(mesh.nFz)[fzm], np.arange(mesh.nFz)[fzp]])
+
+        xBC_x = np.r_[xBC_xm, xBC_xp][sortindsfx]
+        xBC_y = np.r_[xBC_ym, xBC_yp][sortindsfy]
+        xBC_z = np.r_[xBC_zm, xBC_zp][sortindsfz]
+
+        yBC_x = np.r_[yBC_xm, yBC_xp][sortindsfx]
+        yBC_y = np.r_[yBC_ym, yBC_yp][sortindsfy]
+        yBC_z = np.r_[yBC_zm, yBC_zp][sortindsfz]
+
+        xBC = np.r_[xBC_x, xBC_y, xBC_z]
+        yBC = np.r_[yBC_x, yBC_y, yBC_z]
+
+    return xBC, yBC
+
+class Test1D_InhomogeneousMixed(Tests.OrderTest):
+    name = "1D - Mixed"
+    meshTypes = MESHTYPES
+    meshDimension = 1
+    expectedOrders = 2
+    meshSizes = [4, 8, 16, 32]
+
+    def getError(self):
+        #Test function
+        phi_fun = lambda x: np.cos(np.pi*x)
+        j_fun = lambda x: np.pi*np.sin(np.pi*x)
+        phi_deriv = lambda x: -j_fun(x)
+        q_fun = lambda x: (np.pi**2)*np.cos(np.pi*x)
+
+        xc_ana = phi_fun(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        j_ana = j_fun(self.M.gridFx)
+
+        # Get boundary locations
+        vecN = self.M.vectorNx
+        vecC = self.M.vectorCCx
+
+        # Setup Mixed B.C (alpha, beta, gamma)
+        alpha_xm, alpha_xp = 1., 1.
+        beta_xm, beta_xp =  1., 1.
+        alpha = np.r_[alpha_xm, alpha_xp]
+        beta = np.r_[beta_xm, beta_xp]
+        vecN = self.M.vectorNx
+        vecC = self.M.vectorCCx
+        phi_bc = phi_fun(vecN[[0,-1]])
+        phi_deriv_bc = phi_deriv(vecN[[0,-1]])
+        gamma = alpha*phi_bc + beta*phi_deriv_bc
+        x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma)
+
+
+        sigma = np.ones(self.M.nC)
+        Mfrho = self.M.getFaceInnerProduct(1./sigma)
+        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
+        V = Utils.sdiag(self.M.vol)
+        Div = V*self.M.faceDiv
+        P_BC, B = self.M.getBCProjWF_simple()
+        q = q_fun(self.M.gridCC)
+        M = B*self.M.aveCC2F
+        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
+        # Mrhoj = D.T V phi + P_BC*Utils.sdiag(y_BC)*M phi - P_BC*x_BC
+        rhs = V*q + Div*MfrhoI*P_BC*x_BC
+        A = Div*MfrhoI*G
+
+        if self.myTest == 'xc':
+            #TODO: fix the null space
+            Ainv = Solver(A)
+            xc = Ainv*rhs
+            err = np.linalg.norm((xc-xc_ana), np.inf)
+        else:
+            NotImplementedError
+        return err
+
+
+    def test_order(self):
+        print "==== Testing Mixed boudary conduction for CC-problem ===="
+        self.name = "1D"
+        self.myTest = 'xc'
+        self.orderTest()
+
+class Test2D_InhomogeneousMixed(Tests.OrderTest):
+    name = "2D - Mixed"
+    meshTypes = MESHTYPES
+    meshDimension = 2
+    expectedOrders = 2
+    meshSizes = [4, 8, 16, 32]
+
+    def getError(self):
+        #Test function
+        phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
+        j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])
+        j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])
+        phideriv_funX = lambda x: -j_funX(x)
+        phideriv_funY = lambda x: -j_funY(x)
+        q_fun = lambda x: +2*(np.pi**2)*phi_fun(x)
+
+        xc_ana = phi_fun(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        jX_ana = j_funX(self.M.gridFx)
+        jY_ana = j_funY(self.M.gridFy)
+        j_ana = np.r_[jX_ana,jY_ana]
+
+        # Get boundary locations
+        fxm,fxp,fym,fyp = self.M.faceBoundaryInd
+        gBFxm = self.M.gridFx[fxm,:]
+        gBFxp = self.M.gridFx[fxp,:]
+        gBFym = self.M.gridFy[fym,:]
+        gBFyp = self.M.gridFy[fyp,:]
+
+        # Setup Mixed B.C (alpha, beta, gamma)
+        alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        beta_xm, beta_xp =  np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+        beta_ym, beta_yp =  np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+
+        phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
+        phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
+
+        phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
+        phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
+
+        gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
+        gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
+        gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
+        gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
+        gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp)
+
+        alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
+        beta = [beta_xm, beta_xp, beta_ym, beta_yp]
+        gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
+
+        x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma)
+
+
+        sigma = np.ones(self.M.nC)
+        Mfrho = self.M.getFaceInnerProduct(1./sigma)
+        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
+        V = Utils.sdiag(self.M.vol)
+        Div = V*self.M.faceDiv
+        P_BC, B = self.M.getBCProjWF_simple()
+        q = q_fun(self.M.gridCC)
+        M = B*self.M.aveCC2F
+        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
+        rhs = V*q + Div*MfrhoI*P_BC*x_BC
+        A = Div*MfrhoI*G
+
+        if self.myTest == 'xc':
+            Ainv = Solver(A)
+            xc = Ainv*rhs
+            err = np.linalg.norm((xc-xc_ana), np.inf)
+        else:
+            NotImplementedError
+        return err
+
+
+    def test_order(self):
+        print "==== Testing Mixed boudary conduction for CC-problem ===="
+        self.name = "2D"
+        self.myTest = 'xc'
+        self.orderTest()
+
+class Test3D_InhomogeneousMixed(Tests.OrderTest):
+    name = "3D - Mixed"
+    meshTypes = MESHTYPES
+    meshDimension = 3
+    expectedOrders = 2
+    meshSizes = [4, 8, 16]
+
+    def getError(self):
+        #Test function
+        phi_fun = lambda x: np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
+        j_funX = lambda x: +np.pi*np.sin(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
+        j_funY = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.sin(np.pi*x[:,1])*np.cos(np.pi*x[:,2])
+        j_funZ = lambda x: +np.pi*np.cos(np.pi*x[:,0])*np.cos(np.pi*x[:,1])*np.sin(np.pi*x[:,2])
+
+        phideriv_funX = lambda x: -j_funX(x)
+        phideriv_funY = lambda x: -j_funY(x)
+        phideriv_funZ = lambda x: -j_funZ(x)
+
+        q_fun = lambda x: 3*(np.pi**2)*phi_fun(x)
+
+        xc_ana = phi_fun(self.M.gridCC)
+        q_ana = q_fun(self.M.gridCC)
+        jX_ana = j_funX(self.M.gridFx)
+        jY_ana = j_funY(self.M.gridFy)
+        j_ana = np.r_[jX_ana,jY_ana,jY_ana]
+
+        # Get boundary locations
+        fxm,fxp,fym,fyp,fzm,fzp = self.M.faceBoundaryInd
+        gBFxm = self.M.gridFx[fxm,:]
+        gBFxp = self.M.gridFx[fxp,:]
+        gBFym = self.M.gridFy[fym,:]
+        gBFyp = self.M.gridFy[fyp,:]
+        gBFzm = self.M.gridFz[fzm,:]
+        gBFzp = self.M.gridFz[fzp,:]
+
+        # Setup Mixed B.C (alpha, beta, gamma)
+        alpha_xm, alpha_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        beta_xm, beta_xp =  np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
+        alpha_ym, alpha_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+        beta_ym, beta_yp =  np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
+        alpha_zm, alpha_zp = np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1])
+        beta_zm, beta_zp =  np.ones_like(gBFzm[:,1]), np.ones_like(gBFzp[:,1])
+
+
+        phi_bc_xm, phi_bc_xp = phi_fun(gBFxm), phi_fun(gBFxp)
+        phi_bc_ym, phi_bc_yp = phi_fun(gBFym), phi_fun(gBFyp)
+        phi_bc_zm, phi_bc_zp = phi_fun(gBFzm), phi_fun(gBFzp)
+
+        phiderivX_bc_xm, phiderivX_bc_xp = phideriv_funX(gBFxm), phideriv_funX(gBFxp)
+        phiderivY_bc_ym, phiderivY_bc_yp = phideriv_funY(gBFym), phideriv_funY(gBFyp)
+        phiderivY_bc_zm, phiderivY_bc_zp = phideriv_funY(gBFzm), phideriv_funY(gBFzp)
+
+        gamma_fun = lambda alpha, beta, phi, phi_deriv: alpha*phi + beta*phi_deriv
+        gamma_xm = gamma_fun(alpha_xm, beta_xm, phi_bc_xm, phiderivX_bc_xm)
+        gamma_xp = gamma_fun(alpha_xp, beta_xp, phi_bc_xp, phiderivX_bc_xp)
+        gamma_ym = gamma_fun(alpha_ym, beta_ym, phi_bc_ym, phiderivY_bc_ym)
+        gamma_yp = gamma_fun(alpha_yp, beta_yp, phi_bc_yp, phiderivY_bc_yp)
+        gamma_zm = gamma_fun(alpha_zm, beta_zm, phi_bc_zm, phiderivY_bc_zm)
+        gamma_zp = gamma_fun(alpha_zp, beta_zp, phi_bc_zp, phiderivY_bc_zp)
+
+        alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
+        beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
+        gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
+
+        x_BC, y_BC = getxBCyBC(self.M, alpha, beta, gamma)
+
+
+        sigma = np.ones(self.M.nC)
+        Mfrho = self.M.getFaceInnerProduct(1./sigma)
+        MfrhoI = self.M.getFaceInnerProduct(1./sigma, invMat=True)
+        V = Utils.sdiag(self.M.vol)
+        Div = V*self.M.faceDiv
+        P_BC, B = self.M.getBCProjWF_simple()
+        q = q_fun(self.M.gridCC)
+        M = B*self.M.aveCC2F
+        G = Div.T - P_BC*Utils.sdiag(y_BC)*M
+        rhs = V*q + Div*MfrhoI*P_BC*x_BC
+        A = Div*MfrhoI*G
+
+        if self.myTest == 'xc':
+            #TODO: fix the null space
+            Ainv = Solver(A)
+            xc = Ainv*rhs
+            err = np.linalg.norm((xc-xc_ana), np.inf)
+        else:
+            NotImplementedError
+        return err
+
+
+    def test_order(self):
+        print "==== Testing Mixed boudary conduction for CC-problem ===="
+        self.name = "3D"
+        self.myTest = 'xc'
+        self.orderTest()
+
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/mt/test_Problem3D_againstAnalytic.py

@@ -242,6 +242,9 @@ 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))