diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index c5754c93..793b2581 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -565,58 +565,7 @@ class DiffOperators(object): return Pbc, Pin, Pout - - def unitCellGradx(): - doc = """Cell centered Gradient in the x dimension used for - regularization. The gradient operator is square (nC-by-nC)""" - def fget(self): - if self.dim < 3: return None - if getattr(self, '_unitCellGradx', None) is None: - - n = self.vnC - gx = ddx(n[0]-1) - gx_square = sp.vstack((gx,gx[-1,:]*-1), format="csr") - - self._unitCellGradx = kron3(speye(n[2]), speye(n[1]), gx_square) - - return self._unitCellGradx - return locals() - unitCellGradx = property(**unitCellGradx()) - - def unitCellGrady(): - doc = """Cell centered Gradient in they dimension used for - regularization. The gradient operator is square (nC-by-nC)""" - def fget(self): - if self.dim < 3: return None - if getattr(self, '_unitCellGrady', None) is None: - - n = self.vnC - gy = ddx(n[1]-1) - gy_square = sp.vstack((gy,gy[-1,:]*-1), format="csr") - - self._unitCellGrady = kron3(speye(n[2]), gy_square, speye(n[0])) - - return self._unitCellGrady - return locals() - unitCellGrady = property(**unitCellGrady()) - - def unitCellGradz(): - doc = """Cell centered Gradient in they dimension used for - regularization. The gradient operator is square (nC-by-nC)""" - def fget(self): - if self.dim < 3: return None - if getattr(self, '_unitCellGradz', None) is None: - - n = self.vnC - gz = ddx(n[2]-1) - gz_square = sp.vstack((gz,gz[-1,:]*-1), format="csr") - - self._unitCellGradz = kron3( gz_square , speye(n[1]), speye(n[0])) - - return self._unitCellGradz - return locals() - unitCellGradz = property(**unitCellGradz()) - + # --------------- Averaging --------------------- @property @@ -797,4 +746,3 @@ class DiffOperators(object): kron3(av(n[2]), speye(n[1]+1), av(n[0])), kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr") return self._aveN2F - diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 95bb40cc..dd74b0b5 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -281,212 +281,3 @@ class Tikhonov(BaseRegularization): r = self.W * ( self.mapping * (m - self.mref) ) 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") - - def __init__(self, mesh, mapping=None, **kwargs): - BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs) - - @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.mesh.unitCellGradx - 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.mesh.unitCellGrady - 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.mesh.unitCellGradz - return self._Wz - - @property - def Wsmooth(self): - """Full smoothness regularization matrix W""" - if getattr(self, '_Wsmooth', None) is None: - wlist = (self.Wx,) - if self.mesh.dim > 1: - wlist += (self.Wy,) - if self.mesh.dim > 2: - wlist += (self.Wz,) - 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 - -class SparseRegularization(Simple): - - eps = 1e-1 - m = None - gamma = 1. - p = 0. - qx = 2. - qy = 2. - qz = 2. - - def __init__(self, mesh, mapping=None, **kwargs): - Simple.__init__(self, mesh, mapping=mapping, **kwargs) - - @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 - - @property - def Ws(self): - """Regularization matrix Ws""" - if getattr(self, 'm', None) is None: - self.Rs = Utils.speye(self.mesh.nC) - - else: - f_m = self.m - self.rs = self.R(f_m , self.p, self.eps) - #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs)) - self.Rs = Utils.sdiag( self.rs ) - - self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs - - return self._Ws - - @property - def Wx(self): - """Regularization matrix Wx""" - - if getattr(self, 'm', None) is None: - self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0]) - - else: - f_m = self.mesh.unitCellGradx * self.m - self.rx = self.R( f_m , self.qx, self.eps) - self.Rx = Utils.sdiag( self.rx ) - - if getattr(self, '_Wx', None) is None: - self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx - return self._Wx - - @property - def Wy(self): - """Regularization matrix Wy""" - - if getattr(self, 'm', None) is None: - self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0]) - - else: - f_m = self.mesh.unitCellGrady * self.m - self.ry = self.R( f_m , self.qy, self.eps) - self.Ry = Utils.sdiag( self.ry ) - - if getattr(self, '_Wy', None) is None: - self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady - return self._Wy - - @property - def Wz(self): - """Regularization matrix Wz""" - - if getattr(self, 'm', None) is None: - self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0]) - - else: - f_m = self.mesh.unitCellGradz * self.m - self.rz = self.R( f_m , self.qz, self.eps) - self.Rz = Utils.sdiag( self.rz ) - - if getattr(self, '_Wz', None) is None: - self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz - return self._Wz - - - def R(self, f_m , p, dec): - - eta = (self.eps**(1-p/2.))**0.5 - r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.) - - return r