Update directives

Add IRLS example

SimPEG/Directives.py

@@ -237,39 +237,41 @@ 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 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):
 
@@ -295,6 +297,8 @@ class update_IRLS(InversionDirective):
 
         # 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
@@ -320,22 +324,39 @@ class update_IRLS(InversionDirective):
         
          # Update the model used for the IRLS weights
         self.reg.curModel = 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.reg.curModel.size))**2.
-        PC     = Utils.sdiag(diagA**-1.)
-        self.opt.approxHinv = PC
          
         # Temporarely set gamma to 1.
         self.reg.gamma = 1.
         
         # Compute change in model objective function and update scaling
         phim_new = self.reg.eval(self.invProb.curModel)
+
         self.reg.gamma = self.phi_m_last / phim_new
         
-        # 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
         
+        # TO DO: Check optimization class, data misfit not matching reality        
+        #dpred = self.prob.fields(self.invProb.curModel)
+        #phid = self.invProb.dmisfit.eval(self.invProb.curModel)
+        #print self.survey.std[0]
+        #print phid
+        #print self.invProb.phi_d
+        #print self.invProb.phi_d_last
+        
+        self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
+
+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
+            print 'Updated pre-cond'
+    
 #==============================================================================
 #          import pylab as plt
 #          plt.figure()

SimPEG/Examples/Inversion_IRLS.py

@@ -1,7 +1,7 @@
 from SimPEG import *
 
 
-def run(N=100, plotIt=True):
+def run(N=200, plotIt=True):
     """
         Inversion: Linear Problem
         =========================
@@ -9,39 +9,21 @@ def run(N=100, plotIt=True):
         Here we go over the basics of creating a linear problem and inversion.
 
     """
-
-    class LinearSurvey(Survey.BaseSurvey):
-        def projectFields(self, u):
-            return u
-
-    class LinearProblem(Problem.BaseProblem):
-
-        surveyPair = LinearSurvey
-
-        def __init__(self, mesh, G, **kwargs):
-            Problem.BaseProblem.__init__(self, mesh, **kwargs)
-            self.G = G
-
-        def fields(self, m, u=None):
-            return self.G.dot(m)
-
-        def Jvec(self, m, v, u=None):
-            return self.G.dot(v)
-
-        def Jtvec(self, m, v, u=None):
-            return self.G.T.dot(v)
-
-
+    
+    
     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.
 
-    nk = 20
-    jk = np.linspace(1.,20.,nk)
-    p = -0.25
-    q = 0.25
-
-    g = lambda k: np.exp(p*jk[k]*mesh.vectorCCx)*np.cos(2*np.pi*q*jk[k]*mesh.vectorCCx)
+    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))
 
@@ -52,25 +34,43 @@ def run(N=100, plotIt=True):
     mtrue[mesh.vectorCCx > 0.3] = 1.
     mtrue[mesh.vectorCCx > 0.45] = -0.5
     mtrue[mesh.vectorCCx > 0.6] = 0
+  
 
-    prob = LinearProblem(mesh, G)
-    survey = LinearSurvey()
+    prob = Problem.LinearProblem(mesh, G)
+    survey = Survey.LinearSurvey()
     survey.pair(prob)
-    survey.makeSyntheticData(mtrue, std=0.01)
-
-    M = prob.mesh
-
-    reg = Regularization.Tikhonov(mesh)
+    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) )**0
+    
+    reg = Regularization.Simple(mesh)
+    reg.wght = wr
+    
     dmis = DataMisfit.l2_DataMisfit(survey)
-    opt = Optimization.ProjectedGNCG(maxIter=20,lower=-1.,upper=1., maxIterCG= 20, tolCG = 1e-3)
+    dmis.Wd = 1./wd
+    
+    opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4)
     invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
-    beta = Directives.BetaSchedule()
+    invProb.curModel = m0
+    
+    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
+    target = Directives.TargetMisfit()
+
     betaest = Directives.BetaEstimate_ByEig()
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest])
-    m0 = np.zeros_like(survey.mtrue)
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
+    
 
     mrec = inv.run(m0)
-
+    
+    print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) 
+    
     if plotIt:
         import matplotlib.pyplot as plt
 
@@ -78,12 +78,54 @@ def run(N=100, plotIt=True):
         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(M.vectorCCx, survey.mtrue, 'b-')
-        axes[1].plot(M.vectorCCx, mrec, 'r-')
-        axes[1].legend(('True Model', 'Recovered Model'))
+        
+        axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
+        axes[1].plot(mesh.vectorCCx, mrec, 'r-')
+        #axes[1].legend(('True Model', 'Recovered Model'))
+        axes[1].set_ylim(-1.0,1.25)
         plt.show()
+        
+    # 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 = 2e-3
+    reg.eps_q = 2e-3
+    reg.norms   = [0., 0., 2., 2.]
+    reg.wght = wr
+    
+    opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterCG= 200, 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:
+        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)
+        
     return prob, survey, mesh, mrec
 
 if __name__ == '__main__':

SimPEG/Regularization.py

@@ -643,13 +643,11 @@ class Simple(Tikhonov):
 class Sparse(Simple):
 
     # set default values
-    eps      = 1e-1
+    eps_p      = 1e-1
+    eps_q      = 1e-1
     curModel = None # use a model to compute the weights
     gamma    = 1.
-    norms    = [0., .2, 2., 2., 1.]
-    qx       = 2.
-    qy       = 2.
-    qz       = 2.
+    norms    = [0., 2., 2., 2.]
     wght     = 1.
 
     def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
@@ -666,7 +664,7 @@ class Sparse(Simple):
 
         else:
             f_m = self.curModel - self.reg.mref
-            self.rs = self.R(f_m , self.norms[0])
+            self.rs = self.R(f_m , self.eps_p, self.norms[0])
             #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
             self.Rs = Utils.sdiag( self.rs )
 
@@ -682,7 +680,7 @@ class Sparse(Simple):
 
         else:
             f_m = self.regmesh.cellDiffxStencil * self.curModel
-            self.rx = self.R( f_m , self.qx)
+            self.rx = self.R( f_m , self.eps_q, self.norms[1])
             self.Rx = Utils.sdiag( self.rx )
 
         return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
@@ -696,7 +694,7 @@ class Sparse(Simple):
 
         else:
             f_m = self.regmesh.cellDiffyStencil * self.curModel
-            self.ry = self.R( f_m , self.qy)
+            self.ry = self.R( f_m , self.eps_q, self.norms[2])
             self.Ry = Utils.sdiag( self.ry )
 
         return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
@@ -710,7 +708,7 @@ class Sparse(Simple):
 
         else:
             f_m = self.regmesh.cellDiffzStencil * self.curModel
-            self.rz = self.R( f_m , self.qz)
+            self.rz = self.R( f_m , self.eps_q, self.norms[3])
             self.Rz = Utils.sdiag( self.rz )
 
         return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
@@ -735,9 +733,9 @@ class Sparse(Simple):
         #self._W = sp.vstack(wlist)
         return sp.vstack(wlist)
 
-    def R(self, f_m , exponent):
+    def R(self, f_m , eps, exponent):
 
-        eta = (self.eps**(1-exponent/2.))**0.5
-        r = eta / (f_m**2.+self.eps**2.)**((1-exponent/2.)/2.)
+        eta = (eps**(1-exponent/2.))**0.5
+        r = eta / (f_m**2.+ eps**2.)**((1-exponent/2.)/2.)
 
         return r

SimPEG/Survey.py

@@ -372,7 +372,7 @@ class BaseSurvey(object):
         self.dtrue = self.dpred(m, u=u)
         noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape)
         self.dobs = self.dtrue+noise
-        self.std = self.dobs*0 + std
+        self.std = self.dobs*0. + std
         return self.dobs
 
 class LinearSurvey(BaseSurvey):