wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-03 04:20:35 +08:00
Code Issues Packages Projects Releases Wiki Activity

Move cell-based weights (i.e. distance weighting) inside regularization.

Browse Source 
Fix gamma parameter update
TO DO: Check inversion print screen -> values don't match reality.
This commit is contained in:
D Fournier
2016-03-11 11:40:47 -08:00
parent 838035adae
commit 38b4079f0b
2 changed files with 48 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Directives.py
+5 -3
View File
@@ -298,6 +298,7 @@ class update_IRLS(InversionDirective):
phim_new = self.reg.eval(self.invProb.curModel)
self.gamma = self.phi_m_last / phim_new
self.reg.curModel = self.invProb.curModel
self.reg.gamma = self.gamma
def endIter(self):
@@ -315,13 +316,14 @@ class update_IRLS(InversionDirective):
self.reg.curModel = self.invProb.curModel
# Update the pre-conditioner
diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.m.size))**2.
diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.curModel.size))**2.
PC = Utils.sdiag(diagA**-1.)
self.opt.approxHinv = PC
self.reg.gamma = 1.
phim_new = self.reg.eval(self.invProb.curModel)
self.reg.gamma = self.reg.gamma * self.invProb.phi_m_last / phim_new
self.reg.gamma = self.invProb.phi_m_last / phim_new
#==============================================================================
# import pylab as plt
SimPEG/Regularization.py
+43 -17
View File
@@ -540,36 +540,40 @@ class Simple(BaseRegularization):
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")
wght = 1.
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
if isinstance(self.wght,float):
self.wght = np.ones(self.regmesh.nC) * self.wght
@property
def Ws(self):
"""Regularization matrix Ws"""
if getattr(self,'_Ws', None) is None:
self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5)
self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
return self._Ws
@property
def Wx(self):
"""Regularization matrix Wx"""
if getattr(self, '_Wx', None) is None:
self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x)**0.5)*self.regmesh.cellDiffxStencil
self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil
return self._Wx
@property
def Wy(self):
"""Regularization matrix Wy"""
if getattr(self, '_Wy', None) is None:
self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y)**0.5)*self.regmesh.cellDiffyStencil
self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil
return self._Wy
@property
def Wz(self):
"""Regularization matrix Wz"""
if getattr(self, '_Wz', None) is None:
self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z)**0.5)*self.regmesh.cellDiffzStencil
self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil
return self._Wz
@property
@@ -648,10 +652,13 @@ class Sparse(Simple):
qx = 2.
qy = 2.
qz = 2.
wght = 1.
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
if isinstance(self.wght,float):
self.wght = np.ones(self.regmesh.nC) * self.wght
@property
def Ws(self):
@@ -661,26 +668,26 @@ class Sparse(Simple):
else:
f_m = self.curModel
self.rs = self.R(f_m , self.p, self.eps)
self.rs = self.R(f_m , self.p)
#print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
self.Rs = Utils.sdiag( self.rs )
return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs
return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
@property
def Wx(self):
"""Regularization matrix Wx"""
if getattr(self, 'curModel', None) is None:
self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
else:
f_m = self.regmesh.cellDiffxStencil * self.curModel
self.rx = self.R( f_m , self.qx, self.eps)
self.rx = self.R( f_m , self.qx)
self.Rx = Utils.sdiag( self.rx )
return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma)**0.5)*self.Rx*self.regmesh.cellDiffxStencil
return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
@property
def Wy(self):
@@ -691,10 +698,10 @@ class Sparse(Simple):
else:
f_m = self.regmesh.cellDiffyStencil * self.curModel
self.ry = self.R( f_m , self.qy, self.eps)
self.ry = self.R( f_m , self.qy)
self.Ry = Utils.sdiag( self.ry )
return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma)**0.5)*self.Ry*self.regmesh.cellDiffyStencil
return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
@property
def Wz(self):
@@ -705,12 +712,31 @@ class Sparse(Simple):
else:
f_m = self.regmesh.cellDiffzStencil * self.curModel
self.rz = self.R( f_m , self.qz, self.eps)
self.rz = self.R( f_m , self.qz)
self.Rz = Utils.sdiag( self.rz )
return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma)**0.5)*self.Rz*self.regmesh.cellDiffzStencil
return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
@property
def Wsmooth(self):
"""Full smoothness regularization matrix W"""
#if getattr(self, '_Wsmooth', None) is None:
wlist = (self.Wx,)
if self.regmesh.dim > 1:
wlist += (self.Wy,)
if self.regmesh.dim > 2:
wlist += (self.Wz,)
#self._Wsmooth = sp.vstack(wlist)
return sp.vstack(wlist)
@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 sp.vstack(wlist)
def R(self, f_m , exponent):
eta = (self.eps**(1-exponent/2.))**0.5