From 38b4079f0bdb23eaa7b5810ed925e2a4f32102ce Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 11 Mar 2016 11:40:47 -0800 Subject: [PATCH] Move cell-based weights (i.e. distance weighting) inside regularization. Fix gamma parameter update TO DO: Check inversion print screen -> values don't match reality. --- SimPEG/Directives.py | 8 ++++-- SimPEG/Regularization.py | 60 ++++++++++++++++++++++++++++------------ 2 files changed, 48 insertions(+), 20 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 46b0b087..dbb525f7 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -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 diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 6ad94caf..6c6b4b4b 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -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