mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-01 15:45:17 +08:00
Create a Simple and a SparseRegularization class.
The SparseRegularization class allows implementation of p-q norms.
This commit is contained in:
Rowan Cockett
@@ -282,3 +282,213 @@ 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")
|
||||
|
||||
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
|
||||
|
||||
Reference in New Issue
Block a user