Merge branch 'master' into parallel

SimPEG/Maps.py

@@ -285,6 +285,37 @@ class LogMap(IdentityMap):
     def inverse(self, m):
         return np.exp(Utils.mkvc(m))
 
+class FullMap(IdentityMap):
+    """
+    FullMap
+
+    Given a scalar, the FullMap maps the value to the 
+    full model space.
+    """
+
+    def __init__(self,mesh,**kwargs):
+        IdentityMap.__init__(self, mesh,**kwargs)
+
+    @property
+    def nP(self):
+        return 1
+
+    def _transform(self, m):
+        """
+            :param m: model (scalar)
+            :rtype: numpy.array
+            :return: transformed model
+        """
+        return np.ones(self.mesh.nC)*m
+
+    def deriv(self, m):
+        """
+            :param numpy.array m: model
+            :rtype: numpy.array
+            :return: derivative of transformed model
+        """
+        return np.ones([self.mesh.nC,1])     
+      
 
 class Vertical1DMap(IdentityMap):
     """Vertical1DMap

SimPEG/PropMaps.py

@@ -239,6 +239,7 @@ class PropMap(object):
             setattr(self, '%sMap'%name, mapping)
             setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP)))
             nP += mapping.nP
+        self.nP = nP 
 
     @property
     def defaultInvProp(self):

SimPEG/Regularization.py

@@ -261,22 +261,30 @@ class Tikhonov(BaseRegularization):
         return self._Wzz
 
     @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Ws, self.Wx, self.Wxx)
+    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.W * ( self.mapping * (m) )
-            r2 = self.Ws * ( self.mapping * (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.smoothModel == False:
             r = self.W * ( self.mapping * (m - self.mref) )
@@ -302,12 +310,12 @@ class Tikhonov(BaseRegularization):
         """
         if self.smoothModel == True:
             mD1 = self.mapping.deriv(m)
-            mD2 = self.mapping.deriv(self.mref)
-            r1 = self.W * ( self.mapping * (m) )
-            r2 = self.Ws * ( self.mapping * (self.mref) )
-            out1 = mD1.T * ( self.W.T * r1 )
+            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
+            out = out1+out2
         elif self.smoothModel == False:
             mD = self.mapping.deriv(m - self.mref)
             r = self.W * ( self.mapping * (m - self.mref) )

SimPEG/Tests/test_maps.py

@@ -6,8 +6,8 @@ from scipy.sparse.linalg import dsolve
 
 TOL = 1e-14
 
-MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting"]
-MAPS_TO_TEST_3D = [             "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting"]
+MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting", "FullMap"]
+MAPS_TO_TEST_3D = [             "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting", "FullMap"]
 
 class MapTests(unittest.TestCase):
 

SimPEG/Tests/test_regularization.py

@@ -22,7 +22,7 @@ class RegularizationTests(unittest.TestCase):
             mapping = r.mapPair(self.mesh2)
             reg = r(self.mesh2, mapping=mapping)
             m = np.random.rand(mapping.nP)
-            reg.mref = m[:]*0
+            reg.mref = m[:]*np.mean(m)
             passed = checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
             self.assertTrue(passed)
 

SimPEG/Utils/ModelBuilder.py

@@ -3,10 +3,27 @@ import scipy.ndimage as ndi
 import scipy.sparse as sp
 from matutils import mkvc
 
-
-def getIndecesBlock(p0,p1,ccMesh):
+def addBlock(gridCC, modelCC, p0, p1, blockProp):
     """
-        Creates a vector containing the block indexes in the cell centerd mesh.
+        Add a block to an exsisting cell centered model, modelCC
+
+        :param numpy.array, gridCC: mesh.gridCC is the cell centered grid
+        :param numpy.array, modelCC: cell centered model
+        :param numpy.array, p0: bottom, southwest corner of block
+        :param numpy.array, p1: top, northeast corner of block
+        :blockProp float, blockProp: property to assign to the model 
+
+        :return numpy.array, modelBlock: model with block
+    """
+    ind = getIndicesBlock(p0, p1, gridCC)
+    modelBlock = modelCC.copy()
+    modelBlock[ind] = blockProp
+    return modelBlock
+
+
+def getIndicesBlock(p0,p1,ccMesh):
+    """
+        Creates a vector containing the block indices in the cell centers mesh.
         Returns a tuple
 
         The block is defined by the points
@@ -78,7 +95,7 @@ def defineBlock(ccMesh,p0,p1,vals=[0,1]):
         vals[1]  conductivity of the ground
     """
     sigma = np.zeros(ccMesh.shape[0]) + vals[1]
-    ind   = getIndecesBlock(p0,p1,ccMesh)
+    ind   = getIndicesBlock(p0,p1,ccMesh)
 
     sigma[ind] = vals[0]
 
@@ -132,7 +149,7 @@ def defineTwoLayers(ccMesh,depth,vals=[0,1]):
     # The depth is always defined on the last one.
     p1[len(p1)-1] -= depth
 
-    ind   = getIndecesBlock(p0,p1,ccMesh)
+    ind   = getIndicesBlock(p0,p1,ccMesh)
 
     sigma[ind] = vals[0];