wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-02 12:32:09 +08:00
Code Issues Packages Projects Releases Wiki Activity

Remove second derivatives in simple regularization

Browse Source
This commit is contained in:
Rowan Cockett
2016-02-16 13:01:32 -08:00
parent e4c30705ff
commit 1d06b768e2
1 changed files with 35 additions and 63 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Regularization.py
+35 -63
View File
@@ -285,22 +285,17 @@ class Tikhonov(BaseRegularization):
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")
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"""
@@ -329,36 +324,15 @@ class Simple(BaseRegularization):
self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz
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)
wlist = (self.Wx,)
if self.mesh.dim > 1:
wlist += (self.Wy, self.Wyy)
wlist += (self.Wy,)
if self.mesh.dim > 2:
wlist += (self.Wz, self.Wzz)
wlist += (self.Wz,)
self._Wsmooth = sp.vstack(wlist)
return self._Wsmooth
@@ -413,20 +387,18 @@ class Simple(BaseRegularization):
return out
class SparseRegularization(Simple):
eps = 1e-1
m = None
eps = 1e-1
m = None
gamma = 1.
p = 0.
qx = 2.
qy = 2.
qz = 2.
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"""
@@ -446,35 +418,35 @@ class SparseRegularization(Simple):
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)
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])
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
@@ -482,15 +454,15 @@ class SparseRegularization(Simple):
@property
def Wy(self):
"""Regularization matrix Wy"""
if getattr(self, 'm', None) is None:
self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0])
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
@@ -500,18 +472,18 @@ class SparseRegularization(Simple):
"""Regularization matrix Wz"""
if getattr(self, 'm', None) is None:
self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0])
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
return self._Wz
def R(self, f_m , p, dec):
eta = (self.eps**(1-p/2.))**0.5