From 85b55139e898475bf72aa512e1b327a6766570d0 Mon Sep 17 00:00:00 2001
From: D Fournier <fourndo@gmail.com>
Date: Thu, 28 Jan 2016 18:36:34 -0800
Subject: [PATCH] Implement simple regularization Modified the
 Optimization.ProjectedGNCG to allow active cells back in. Fix problem
 regarding the Directive.TargetMisfit --> Survey.Linear had wrong nD value

---
 SimPEG/Optimization.py   |  16 +++-
 SimPEG/Regularization.py | 160 +++++++++++++++++++++++++++++++++++++++
 SimPEG/Survey.py         |   2 +-
 3 files changed, 176 insertions(+), 2 deletions(-)

diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py
index 4f2cb062..20a74d0d 100644
--- a/SimPEG/Optimization.py
+++ b/SimPEG/Optimization.py
@@ -989,5 +989,19 @@ 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
\ No newline at end of file
diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py
index 6c5b5f97..42c5b359 100644
--- a/SimPEG/Regularization.py
+++ b/SimPEG/Regularization.py
@@ -243,3 +243,163 @@ 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 Gx(self):
+        
+        n = self.mesh.vnC
+        gx = Utils.ddx(n[0]-1)
+        gx_square = sp.vstack((gx,gx[-1,:]*-1), format="csr")
+        
+        self._Gx = Utils.kron3(Utils.speye(n[2]), Utils.speye(n[1]), gx_square)
+        
+        return self._Gx
+        
+    @property
+    def Gy(self):
+        
+        n = self.mesh.vnC
+        gy = Utils.ddx(n[1]-1)
+        gy_square = sp.vstack((gy,gy[-1,:]*-1), format="csr")
+        
+        self._Gy = Utils.kron3(Utils.speye(n[2]), gy_square, Utils.speye(n[0]))
+        
+        return self._Gy
+        
+    @property
+    def Gz(self):
+        
+        n = self.mesh.vnC
+        gz = Utils.ddx(n[2]-1)
+        gz_square = sp.vstack((gz,gz[-1,:]*-1), format="csr")
+        
+        self._Gz = Utils.kron3( gz_square , Utils.speye(n[1]), Utils.speye(n[0]))
+        
+        return self._Gz
+        
+    @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.Gx
+        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.Gy
+        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.Gz
+        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
diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py
index 98a44ee3..f8bdec47 100644
--- a/SimPEG/Survey.py
+++ b/SimPEG/Survey.py
@@ -380,4 +380,4 @@ class LinearSurvey(BaseSurvey):
     
     @property
     def nD(self):
-        return self.prob.G.shape[1]
+        return self.prob.G.shape[0]