Minor changes...test runs fine from here

SimPEG/Directives.py

@@ -271,7 +271,8 @@ class Update_IRLS(InversionDirective):
 
         self.reg.curModel = self.invProb.curModel
         self.reg.gamma = self.gamma
-
+        
+        self.reg._W = None
         if getattr(self, 'phi_d_last', None) is None:
             self.phi_d_last = self.invProb.phi_d
 
@@ -287,7 +288,8 @@ class Update_IRLS(InversionDirective):
 
         # Get phi_m at the end of current iteration
         self.phi_m_last = self.invProb.phi_m_last
-
+        
+        self.reg._W = None
          # Update the model used for the IRLS weights
         self.reg.curModel = self.invProb.curModel
 
@@ -296,24 +298,37 @@ class Update_IRLS(InversionDirective):
 
         # Compute new model objective function value
         phim_new = self.reg.eval(self.invProb.curModel)
-
+        
         # Update gamma to scale the regularization between IRLS iterations
         self.reg.gamma = self.phi_m_last / phim_new
-
+        
         # Set the weighting matrix to None so that it is recomputed next time
         # it is called in the inversion
         self.reg._W = None
 
 class Update_lin_PreCond(InversionDirective):
-
-
+    """
+    Create a Jacobi preconditioner for the linear problem
+    """
+    onlyOnStart=False
+    
+    def initialize(self):
+    
+        if getattr(self.opt, 'approxHinv', None) is None:
+            # Update the pre-conditioner
+            diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
+            self.opt.approxHinv = PC
+            
     def endIter(self):
         # Cool the threshold parameter
-
+        if self.onlyOnStart==True:
+            return
+            
         if getattr(self.opt, 'approxHinv', None) is not None:
             # Update the pre-conditioner
             diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
-            PC     = Utils.sdiag(diagA**-1.)
+            PC     = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.)
             self.opt.approxHinv = PC
 
 

SimPEG/Regularization.py

@@ -663,9 +663,8 @@ class Sparse(Simple):
             self.Rs = Utils.speye(self.regmesh.nC)
 
         else:
-            f_m = self.curModel - self.reg.mref
+            f_m = self.mapping * (self.curModel - self.reg.mref)
             self.rs = self.R(f_m , self.eps_p, self.norms[0])
-            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
             self.Rs = Utils.sdiag( self.rs )
 
         return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
@@ -679,7 +678,7 @@ class Sparse(Simple):
             self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
 
         else:
-            f_m = self.regmesh.cellDiffxStencil * self.curModel
+            f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel)
             self.rx = self.R( f_m , self.eps_q, self.norms[1])
             self.Rx = Utils.sdiag( self.rx )
 
@@ -693,7 +692,7 @@ class Sparse(Simple):
             self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
 
         else:
-            f_m = self.regmesh.cellDiffyStencil * self.curModel
+            f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel)
             self.ry = self.R( f_m , self.eps_q, self.norms[2])
             self.Ry = Utils.sdiag( self.ry )
 
@@ -707,7 +706,7 @@ class Sparse(Simple):
             self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
 
         else:
-            f_m = self.regmesh.cellDiffzStencil * self.curModel
+            f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel)
             self.rz = self.R( f_m , self.eps_q, self.norms[3])
             self.Rz = Utils.sdiag( self.rz )
 
@@ -729,6 +728,7 @@ class Sparse(Simple):
     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