Initial work on the octree creation example.

Em/dev

SimPEG/DataMisfit.py

@@ -59,20 +59,6 @@ class BaseDataMisfit(object):
         """
         raise NotImplementedError('This method should be overwritten.')
 
-    # TODO: implement target misfit as a property, or possibly as an inversion directive.
-
-    # def target(self, forward):
-    #     """target(forward)
-
-    #         Target for data misfit. By default this is the number of data,
-    #         which satisfies the Discrepancy Principle.
-
-    #         :rtype: float
-    #         :return: data misfit target
-
-    #     """
-    #     prob, survey = self.splitForward(forward)
-    #     return survey.nD
 
 
 class l2_DataMisfit(BaseDataMisfit):
@@ -103,10 +89,18 @@ class l2_DataMisfit(BaseDataMisfit):
         """
 
         if getattr(self, '_Wd', None) is None:
-            print 'SimPEG.l2_DataMisfit is creating default weightings for Wd.'
+            
             survey = self.survey
-            eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
-            self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+eps))
+
+            if getattr(survey,'std', None) is None:
+                print 'SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%'
+                survey.std = 0.05
+
+            if getattr(survey, 'eps', None) is None:
+                print 'SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||'
+                survey.eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
+
+            self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+survey.eps))
         return self._Wd
 
     @Wd.setter

SimPEG/Directives.py

@@ -239,51 +239,14 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
 
 
 
-class update_IRLS(InversionDirective):
+# class UpdateReferenceModel(Parameter):
 
-     m = None
-     eps_min = None
-     factor = None
-     gamma = None
-     phi_m_last = None
-     
-     def initialize(self):
-        
-         # Scale the regularization for changes in norm
-         if getattr(self, 'phi_m_last', None) is not None:
-             self.reg.gamma = 1.
-             phim_new = self.reg.eval(self.invProb.curModel)
-             self.gamma = self.phi_m_last / phim_new
-         
-         self.reg.gamma = self.gamma
-         
-     def endIter(self):
-         # Cool the threshold parameter 
-         if getattr(self, 'factor', None) is not None:
-             eps = self.reg.eps / self.factor
-             
-             if getattr(self, 'eps_min', None) is not None:
-                 self.reg.eps = np.max([self.eps_min,eps])
-             else:
-                 self.reg.eps = eps
-                 
-         
-         # Update the model used for the IRLS weights
-         if getattr(self, 'm', None) is None:
-             self.reg.m = self.invProb.curModel
-         
-         # 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.prob.mesh.nC))**2.
-         PC     = Utils.sdiag(diagA**-1.)
+#     mref0 = None
 
-         self.opt.approxHinv = PC
-         
-         phim_new = self.reg.eval(self.invProb.curModel)
-         self.reg.gamma = self.reg.gamma * self.invProb.phi_m_last / phim_new
-
-#==============================================================================
-#          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))
-#==============================================================================
+#     def nextIter(self):
+#         mref = getattr(self, 'm_prev', None)
+#         if mref is None:
+#             if self.debug: print 'UpdateReferenceModel is using mref0'
+#             mref = self.mref0
+#         self.m_prev = self.invProb.m_current
+#         return mref

SimPEG/EM/FDEM/FDEM.py

@@ -53,13 +53,13 @@ class BaseFDEMProblem(BaseEMProblem):
             Ainv.clean()
         return F
 
-    def Jvec(self, m, v, f=None):
+    def Jvec(self, m, v, u=None):
         """
             Sensitivity times a vector
         """
 
-        if f is None:
-           f = self.fields(m)
+        if u is None:
+           u = self.fields(m)
 
         self.curModel = m
 
@@ -71,34 +71,34 @@ class BaseFDEMProblem(BaseEMProblem):
 
             for src in self.survey.getSrcByFreq(freq):
                 ftype = self._fieldType + 'Solution'
-                u_src = f[src, ftype]
+                u_src = u[src, ftype]
                 dA_dm = self.getADeriv_m(freq, u_src, v)
                 dRHS_dm = self.getRHSDeriv_m(freq, src, v) 
                 du_dm = Ainv * ( - dA_dm + dRHS_dm )
                 
                 for rx in src.rxList:
-                    df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
+                    df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
                     df_dudu_dm = df_duFun(src, du_dm, adjoint=False)
 
-                    df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
+                    df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
                     df_dm = df_dmFun(src, v, adjoint=False)
 
                     Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex)
 
-                    P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m
+                    P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m
 
                     Jv[src, rx] = P(Df_Dm)
 
             Ainv.clean()
         return Utils.mkvc(Jv)
 
-    def Jtvec(self, m, v, f=None):
+    def Jtvec(self, m, v, u=None):
         """
             Sensitivity transpose times a vector
         """
 
-        if f is None:
-            f = self.fields(m)
+        if u is None:
+            u = self.fields(m)
 
         self.curModel = m
 
@@ -114,12 +114,12 @@ class BaseFDEMProblem(BaseEMProblem):
 
             for src in self.survey.getSrcByFreq(freq):
                 ftype = self._fieldType + 'Solution'
-                u_src = f[src, ftype]
+                u_src = u[src, ftype]
 
                 for rx in src.rxList:
-                    PTv = rx.projectFieldsDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
+                    PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
 
-                    df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None)
+                    df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
                     df_duT = df_duTFun(src, PTv, adjoint=True)
                     
                     ATinvdf_duT = ATinv * df_duT
@@ -128,7 +128,7 @@ class BaseFDEMProblem(BaseEMProblem):
                     dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True)
                     du_dmT = -dA_dmT + dRHS_dmT
 
-                    df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None)
+                    df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
                     dfT_dm = df_dmFun(src, PTv, adjoint=True)
 
                     du_dmT += dfT_dm
@@ -142,7 +142,8 @@ class BaseFDEMProblem(BaseEMProblem):
                         raise Exception('Must be real or imag')
             
             ATinv.clean()
-        return Jtv
+
+        return Utils.mkvc(Jtv)
 
     def getSourceTerm(self, freq):
         """

SimPEG/EM/FDEM/SrcFDEM.py

@@ -277,10 +277,10 @@ class CircularLoop(BaseSrc):
             if not prob.mesh.isSymmetric:
                 # TODO ?
                 raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', self.radius, mu=self.mu)
+            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
-            srcfct = MagneticLoopVectorPotential
+            srcfct = MagneticDipoleVectorPotential
             ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
             ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
             az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -26,13 +26,6 @@ class Rx(SimPEG.Survey.BaseRx):
                     'byi':['b', 'Fy', 'imag'],
                     'bzi':['b', 'Fz', 'imag'],
 
-                    'bxr_sec':['bSecondary', 'Fx', 'real'],
-                    'byr_sec':['bSecondary', 'Fy', 'real'],
-                    'bzr_sec':['bSecondary', 'Fz', 'real'],
-                    'bxi_sec':['bSecondary', 'Fx', 'imag'],
-                    'byi_sec':['bSecondary', 'Fy', 'imag'],
-                    'bzi_sec':['bSecondary', 'Fz', 'imag'],
-
                     'jxr':['j', 'Fx', 'real'],
                     'jyr':['j', 'Fy', 'real'],
                     'jzr':['j', 'Fz', 'real'],
@@ -113,7 +106,7 @@ class Survey(SimPEG.Survey.BaseSurvey):
         SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
 
         _freqDict = {}
-        for src in self.srcList:
+        for src in srcList:
             if src.freq not in _freqDict:
                 _freqDict[src.freq] = []
             _freqDict[src.freq] += [src]

SimPEG/EM/TDEM/SurveyTDEM.py

@@ -79,32 +79,12 @@ class SrcTDEM(Survey.BaseSrc):
 
 class SrcTDEM_VMD_MVP(SrcTDEM):
 
-    def __init__(self,rxList,loc,waveformType="STEPOFF"):
+    def __init__(self,rxList,loc):
         self.loc = loc
-        self.waveformType = waveformType
         SrcTDEM.__init__(self,rxList)
 
     def getInitialFields(self, mesh):
         """Vertical magnetic dipole, magnetic vector potential"""
-        if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"        
-            if mesh._meshType is 'CYL':
-                if mesh.isSymmetric:
-                    MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
-                else:
-                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            elif mesh._meshType is 'TENSOR':
-                MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
-            else:
-                raise Exception('Unknown mesh for VMD') 
-            return {"b": mesh.edgeCurl*MVP}               
-        elif self.waveformType == "GENERAL":
-            print ">> General waveform: Zero initial condition"
-            return {"b": np.zeros(mesh.nF)}
-        else:
-            raise NotImplementedError("Only use STEPOFF or GENERAL")
-
-    def getMeS(self, mesh, MfMui):
         if mesh._meshType is 'CYL':
             if mesh.isSymmetric:
                 MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
@@ -113,12 +93,13 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
         elif mesh._meshType is 'TENSOR':
             MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
         else:
-            raise Exception('Unknown mesh for VMD')            
-        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
+            raise Exception('Unknown mesh for VMD')
+
+        return {"b": mesh.edgeCurl*MVP}
 
 
 class SrcTDEM_CircularLoop_MVP(SrcTDEM):
-    def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
+    def __init__(self,rxList,loc,radius,waveformType):
         self.loc = loc
         self.radius = radius
         self.waveformType = waveformType

SimPEG/EM/__init__.py

@@ -1,6 +1,6 @@
-# from EM import *
 import TDEM
 import FDEM
 import Base
 import Analytics
 import Utils
+from scipy.constants import mu_0, epsilon_0

SimPEG/Examples/DC_PseudoSection_Simulation.py

@@ -1,179 +0,0 @@
-from SimPEG import *
-import simpegDCIP as DC
-import scipy.interpolate as interpolation
-import matplotlib.pyplot as plt
-import time
-import re
-
-def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi=np.r_[25.,25.], param = np.r_[30.,30.,5], stype = 'dpdp',  plotIt=True):
-    """
-    DC Forward Simulation
-    
-    Forward model conductive spheres in a half-space and plot a pseudo-section
-    
-    Created on Mon Feb 01 19:28:06 2016
-    
-    @fourndo
-    """
-    
-    # First we need to create a mesh and a model.
-    
-    # This is our mesh
-    dx    = 5.
-    
-    hxind = [(dx,15,-1.3), (dx, 75), (dx,15,1.3)]
-    hyind = [(dx,15,-1.3), (dx, 10), (dx,15,1.3)]
-    hzind = [(dx,15,-1.3),(dx, 15)]
-    
-    mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
-    
-    
-    # Set background conductivity
-    model = np.ones(mesh.nC) * sig[0]
-    
-    # First anomaly
-    ind = Utils.ModelBuilder.getIndicesSphere(loc[:,0],radi[0],mesh.gridCC)
-    model[ind] = sig[1]
-    
-    # Second anomaly
-    ind = Utils.ModelBuilder.getIndicesSphere(loc[:,1],radi[1],mesh.gridCC)
-    model[ind] = sig[2]
-    
-    # Get index of the center
-    indy = int(mesh.nCy/2)
-
-    
-    # Plot the model for reference
-    # Define core mesh extent
-    xlim = 200
-    zlim = 125
-    
-    # Specify the survey type: "pdp" | "dpdp"
-    
-    
-    # Then specify the end points of the survey. Let's keep it simple for now and survey above the anomalies, top of the mesh
-    ends = [(-175,0),(175,0)]
-    ends = np.c_[np.asarray(ends),np.ones(2).T*mesh.vectorNz[-1]]
-    
-    # Snap the endpoints to the grid. Easier to create 2D section.
-    indx = Utils.closestPoints(mesh, ends )
-    locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]]
-    
-    # We will handle the geometry of the survey for you and create all the combination of tx-rx along line
-    [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2])
-    
-    # Define some global geometry
-    dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) ) 
-    dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len
-    dl_y = ( Tx[-1][1,1] - Tx[0][1,0]  ) / dl_len
-    azm =  np.arctan(dl_y/dl_x)
-
-    #Set boundary conditions
-    mesh.setCellGradBC('neumann')
-    
-    # Define the differential operators needed for the DC problem
-    Div = mesh.faceDiv
-    Grad = mesh.cellGrad
-    Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
-    
-    A = Div*Msig*Grad
-    
-    # Change one corner to deal with nullspace
-    A[0,0] = 1
-    A = sp.csc_matrix(A)
-    
-    # We will solve the system iteratively, so a pre-conditioner is helpful
-    # This is simply a Jacobi preconditioner (inverse of the main diagonal)
-    dA = A.diagonal()
-    P = sp.spdiags(1/dA,0,A.shape[0],A.shape[0])
-
-    # Now we can solve the system for all the transmitters
-    # We want to store the data
-    data = []
-    
-    # There is probably a more elegant way to do this, but we can just for-loop through the transmitters
-    for ii in range(len(Tx)):
-        
-        start_time = time.time() # Let's time the calculations
-        
-        #print("Transmitter %i / %i\r" % (ii+1,len(Tx)))
-        
-        # Select dipole locations for receiver
-        rxloc_M = np.asarray(Rx[ii][:,0:3])
-        rxloc_N = np.asarray(Rx[ii][:,3:])
-
-      
-        # For usual cases "dpdp" or "gradient"
-        if not re.match(stype,'pdp'): 
-            inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T )
-            RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] )   
-            
-        else: 
-            
-            # Create an "inifinity" pole
-            tx =  np.squeeze(Tx[ii][:,0:1])
-            tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
-            inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T)
-            RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1] / mesh.vol[inds] )  
-    
-    
-        # Iterative Solve
-        Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
-    
-        # We now have the potential everywhere
-        phi = mkvc(Ainvb[0])
-        
-        # Solve for phi on pole locations
-        P1 = mesh.getInterpolationMat(rxloc_M, 'CC')
-        P2 = mesh.getInterpolationMat(rxloc_N, 'CC')
-            
-        # Compute the potential difference
-        dtemp = (P1*phi - P2*phi)*np.pi
-        
-        data.append( dtemp )     
-        print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time),
-        
-    print 'Transmitter {0} of {1}'.format(ii,len(Tx))
-    print 'Forward completed'
-    
-    
-    # Let's just convert the 3D format into 2D (distance along line) and plot
-    [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(Tx,Rx)
-    
-
-    # Here is an example for the first tx-rx array
-    if plotIt:
-        fig = plt.figure()
-        ax = plt.subplot(2,1,1, aspect='equal')
-        mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True)
-        ax.set_title('E-W section at '+str(mesh.vectorCCy[indy])+' m')
-        plt.gca().set_aspect('equal', adjustable='box')
-
-        plt.scatter(Tx[0][0,:],Tx[0][2,:],s=40,c='g', marker='v')
-        plt.scatter(Rx[0][:,0::3],Rx[0][:,2::3],s=40,c='y')
-        plt.xlim([-xlim,xlim])
-        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
-        
-        
-        ax = plt.subplot(2,1,2, aspect='equal')   
-        
-        # Plot the location of the spheres for reference
-        circle1=plt.Circle((loc[0,0]-Tx[0][0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
-        circle2=plt.Circle((loc[0,1]-Tx[0][0,0],loc[2,1]),radi[1],color='k',fill=False, lw=3)
-        ax.add_artist(circle1)
-        ax.add_artist(circle2)  
-        
-        # Add the speudo section
-        DC.plot_pseudoSection(Tx2d,Rx2d,data,mesh.vectorNz[-1],stype)
-        
-        plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
-        plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y')
-        plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k')
-        plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
-        
-        plt.show()
-        
-        return fig, ax
-        
-if __name__ == '__main__':
-    run()

SimPEG/Examples/EM_FDEM_1D_Inversion.py

@@ -0,0 +1,116 @@
+from SimPEG import *
+import SimPEG.EM as EM
+from SimPEG.EM import mu_0
+
+
+def run(plotIt=True):
+    """
+        EM: FDEM: 1D: Inversion
+        =======================
+
+        Here we will create and run a FDEM 1D inversion.
+
+    """
+
+    cs, ncx, ncz, npad = 5., 25, 15, 15
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hz], '00C')
+
+    layerz = -100.
+
+    active = mesh.vectorCCz<0.
+    layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=layerz)
+    actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
+    sig_half = 2e-2
+    sig_air = 1e-8
+    sig_layer = 1e-2
+    sigma = np.ones(mesh.nCz)*sig_air
+    sigma[active] = sig_half
+    sigma[layer] = sig_layer
+    mtrue = np.log(sigma[active])
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        fig, ax = plt.subplots(1,1, figsize = (3, 6))
+        plt.semilogx(sigma[active], mesh.vectorCCz[active])
+        ax.set_ylim(-500, 0)
+        ax.set_xlim(1e-3, 1e-1)
+        ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
+        ax.set_ylabel('Depth (m)', fontsize = 14)
+        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
+
+
+    rxOffset=10. 
+    bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi')
+
+    freqs = np.logspace(1,3,10)
+    srcLoc = np.array([0., 0., 10.])
+
+    srcList = []
+    [srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
+
+    survey = EM.FDEM.Survey(srcList)
+    prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
+
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+        prb.Solver = SolverLU
+
+    prb.pair(survey)
+
+    std = 0.05
+    survey.makeSyntheticData(mtrue, std)
+
+    survey.std = std
+    survey.eps = np.linalg.norm(survey.dtrue)*1e-5
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        fig, ax = plt.subplots(1,1, figsize = (6, 6))
+        ax.semilogx(freqs,survey.dtrue[:freqs.size], 'b.-')
+        ax.semilogx(freqs,survey.dobs[:freqs.size], 'r.-')
+        ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16)
+        ax.set_xlabel('Time (s)', fontsize = 14)
+        ax.set_ylabel('$B_z$ (T)', fontsize = 16)
+        ax.set_xlabel('Time (s)', fontsize = 14)
+        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
+
+    dmisfit = DataMisfit.l2_DataMisfit(survey)
+    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
+    reg = Regularization.Tikhonov(regMesh)
+    opt = Optimization.InexactGaussNewton(maxIter = 6)
+    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
+
+    # Create an inversion object
+    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
+    betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaest])
+    m0 = np.log(np.ones(mtrue.size)*sig_half)
+    reg.alpha_s = 1e-3
+    reg.alpha_x = 1.
+    prb.counter = opt.counter = Utils.Counter()
+    opt.LSshorten = 0.5
+    opt.remember('xc')
+
+    mopt = inv.run(m0)
+
+    if plotIt:
+        import matplotlib.pyplot as plt
+        fig, ax = plt.subplots(1,1, figsize = (3, 6))
+        plt.semilogx(sigma[active], mesh.vectorCCz[active])
+        plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])
+        ax.set_ylim(-500, 0)
+        ax.set_xlim(1e-3, 1e-1)
+        ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
+        ax.set_ylabel('Depth (m)', fontsize = 14)
+        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
+        plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'],loc='best')
+        plt.show()
+
+
+if __name__ == '__main__':
+    run()

SimPEG/Examples/EM_TDEM_1D_Inversion.py

@@ -1,6 +1,6 @@
 from SimPEG import *
 import SimPEG.EM as EM
-from scipy.constants import mu_0
+from SimPEG.EM import mu_0
 
 
 def run(plotIt=True):
@@ -50,20 +50,18 @@ def run(plotIt=True):
     prb.Solver = SolverLU
     prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)]
     prb.pair(survey)
-    dtrue = survey.dpred(mtrue)
 
-
-    survey.dtrue = dtrue
+    # create observed data
     std = 0.05
-    noise = std*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
-    survey.dobs = survey.dtrue+noise
-    survey.std = survey.dobs*0 + std
-    survey.Wd = 1/(abs(survey.dobs)*std)
+    
+    survey.dobs = survey.makeSyntheticData(mtrue,std)
+    survey.std = std 
+    survey.eps = 1e-5*np.linalg.norm(survey.dobs)
 
     if plotIt:
         import matplotlib.pyplot as plt
         fig, ax = plt.subplots(1,1, figsize = (10, 6))
-        ax.loglog(rx.times, dtrue, 'b.-')
+        ax.loglog(rx.times, survey.dtrue, 'b.-')
         ax.loglog(rx.times, survey.dobs, 'r.-')
         ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16)
         ax.set_xlabel('Time (s)', fontsize = 14)
@@ -76,6 +74,7 @@ def run(plotIt=True):
     reg = Regularization.Tikhonov(regMesh)
     opt = Optimization.InexactGaussNewton(maxIter = 5)
     invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
+
     # Create an inversion object
     beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
     betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)

SimPEG/Examples/Mesh_OcTree_Creation.py

@@ -0,0 +1,125 @@
+
+if __name__ == '__main__':
+
+
+    import matplotlib.pyplot as plt
+    import matplotlib
+    from mpl_toolkits.mplot3d import Axes3D
+    import matplotlib.colors as colors
+    import matplotlib.cm as cmx
+
+    def topo(x):
+        return np.sin(x*(2.*np.pi))*0.3 + 0.5
+
+    def function(cell):
+        r = cell.center - np.array([0.5]*len(cell.center))
+        dist = np.sqrt(r.dot(r))
+        # dist2 = np.abs(cell.center[-1] - topo(cell.center[0]))
+
+        # dist = min([dist1,dist2])
+        # if dist < 0.05:
+        #     return 5
+        if dist < 0.1:
+            return 5
+        if dist < 0.2:
+            return 4
+        if dist < 0.4:
+            return 3
+        return 2
+
+    # T = TreeMesh([[(1,128)],[(1,128)],[(1,128)]],levels=7)
+    # T = TreeMesh([128,128,128])
+    # T = TreeMesh([64,64],levels=6)
+    T = TreeMesh([8,8])
+    # T = TreeMesh([[(1,128)],[(1,128)]],levels=7)
+    # T.refine(lambda xc:2, balance=False)
+    # T._index([0,0,0])
+    # T._pointer(0)
+
+
+    # tic = time.time()
+    T.refine(function)#, balance=False)
+    # print time.time() - tic
+    # print T.nC
+    # T.plotSlice(np.log(T.vol))#np.random.rand(T.nC))
+    T.plotGrid()
+    # print [c for c in T]
+    c = T[0]
+    plt.plot(c.center[0],c.center[1],'r.')
+    nodes = c.nodes
+    for n in nodes:
+        _ = T._gridN[n,:]
+        plt.plot(_[0],_[1],'gs')
+    plt.show()
+    blah
+
+    # T.plotImage(np.arange(len(T.vol)),showIt=True)
+
+    # print T.getFaceInnerProduct()
+    # print T.gridFz
+
+
+    # T._refineCell([8,0,1])
+    # T._refineCell([8,0,2])
+    # T._refineCell([12,0,2])
+    # T._refineCell([8,4,2])
+    # T._refineCell([6,0,3])
+    # T._refineCell([8,8,1])
+    # T._refineCell([0,0,0,1])
+    # T.__dirty__ = True
+
+
+    # print T.gridFx.shape[0], T.nFx
+
+
+
+    ax = plt.subplot(211)
+    ax.spy(T.edgeCurl)
+
+    # print Mesh.TensorMesh([2,2,2]).edgeCurl.todense()
+    # print T.edgeCurl.todense()
+    # print Mesh.TensorMesh([2,2,2]).edgeCurl.todense() - T.edgeCurl.todense()
+    # print T.gridEy - Mesh.TensorMesh([2,2,2]).gridEy
+
+    # print T.edge
+    # T.plotGrid(ax=ax)
+
+    # R = deflationMatrix(T._facesX, T._hangingFx, T._fx2i)
+    # print R
+
+    ax = plt.subplot(212)#, projection='3d')
+    ax.spy(Mesh.TensorMesh([2,2,2]).edgeCurl)
+
+    # ax = plt.subplot(313)
+    # ax.spy(T.faceDiv[:,:T.nFx] * R)
+
+
+    # T.balance()
+    # T.plotGrid(ax=ax)
+
+    # cx = T._getNextCell([0,0,1],direction=0,positive=True)
+    # print cx
+    # # print [T._asPointer(_) for _ in cx]
+    # cx = T._getNextCell([8,0,3],direction=0,positive=False)
+    # print T._asPointer(cx)
+    # cx = T._getNextCell([8,8,1],direction=1,positive=False)
+    # print cx, #[T._asPointer(_) for _ in cx]
+    # cm = T._getNextCell([64,80,4],direction=0,positive=False)
+    # cy = T._getNextCell([64,80,4],direction=1,positive=True)
+    # cp = T._getNextCell([64,80,4],direction=1,positive=False)
+
+    # ax.plot( T._cellN([4,0,1])[0],T._cellN([4,0,1])[1], 'yd')
+    # ax.plot( T._cellN(cx)[0],T._cellN(cx)[1], 'ys')
+    # ax.plot( T._cellN(cm)[0],T._cellN(cm)[1], 'ys')
+    # ax.plot( T._cellN(cy)[0],T._cellN(cy)[1], 'ys')
+    # ax.plot( T._cellN(cp[0])[0],T._cellN(cp[0])[1], 'ys')
+    # ax.plot( T._cellN(cp[1])[0],T._cellN(cp[1])[1], 'ys')
+
+
+
+
+
+    # print T.nN
+
+    plt.show()
+

SimPEG/Examples/__init__.py

@@ -1,6 +1,7 @@
 # Run this file to add imports.
 
 ##### AUTOIMPORTS #####
+import EM_FDEM_1D_Inversion
 import EM_FDEM_Analytic_MagDipoleWholespace
 import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
@@ -14,7 +15,7 @@ import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
 
-__examples__ = ["EM_FDEM_Analytic_MagDipoleWholespace", "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"]
+__examples__ = ["EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "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"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Exceptions.py

@@ -1,13 +0,0 @@
-
-
-class SimPEGException(Exception):
-
-    def __init__(self, reason=''):
-        self.reason = reason
-
-    def __str__(self):
-        return '%s: %s' %(self.__class__.__name__,  self.reason)
-
-
-class PairingException(SimPEGException):
-    pass

SimPEG/Maps.py

@@ -4,7 +4,6 @@ from Tests import checkDerivative
 from PropMaps import PropMap, Property
 from numpy.polynomial import polynomial
 from scipy.interpolate import UnivariateSpline
-from scipy.spatial import cKDTree
 
 class IdentityMap(object):
     """
@@ -122,7 +121,7 @@ class IdentityMap(object):
             if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]:
                 raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape)))
             return self._transform(val)
-        raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray! Not a %s'%type(val))
+        raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!')
 
     def __str__(self):
         return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
@@ -459,89 +458,6 @@ class Mesh2Mesh(IdentityMap):
         return self.P
 
 
-class Mesh2MeshTopo(IdentityMap):
-    """
-        Takes a model on one mesh are translates it to another mesh
-        with consideration of topography
-
-    """
-    tree = None
-    nIterpPts = 6
-    P = None #: The CSR projection matrix.
-
-    def __init__(self, meshes, actinds, **kwargs):
-        Utils.setKwargs(self, **kwargs)
-
-        assert type(meshes) is list, "meshes must be a list of two meshes"
-        assert len(meshes) == 2, "meshes must be a list of two meshes"
-        assert type(actinds) is list, "actinds must be a list of two meshes"
-        assert len(actinds) == 2, "actinds must be a list of two meshes"
-        assert meshes[0].dim == meshes[1].dim, """The two meshes must be the same dimension"""
-
-        self.mesh  = meshes[0]
-        self.mesh2 = meshes[1]
-        self.actind = actinds[0]
-        self.actind2 = actinds[1]
-        self._createProjection()
-
-        # Old version using SimPEG interpolation
-        # self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC,'CC',zerosOutside=True)
-
-    def genActiveindfromTopo(mesh, xyztopo):
-        #TODO: This possibly needs to be improved use vtk(?)
-        if mesh.dim==3:
-            nCxy = mesh.nCx*mesh.nCy
-            Zcc = mesh.gridCC[:,2].reshape((nCxy, mesh.nCz), order='F')
-            Ftopo = NearestNDInterpolator(xyztopo[:,:2], xyztopo[:,2])
-            XY = Utils.ndgrid(mesh.vectorCCx, mesh.vectorCCy)
-            XY.shape
-            topo = Ftopo(XY)
-            actind = []
-            for ixy in range(nCxy):
-                actind.append(topo[ixy] <= Zcc[ixy,:])
-        else:
-            raise NotImplementedError("Only 3D is working")
-
-        return Utils.mkvc(np.vstack(actind))
-
-    #Question .. is it only generated once?
-    def _createProjection(self):
-        """
-            KD Tree interpolation onto the active cells.
-        """
-        if self.tree==None:
-            self.tree = cKDTree(zip(self.mesh.gridCC[self.actind,0], self.mesh.gridCC[self.actind,1], self.mesh.gridCC[self.actind,2]))
-        d, inds = self.tree.query(zip(self.mesh2.gridCC[self.actind2,0],self.mesh2.gridCC[self.actind2,1],self.mesh2.gridCC[self.actind2,2]), k=self.nIterpPts)
-        # Not sure consideration of the volume ...
-        # vol = np.zeros((self.actind2.sum(), self.nIterpPts))
-        # for i in range(self.nIterpPts):
-        #     vol[:,i] = self.mesh.vol[inds[:,i]]
-        w = 1. / d**2
-        w = Utils.sdiag(1./np.sum(w, axis=1)) * (w)
-        I = Utils.mkvc(np.arange(inds.shape[0]).reshape([-1,1]).repeat(self.nIterpPts, axis=1))
-        J = Utils.mkvc(inds)
-        P = sp.coo_matrix( (Utils.mkvc(w),(I, J)), shape=(inds.shape[0], (self.actind).sum()) )
-        # self.P = Utils.sdiag(self.mesh2.vol[self.actind2])*P.tocsc()
-        self.P = P.tocsr()
-
-    @property
-    def shape(self):
-        """Number of parameters in the model."""
-        # return (self.mesh.nC, self.mesh2.nC)
-        return (self.actind2.sum(), self.actind.sum())
-
-    @property
-    def nP(self):
-        """Number of parameters in the model."""
-        # return self.mesh2.nC
-        return self.actind2.sum()
-
-    def _transform(self, m):
-        return self.P*m
-
-    def deriv(self, m):
-        return self.P
-
 class ActiveCells(IdentityMap):
     """
         Active model parameters.

SimPEG/Mesh/DiffOperators.py

@@ -565,58 +565,7 @@ class DiffOperators(object):
 
         return Pbc, Pin, Pout
 
-    
-    def unitCellGradx():
-        doc = """Cell centered Gradient in the x dimension used for 
-                regularization. The gradient operator is square (nC-by-nC)"""
-        def fget(self):
-            if self.dim < 3: return None
-            if getattr(self, '_unitCellGradx', None) is None:
-                
-                n = self.vnC
-                gx = ddx(n[0]-1)
-                gx_square = sp.vstack((gx,gx[-1,:]*-1), format="csr")
-                
-                self._unitCellGradx = kron3(speye(n[2]), speye(n[1]), gx_square)
-        
-            return self._unitCellGradx
-        return locals()
-    unitCellGradx = property(**unitCellGradx())
-    
-    def unitCellGrady():
-        doc = """Cell centered Gradient in they dimension used for 
-                regularization. The gradient operator is square (nC-by-nC)"""
-        def fget(self):
-            if self.dim < 3: return None
-            if getattr(self, '_unitCellGrady', None) is None:
-        
-                n = self.vnC
-                gy = ddx(n[1]-1)
-                gy_square = sp.vstack((gy,gy[-1,:]*-1), format="csr")
-                
-                self._unitCellGrady = kron3(speye(n[2]), gy_square, speye(n[0]))
-            
-            return self._unitCellGrady
-        return locals()
-    unitCellGrady = property(**unitCellGrady())
-        
-    def unitCellGradz():
-        doc = """Cell centered Gradient in they dimension used for 
-                regularization. The gradient operator is square (nC-by-nC)"""
-        def fget(self):
-            if self.dim < 3: return None
-            if getattr(self, '_unitCellGradz', None) is None:
-        
-                n = self.vnC
-                gz = ddx(n[2]-1)
-                gz_square = sp.vstack((gz,gz[-1,:]*-1), format="csr")
-                
-                self._unitCellGradz = kron3( gz_square , speye(n[1]), speye(n[0]))
-            
-            return self._unitCellGradz
-        return locals()
-    unitCellGradz = property(**unitCellGradz())        
-        
+
     # --------------- Averaging ---------------------
 
     @property

SimPEG/Optimization.py

@@ -989,19 +989,5 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
                 if np.logical_or(norm(resid)/normResid0 <= self.tolCG, cgiter == self.maxIterCG):
                     cgFlag = 1
                 # End CG Iterations
-                    
-            # Take a gradient step on the active cells if exist
-            if temp != self.xc.size:
-        
-                rhs_a = (Active) * -self.g
-                
-                dm_i = max( abs( delx ) )
-                dm_a = max( abs(rhs_a) )           
-    
-                delx = delx + rhs_a * dm_i / dm_a /10.
-                    
-            # 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.
 
-            return delx
+            return delx

SimPEG/Problem.py

@@ -1,6 +1,6 @@
 import Utils, Survey, Models, numpy as np, scipy.sparse as sp
 Solver = Utils.SolverUtils.Solver
-import Maps, Mesh, Exceptions
+import Maps, Mesh
 from Fields import Fields, TimeFields
 
 class BaseProblem(object):
@@ -18,13 +18,9 @@ class BaseProblem(object):
     Solver = Solver   #: A SimPEG Solver class.
     solverOpts = {}   #: Sovler options as a kwarg dict
 
-    PropMap = None    #: A SimPEG PropertyMap class.
+    mesh    = None    #: A SimPEG.Mesh instance.
 
-    def __init__(self, mesh, mapping=None, **kwargs):
-        Utils.setKwargs(self, **kwargs)
-        assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
-        self.mesh = mesh
-        self.mapping = mapping or Maps.IdentityMap(mesh)
+    PropMap = None    #: A SimPEG PropertyMap class.
 
     @property
     def mapping(self):
@@ -36,9 +32,14 @@ class BaseProblem(object):
             val._assertMatchesPair(self.mapPair)
             self._mapping = val
         else:
-            self._propMapMapping = val
             self._mapping = self.PropMap(val)
 
+    def __init__(self, mesh, mapping=None, **kwargs):
+        Utils.setKwargs(self, **kwargs)
+        assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
+        self.mesh = mesh
+        self.mapping = mapping or Maps.IdentityMap(mesh)
+
     @property
     def survey(self):
         """
@@ -46,22 +47,13 @@ class BaseProblem(object):
         """
         return getattr(self, '_survey', None)
 
-    def pair(self, survey):
+    def pair(self, d):
         """Bind a survey to this problem instance using pointers."""
-        assert isinstance(survey, self.surveyPair), "Survey must be an instance of a %s class."%(self.surveyPair.__name__)
-        if survey.ispaired:
+        assert isinstance(d, self.surveyPair), "Data object must be an instance of a %s class."%(self.surveyPair.__name__)
+        if d.ispaired:
             raise Exception("The survey object is already paired to a problem. Use survey.unpair()")
-        try:
-            self._survey = survey
-            self._validatePairing()
-        except Exceptions.PairingException, e:
-            self._survey = None
-            raise e
-        survey._prob = self
-
-    def _validatePairing(self):
-        """Called when the pair is done, raise a SimPEG.Exceptions.PairingException if unsuccessful"""
-        pass
+        self._survey = d
+        d._prob = self
 
     def unpair(self):
         """Unbind a survey from this problem instance."""
@@ -166,9 +158,6 @@ class BaseProblem(object):
 
 class BaseTimeProblem(BaseProblem):
     """Sets up that basic needs of a time domain problem."""
-    
-    waveformType = "STEPOFF"    
-    current = None
 
     @property
     def timeSteps(self):
@@ -195,11 +184,6 @@ class BaseTimeProblem(BaseProblem):
         self._timeSteps = Utils.meshTensor(value)
         del self.timeMesh
 
-    def currentwaveform(self, wave):
-        self._timeSteps = np.diff(wave[:,0])
-        self.current = wave[:,1]
-        self.waveformType = "GENERAL"
-
     @property
     def nT(self):
         "Number of time steps."
@@ -230,191 +214,4 @@ class BaseTimeProblem(BaseProblem):
             del self._timeMesh
 
 
-class GlobalProblem(BaseProblem):
-    """
-
-        The GlobalProblem allows you to run a whole bunch of SubProblems,
-        potentially in parallel, potentially of different meshes.
-
-        This is handy for working with lots of sources,
-
-    """
-
-    surveyKwargs = {}
-    probKwargs   = {}
-
-    def __init__(self, SubProblem, globalMesh, mapping=None, **kwargs):
-
-        # assert isclass??(SubProblem, BaseProblem), "SubProblem must be a SimPEG.Problem.BaseProblem object."
-        self.surveyPair = SubProblem.surveyPair
-        self.PropMap = SubProblem.PropMap
-        self.mapPair = SubProblem.mapPair
-        self.SubProblem = SubProblem
-
-        Utils.setKwargs(self, **kwargs)
-        assert isinstance(globalMesh, Mesh.BaseMesh), "globalMesh must be a SimPEG.Mesh object."
-        self.globalMesh = globalMesh
-        self.mapping = mapping or Maps.IdentityMap()
-
-    @property
-    def groups(self):
-        """
-            List of lists/integers to say how the sources are grouped.
-
-            e.g.
-
-            survey.srcList = [s0,s1,s2,s3,s4]
-            groups = [ [0,4], [1,3], 2 ]
-        """
-        if getattr(self, '_groups', None) is None:
-            if not self.ispaired: return None
-            self._groups = range(self.survey.nSrc)
-        return self._groups
-    @groups.setter
-    def groups(self, val):
-        assert type(val) is list, 'This should be an list of groups'
-        if self.ispaired:
-            for g in val:
-                assert type(g) in [int, list], 'Must be an integer or a list'
-                if type(g) is int:
-                    assert g >= 0 and g < self.survey.nSrc, '%d is outside the number of sources in the surveys list'%g
-                if type(g) is list:
-                    for sg in g:
-                        assert type(g) is int, 'Must be an integer or a list'
-                        assert g >= 0 and g < self.survey.nSrc, '%d is outside the number of sources in the surveys list'%g
-            assert len(val) == len(self.survey.srcList), 'The groups must be the same length as the srcList in the survey'
-        self._groups = val
-        self._nGroups = None
-
-    @property
-    def meshes(self):
-        if getattr(self, '_meshes', None) is None:
-            if not self.ispaired: return None
-            self._meshes = [self.globalMesh]*self.nGroups
-        return self._meshes
-    @meshes.setter
-    def meshes(self, val):
-        assert type(val) is list
-        if self.ispaired:
-            assert len(val) == self.nGroups
-        self._meshes = val
-
-    @property
-    def nGroups(self):
-        if getattr(self, '_groups', None) is None:
-            return None
-        return len(self.groups)
-
-    def _validatePairing(self):
-        try:
-            self.groups = self.groups # check the assumptions for the grouping
-        except Exception, e:
-            raise Exceptions.PairingException(reason='The grouping does not match the survey')
-        if self.nGroups is not len(self.meshes):
-            raise Exceptions.PairingException(reason='The meshes are not the the same length as the number of groups')
-
-    def getSubProblemandSubSurvey(self, subMap, ind):
-        #This is a core place that we can proceed parallelization
-        assert self.ispaired, 'You must be paired to a survey'
-        assert type(ind) in [int,long] and ind >= 0 and ind < self.nGroups, 'ind must be an index into the group list'
-
-        subMesh = self.meshes[ind]
-        # subMap  = Maps.IdentityMap(subMesh) # this is probably a mesh2mesh mapping?
-        # subMap  = self.getSubMap(subMesh, ind)
-
-        if self.PropMap is None:
-            prob = self.SubProblem(subMesh, mapping=subMap * self.mapping, **self.probKwargs)
-        else:
-            # This will not work with a fancier propmap...
-            prob = self.SubProblem(subMesh, mapping=subMap * self._propMapMapping, **self.probKwargs)
-
-        survey = self.survey.__class__(srcList=self.survey.srcList[self.groups[ind]], **self.surveyKwargs)
-        prob.pair(survey)
-
-        return prob, survey
-
-    # Not sure we need this here ... 
-
-    def getSubMap(self, subMesh, ind):
-        """The sub"""
-        mesh2mesh = Maps.IdentityMap(subMesh) # this is probably a mesh2mesh mapping?
-
-        if self.PropMap is None:
-            subMap = mesh2mesh * self.mapping
-        else:
-            subMap = mesh2mesh * self._propMapMapping
-
-        return subMap
-
-if __name__ == '__main__':
-
-
-    from SimPEG import *
-    from SimPEG import EM
-    from scipy.constants import mu_0
-    from pymatsolver import MumpsSolver
-
-    cs = 10.
-    ncx, ncy, ncz = 10, 10, 10
-    npad = 4
-    freq = 1e2
-
-    hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
-    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
-    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
-    mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC')
-
-    mapping = Maps.ExpMap(mesh)
-
-    x = np.linspace(-10,10,5)
-    XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
-    rxList = EM.FDEM.Rx(XYZ, 'exi')
-    Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
-    Src1 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
-
-
-    prb0 = EM.FDEM.Problem_b(mesh, mapping=mapping, Solver=MumpsSolver)
-    survey = EM.FDEM.Survey([Src0])
-    prb0.pair(survey)
-    prb1 = EM.FDEM.Problem_b(mesh, mapping=mapping, Solver=MumpsSolver)
-    survey = EM.FDEM.Survey([Src1])
-    prb1.pair(survey)
-
-
-
-    sig = 1e-1
-    sigma = np.ones(mesh.nC)*sig
-    sigma[mesh.gridCC[:,2] > 0] = 1e-8
-    m = np.log(sigma)
-
-    GP = GlobalProblem(EM.FDEM.Problem_b, mesh, mapping=mapping, meshes=[mesh,mesh])
-    survey = EM.FDEM.Survey([Src0, Src1])
-    GP.pair(survey)
-
-    gp1 = GP.getSubProblem(0)
-    gp1.Solver = MumpsSolver
-
-    pu = prb0.fields(m)
-    gpu = gp1.fields(m)
-
-    bfz = mesh.r(pu[Src0, 'b'],'F','Fz','M')
-    bfz = mesh.r(gpu[Src0, 'b'],'F','Fz','M')
-    x = np.linspace(-55,55,12)
-    XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
-    P = mesh.getInterpolationMat(XYZ, 'Fz')
-
-    # an = EM.Analytics.FDEM.hzAnalyticDipoleF(x, Src0.freq, sig)
-
-    # diff = np.log10(np.abs(P*np.imag(pu[Src0, 'b']) - mu_0*np.imag(an)))
-    # diff = np.log10(np.abs(P*np.imag(gpu[Src0, 'b']) - mu_0*np.imag(an)))
-
-    import matplotlib.pyplot as plt
-    plt.plot(x,np.log10(np.abs(P*np.imag(pu[Src0, 'b']))), 'r-s')
-    plt.plot(x,np.log10(np.abs(P*np.imag(gpu[Src0, 'b']))), 'b')
-    # plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r')
-    # plt.plot(x,diff,'g')
-    plt.show()
-
-
-
 

SimPEG/Regularization.py

@@ -282,239 +282,3 @@ class Tikhonov(BaseRegularization):
             out = mD.T * ( self.W.T * r )
         return out
 
-class Simple(BaseRegularization):
-    """
-        Only for tensor mesh
-    """     
-    
-    smoothModel = 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")
-    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, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
-
-    
-        
-    @property
-    def Ws(self):
-        """Regularization matrix Ws"""
-        if getattr(self,'_Ws', None) is None:
-                self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
-        return self._Ws
-
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
-        if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.mesh.unitCellGradx
-        return self._Wx
-
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
-        if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.mesh.unitCellGrady
-        return self._Wy
-
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
-        if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz
-        return self._Wz
-
-    @property
-    def Wxx(self):
-        """Regularization matrix Wxx"""
-        if getattr(self, '_Wxx', None) is None:
-            self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
-        return self._Wxx
-
-    @property
-    def Wyy(self):
-        """Regularization matrix Wyy"""
-        if getattr(self, '_Wyy', None) is None:
-            self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
-        return self._Wyy
-
-    @property
-    def Wzz(self):
-        """Regularization matrix Wzz"""
-        if getattr(self, '_Wzz', None) is None:
-            self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
-        return self._Wzz
-
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx, self.Wxx)
-            if self.mesh.dim > 1:
-                wlist += (self.Wy, self.Wyy)
-            if self.mesh.dim > 2:
-                wlist += (self.Wz, self.Wzz)
-            self._Wsmooth = sp.vstack(wlist)
-        return self._Wsmooth
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Ws, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
-    @Utils.timeIt
-    def eval(self, m):
-        if self.smoothModel == True:
-            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.smoothModel == False:
-            r = self.W * ( self.mapping * (m - self.mref) )
-            return 0.5*r.dot(r)
-
-
-    @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.smoothModel == 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.smoothModel == 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 SparseRegularization(Simple):
-    
-    eps = 1e-1
-
-    m = None
-    gamma = 1.
-    p = 0.
-    qx = 2.
-    qy = 2.
-    qz = 2.
-    
-    def __init__(self, mesh, mapping=None, **kwargs):
-        Simple.__init__(self, mesh, mapping=mapping, **kwargs)
-        
-    
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx, self.Wxx)
-            if self.mesh.dim > 1:
-                wlist += (self.Wy, self.Wyy)
-            if self.mesh.dim > 2:
-                wlist += (self.Wz, self.Wzz)
-            self._Wsmooth = sp.vstack(wlist)
-        return self._Wsmooth
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Ws, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-    
-    @property
-    def Ws(self):
-        """Regularization matrix Ws"""
-        if getattr(self, 'm', None) is None:
-            self.Rs = Utils.speye(self.mesh.nC)            
-            
-        else:
-            f_m = self.m
-            self.rs = self.R(f_m , self.p, self.eps)
-            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
-            self.Rs = Utils.sdiag( self.rs )
-            
-        self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs
-            
-        return self._Ws
-
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
-        
-        if getattr(self, 'm', None) is None:
-            self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0])            
-            
-        else:
-            f_m = self.mesh.unitCellGradx * self.m
-            self.rx = self.R( f_m , self.qx, self.eps)
-            self.Rx = Utils.sdiag( self.rx )
-            
-        if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx
-        return self._Wx
-
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
-        
-        if getattr(self, 'm', None) is None:
-            self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0])            
-            
-        else:
-            f_m = self.mesh.unitCellGrady * self.m
-            self.ry = self.R( f_m , self.qy, self.eps)
-            self.Ry = Utils.sdiag( self.ry )
-            
-        if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady
-        return self._Wy
-
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
-
-        if getattr(self, 'm', None) is None:
-            self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0])            
-            
-        else:
-            f_m = self.mesh.unitCellGradz * self.m
-            self.rz = self.R( f_m , self.qz, self.eps)
-            self.Rz = Utils.sdiag( self.rz )
-        
-        if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz
-        return self._Wz    
-
-    
-    def R(self, f_m , p, dec):
-
-        eta = (self.eps**(1-p/2.))**0.5
-        r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.)
-
-        return r

SimPEG/Survey.py

@@ -1,5 +1,6 @@
 import Utils, numpy as np, scipy.sparse as sp, uuid
 
+
 class BaseRx(object):
     """SimPEG Receiver Object"""
 
@@ -204,6 +205,7 @@ class BaseSurvey(object):
     __metaclass__ = Utils.SimPEGMetaClass
 
     std = None       #: Estimated Standard Deviations
+    eps = None       #: Estimated Noise Floor
     dobs = None      #: Observed data
     dtrue = None     #: True data, if data is synthetic
     mtrue = None     #: True model, if data is synthetic
@@ -222,8 +224,6 @@ class BaseSurvey(object):
 
     @srcList.setter
     def srcList(self, value):
-        if isinstance(value, self.srcPair):
-            value = [value]
         assert type(value) is list, 'srcList must be a list'
         assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of %s' % self.srcPair.__name__
         assert len(set(value)) == len(value), 'The srcList must be unique'
@@ -375,11 +375,3 @@ class BaseSurvey(object):
         self.dobs = self.dtrue+noise
         self.std = self.dobs*0 + std
         return self.dobs
-
-class LinearSurvey(BaseSurvey):
-    def projectFields(self, u):
-        return u
-    
-    @property
-    def nD(self):
-        return self.prob.G.shape[0]

SimPEG/Utils/ModelBuilder.py

@@ -118,44 +118,6 @@ def defineElipse(ccMesh, center=[0,0,0], anisotropy=[1,1,1], slope=10., theta=0.
     D = np.sqrt(np.sum(G**2,axis=1))
     return -np.arctan((D-1)*slope)*(2./np.pi)/2.+0.5
 
-def getIndicesSphere(center,radius,ccMesh):
-    """
-        Creates a vector containing the sphere indices in the cell centers mesh.
-        Returns a tuple
-
-        The sphere is defined by the points
-
-        p0, describe the position of the center of the cell
-
-        r, describe the radius of the sphere.
-
-        ccMesh represents the cell-centered mesh
-
-        The points p0 must live in the the same dimensional space as the mesh.
-
-    """
-
-    # Validation: mesh and point (p0) live in the same dimensional space
-    dimMesh = np.size(ccMesh[0,:])
-    assert len(center) == dimMesh, "Dimension mismatch. len(p0) != dimMesh"
-
-    if dimMesh == 1:
-       # Define the reference points
-        
-        ind  = np.abs(center[0] - ccMesh[:,0]) < radius
-
-    elif dimMesh == 2:
-       # Define the reference points
-
-        ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 ) < radius
-
-    elif dimMesh == 3:
-        # Define the points
-        ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 + ( center[2] - ccMesh[:,2] )**2 ) < radius
-
-    # Return a tuple
-    return ind
-
 def defineTwoLayers(ccMesh,depth,vals=[0,1]):
     """
     Define a two layered model.  Depth of the first layer must be specified.

SimPEG/Utils/matutils.py

@@ -26,6 +26,9 @@ def mkvc(x, numDims=1):
     if hasattr(x, 'tovec'):
         x = x.tovec()
 
+    if isinstance(x, Zero):
+        return x
+    
     assert isinstance(x, np.ndarray), "Vector must be a numpy array"
 
     if numDims == 1:
@@ -37,6 +40,9 @@ def mkvc(x, numDims=1):
 
 def sdiag(h):
     """Sparse diagonal matrix"""
+    if isinstance(h, Zero):
+        return h
+
     return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr")
 
 def sdInv(M):
@@ -417,6 +423,12 @@ class Zero(object):
     def __ge__(self, v):return 0 >= v
     def __gt__(self, v):return 0 > v
 
+    @property 
+    def transpose(self): return Zero()
+    
+    @property
+    def T(self): return Zero()
+
 class Identity(object):
     _positive = True
     def __init__(self, positive=True):

docs/examples/EM_FDEM_1D_Inversion.rst

@@ -0,0 +1,26 @@
+.. _examples_EM_FDEM_1D_Inversion:
+
+.. --------------------------------- ..
+..                                   ..
+..    THIS FILE IS AUTO GENEREATED   ..
+..                                   ..
+..    SimPEG/Examples/__init__.py    ..
+..                                   ..
+.. --------------------------------- ..
+
+
+EM: FDEM: 1D: Inversion
+=======================
+
+Here we will create and run a FDEM 1D inversion.
+
+
+
+.. plot::
+
+    from SimPEG import Examples
+    Examples.EM_FDEM_1D_Inversion.run()
+
+.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py
+    :language: python
+    :linenos:

tests/utils/test_Zero.py

@@ -1,5 +1,5 @@
 import unittest
-from SimPEG.Utils import Zero, Identity, sdiag
+from SimPEG.Utils import Zero, Identity, sdiag, mkvc
 from SimPEG import np, sp
 
 class Tests(unittest.TestCase):
@@ -29,6 +29,11 @@ class Tests(unittest.TestCase):
         assert a == 1
         self.assertRaises(ZeroDivisionError, lambda:3/z)
 
+        assert mkvc(z) == 0
+        assert sdiag(z)*a == 0
+        assert z.T == 0
+        assert z.transpose == 0
+
     def test_mat_zero(self):
         z = Zero()
         S = sdiag(np.r_[2,3])