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

Stash the regularization between each beta

Browse Source
This commit is contained in:
D Fournier
2016-05-11 23:04:14 -07:00
parent e10d6878fb
commit cd2360b815
2 changed files with 58 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Directives.py
+21 -3
View File
@@ -271,6 +271,14 @@ class Update_IRLS(InversionDirective):
self.reg.curModel = self.invProb.curModel
self.reg.gamma = self.gamma
print "Initial gamma ", np.linalg.norm(self.reg.gamma)
# Reset the regularization matrices so that it is
# recalculated with new gamma
self.reg._Wsmall = None
self.reg._Wx = None
self.reg._Wy = None
self.reg._Wz = None
self.reg._W = None
if getattr(self, 'phi_d_last', None) is None:
self.phi_d_last = self.invProb.phi_d
@@ -294,16 +302,26 @@ class Update_IRLS(InversionDirective):
# Temporarely set gamma to 1. to get raw phi_m
self.reg.gamma = 1.
# Reset the regularization matrices so that it is
# recalculated for current model
self.reg._Wsmall = None
self.reg._Wx = None
self.reg._Wy = None
self.reg._Wz = None
self.reg._W = None
# Compute new model objective function value
phim_new = self.reg.eval(self.invProb.curModel)
# Update gamma to scale the regularization between IRLS iterations
self.reg.gamma = self.phi_m_last / phim_new
# Set the weighting matrix to None so that it is recomputed next time
# it is called in the inversion
# Reset the regularization matrices again for new gamma
self.reg._Wsmall = None
self.reg._Wx = None
self.reg._Wy = None
self.reg._Wz = None
self.reg._W = None
self.reg._Wsmooth = None
class Update_lin_PreCond(InversionDirective):
"""
SimPEG/Regularization.py
+37 -30
View File
@@ -468,6 +468,7 @@ class Simple(BaseRegularization):
@property
def W(self):
"""Full regularization matrix W"""
print 'wtf why are we using W'
if getattr(self, '_W', None) is None:
wlist = (self.Wsmall, self.Wx)
if self.regmesh.dim > 1:
@@ -907,60 +908,66 @@ class Sparse(Simple):
@property
def Wsmall(self):
"""Regularization matrix Wsmall"""
if getattr(self, 'curModel', None) is None:
self.Rs = Utils.speye(self.regmesh.nC)
if getattr(self,'_Wsmall', None) is None:
if getattr(self, 'curModel', None) is None:
self.Rs = Utils.speye(self.regmesh.nC)
else:
f_m = self.mapping * (self.curModel - self.reg.mref)
self.rs = self.R(f_m , self.eps_p, self.norms[0])
#print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
self.Rs = Utils.sdiag( self.rs )
else:
f_m = self.mapping * (self.curModel - self.reg.mref)
self.rs = self.R(f_m , self.eps_p, self.norms[0])
self.Rs = Utils.sdiag( self.rs )
return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs
self._Wsmall = Utils.sdiag((self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs
return self._Wsmall
@property
def Wx(self):
"""Regularization matrix Wx"""
if getattr(self,'_Wx', None) is None:
if getattr(self, 'curModel', None) is None:
self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
if getattr(self, 'curModel', None) is None:
self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
else:
f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel)
self.rx = self.R( f_m , self.eps_q, self.norms[1])
self.Rx = Utils.sdiag( self.rx )
else:
f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel)
self.rx = self.R( f_m , self.eps_q, self.norms[1])
self.Rx = Utils.sdiag( self.rx )
self._Wx = Utils.sdiag(( self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.cell_weights))**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.cell_weights))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
return self._Wx
@property
def Wy(self):
"""Regularization matrix Wy"""
if getattr(self,'_Wy', None) is None:
if getattr(self, 'curModel', None) is None:
self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
if getattr(self, 'curModel', None) is None:
self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
else:
f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel)
self.ry = self.R( f_m , self.eps_q, self.norms[2])
self.Ry = Utils.sdiag( self.ry )
else:
f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel)
self.ry = self.R( f_m , self.eps_q, self.norms[2])
self.Ry = Utils.sdiag( self.ry )
self._Wy = Utils.sdiag((self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.cell_weights))**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.cell_weights))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
return self._Wy
@property
def Wz(self):
"""Regularization matrix Wz"""
if getattr(self,'_Wz', None) is None:
if getattr(self, 'curModel', None) is None:
self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
if getattr(self, 'curModel', None) is None:
self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
else:
f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel)
self.rz = self.R( f_m , self.eps_q, self.norms[3])
self.Rz = Utils.sdiag( self.rz )
else:
f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel)
self.rz = self.R( f_m , self.eps_q, self.norms[3])
self.Rz = Utils.sdiag( self.rz )
return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
self._Wz = Utils.sdiag((self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
return self._Wz
def R(self, f_m , eps, exponent):