diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 74a92efd..6dca13dd 100644 --- a/SimPEG/Maps.py +++ b/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 diff --git a/SimPEG/PropMaps.py b/SimPEG/PropMaps.py index 995216f7..527a6f7e 100644 --- a/SimPEG/PropMaps.py +++ b/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): diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 9e90a2e9..e3506290 100644 --- a/SimPEG/Regularization.py +++ b/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) ) diff --git a/SimPEG/Tests/test_maps.py b/SimPEG/Tests/test_maps.py index a80bc40b..a4ebd80c 100644 --- a/SimPEG/Tests/test_maps.py +++ b/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): diff --git a/SimPEG/Tests/test_regularization.py b/SimPEG/Tests/test_regularization.py index 2197d27b..48846d2f 100644 --- a/SimPEG/Tests/test_regularization.py +++ b/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) diff --git a/SimPEG/Utils/ModelBuilder.py b/SimPEG/Utils/ModelBuilder.py index f17dac89..35ff3a5b 100644 --- a/SimPEG/Utils/ModelBuilder.py +++ b/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]; diff --git a/docs/InversionWorkflow.png b/docs/InversionWorkflow.png index b1081975..ee8f8146 100644 Binary files a/docs/InversionWorkflow.png and b/docs/InversionWorkflow.png differ