From fbb8cf2731e10000f71e2056215a363ffc0d4023 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 4 May 2016 23:17:01 -0700 Subject: [PATCH 01/17] modularizing regularization --- SimPEG/Regularization.py | 404 ++++++++++++++++++++++++++------------- 1 file changed, 269 insertions(+), 135 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index fc101a61..fd974130 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -403,7 +403,203 @@ class BaseRegularization(object): return mD.T * ( self.W.T * ( self.W * ( mD * v) ) ) -class Tikhonov(BaseRegularization): +class Simple(BaseRegularization): + """ + Simple regularization that does not include length scales in the derivatives. + """ + + mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options + alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Wsmall'], "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") + 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 Wsmall(self): + """Regularization matrix Wsmall""" + if getattr(self,'_Wsmall', None) is None: + self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5) + return self._Wsmall + + @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*(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*(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*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil + return self._Wz + + @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 self._Wsmooth + + @property + def W(self): + """Full regularization matrix W""" + if getattr(self, '_W', None) is None: + wlist = (self.Wsmall, self.Wsmooth) + self._W = sp.vstack(wlist) + return self._W + + @Utils.timeIt + def _evalSmall(self, m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return 0.5 * r.dot(r) + + @property + def W(self): + """Full regularization matrix W""" + if getattr(self, '_W', None) is None: + wlist = (self.Wsmall, self.Wsmooth) + self._W = sp.vstack(wlist) + return self._W + + @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 self._Wsmooth + + + @Utils.timeIt + def _evalSmall(self, m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmallDeriv(self,m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) + + @Utils.timeIt + def _evalSmoothx(self, m): + if self.mrefInSmooth == True: + r = self.Wx * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wx * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmoothy(self, m): + if self.mrefInSmooth == True: + r = self.Wy * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wy * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmoothz(self, m): + if self.mrefInSmooth == True: + r = self.Wz * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wz * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmooth(self, m): + phiSmooth = self._evalSmoothx(m) + if self.regmesh.dim > 1: + phiSmooth += self._evalSmoothy(m) + if self.regmesh.dim > 2: + phiSmooth += self._evalSmoothz(m) + return phiSmooth + + @Utils.timeIt + def _evalSmoothxDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wx * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wx * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wx * ( self.mapping * m ) + return r.T * ( self.Wx * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmoothyDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wy * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wy * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wy * ( self.mapping * m ) + return r.T * ( self.Wy * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmoothzDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wz * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wz * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wz * ( self.mapping * m ) + return r.T * ( self.Wz * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmoothDeriv(self,m): + deriv = self._evalSmoothxDeriv(m) + if self.regmesh.dim > 1: + deriv += self._evalSmoothyDeriv(m) + if self.regmesh.dim > 2: + deriv += self._evalSmoothzDeriv(m) + return deriv + + + @Utils.timeIt + def eval(self, m): + return self._evalSmall(m) + self._evalSmooth(m) + + @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})} + + """ + return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + + + +class Tikhonov(Simple): """ L2 Tikhonov regularization with both smallness and smoothness (first order derivative) contributions. @@ -493,56 +689,92 @@ class Tikhonov(BaseRegularization): self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz return self._Wzz + @property - def Wsmooth(self): + def Wsmooth2(self): """Full smoothness regularization matrix W""" if getattr(self, '_Wsmooth', None) is None: - wlist = (self.Wx, self.Wxx) + wlist = (self.Wxx) if self.regmesh.dim > 1: - wlist += (self.Wy, self.Wyy) + wlist += (self.Wyy) if self.regmesh.dim > 2: - wlist += (self.Wz, self.Wzz) + wlist += (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.Wsmall, self.Wsmooth) - self._W = sp.vstack(wlist) - return self._W - @Utils.timeIt - def _evalSmall(self, m): - r = self.Wsmall * ( self.mapping * (m - self.mref) ) - return 0.5 * r.dot(r) - - @Utils.timeIt - def _evalSmooth(self, m): + def _evalSmoothxx(self, m): if self.mrefInSmooth == True: - r = self.Wsmooth * ( self.mapping * (m - self.mref) ) + r = self.Wxx * ( self.mapping * (m - self.mref) ) elif self.mrefInSmooth == False: - r = self.Wsmooth * ( self.mapping * (m) ) + r = self.Wxx * ( self.mapping * (m) ) return 0.5 * r.dot(r) + @Utils.timeIt + def _evalSmoothyy(self, m): + if self.mrefInSmooth == True: + r = self.Wyy * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wyy * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmoothzz(self, m): + if self.mrefInSmooth == True: + r = self.Wzz * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wzz * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmooth2(self, m): + phiSmooth2 = self._evalSmoothxx(m) + if self.regmesh.dim > 1: + phiSmooth2 += self._evalSmoothyy(m) + if self.regmesh.dim > 2: + phiSmooth2 += self._evalSmoothzz(m) + return phiSmooth2 + + @Utils.timeIt + def _evalSmoothxxDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wxx * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wxx * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wxx * ( self.mapping * m ) + return r.T * ( self.Wxx * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmoothyyDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wyy * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wyy * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wyy * ( self.mapping * m ) + return r.T * ( self.Wyy * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmoothzzDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wzz * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wzz * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wzz * ( self.mapping * m ) + return r.T * ( self.Wzz * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmooth2Deriv(self,m): + deriv = self._evalSmoothxxDeriv(m) + if self.regmesh.dim > 1: + deriv += self._evalSmoothyyDeriv(m) + if self.regmesh.dim > 2: + deriv += self._evalSmoothzzDeriv(m) + return deriv + + @Utils.timeIt def eval(self, m): - return self._evalSmall(m) + self._evalSmooth(m) - - @Utils.timeIt - def _evalSmallDeriv(self,m): - r = self.Wsmall * ( self.mapping * (m - self.mref) ) - return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) - - @Utils.timeIt - def _evalSmoothDeriv(self,m): - if self.mrefInSmooth == True: - r = self.Wsmooth * ( self.mapping * ( m - self.mref ) ) - return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) ) - elif self.mrefInSmooth == False: - r = self.Wsmooth * ( self.mapping * m ) - return r.T * ( self.Wsmooth * self.mapping.deriv(m) ) + return self._evalSmall(m) + self._evalSmooth(m) + self._evalSmooth2(m) @Utils.timeIt def evalDeriv(self, m): @@ -560,88 +792,9 @@ class Tikhonov(BaseRegularization): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmooth2Deriv(m) -class Simple(Tikhonov): - """ - Simple regularization that does not include length scales in the derivatives. - """ - - mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options - alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "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") - 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 Wsmall(self): - """Regularization matrix Wsmall""" - if getattr(self,'_Wsmall', None) is None: - self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5) - return self._Wsmall - - @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*(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*(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*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil - return self._Wz - - @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 self._Wsmooth - - @property - def W(self): - """Full regularization matrix W""" - if getattr(self, '_W', None) is None: - wlist = (self.Wsmall, self.Wsmooth) - self._W = sp.vstack(wlist) - return self._W - - @Utils.timeIt - def _evalSmall(self, m): - r = self.Wsmall * ( self.mapping * (m - self.mref) ) - return 0.5 * r.dot(r) - - @Utils.timeIt - def _evalSmooth(self, m): - if self.mrefInSmooth == True: - r = self.Wsmooth * ( self.mapping * (m - self.mref) ) - elif self.mrefInSmooth == False: - r = self.Wsmooth * ( self.mapping * m) - return 0.5 * r.dot(r) - class Sparse(Simple): @@ -716,25 +869,6 @@ class Sparse(Simple): 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.Wsmall, self.Wsmooth) - self._W = sp.vstack(wlist) - return self._W def R(self, f_m , eps, exponent): From b4ab60c2601b6cb2dd64a65d58863f66f265f458 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Thu, 5 May 2016 11:55:56 -0700 Subject: [PATCH 02/17] Add model mapping to sparse regularization --- SimPEG/Regularization.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index ed1039ec..1b0d180f 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -663,7 +663,7 @@ class Sparse(Simple): self.Rs = Utils.speye(self.regmesh.nC) else: - f_m = self.curModel - self.reg.mref + 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 ) @@ -679,7 +679,7 @@ class Sparse(Simple): self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0]) else: - f_m = self.regmesh.cellDiffxStencil * self.curModel + 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 ) @@ -693,7 +693,7 @@ class Sparse(Simple): self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0]) else: - f_m = self.regmesh.cellDiffyStencil * self.curModel + 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 ) @@ -707,7 +707,7 @@ class Sparse(Simple): self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0]) else: - f_m = self.regmesh.cellDiffzStencil * self.curModel + 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 ) From fb5434695f13388790df41b6f4b293a5425d053d Mon Sep 17 00:00:00 2001 From: D Fournier Date: Tue, 10 May 2016 14:31:45 -0700 Subject: [PATCH 03/17] Alpha_s default to 1.0 --- SimPEG/Regularization.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 7f50ea59..8fdb9fd5 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -409,7 +409,7 @@ class Simple(BaseRegularization): """ mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options - alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Wsmall'], "Smallness weight") + alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "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") From eaa37f42e47b3e1d01ed10fd50d9b7cd8d29796f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 14:51:56 -0700 Subject: [PATCH 04/17] remove duplicate evalSmall --- SimPEG/Regularization.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 8fdb9fd5..1276cbb6 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -469,10 +469,6 @@ class Simple(BaseRegularization): self._W = sp.vstack(wlist) return self._W - @Utils.timeIt - def _evalSmall(self, m): - r = self.Wsmall * ( self.mapping * (m - self.mref) ) - return 0.5 * r.dot(r) @property def W(self): From 3f0c89f10b3783e63b14ccb530b565613d0cbff0 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 14:53:43 -0700 Subject: [PATCH 05/17] remove extra Ws --- SimPEG/Regularization.py | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 1276cbb6..361a3206 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -470,27 +470,6 @@ class Simple(BaseRegularization): return self._W - @property - def W(self): - """Full regularization matrix W""" - if getattr(self, '_W', None) is None: - wlist = (self.Wsmall, self.Wsmooth) - self._W = sp.vstack(wlist) - return self._W - - @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 self._Wsmooth - - @Utils.timeIt def _evalSmall(self, m): r = self.Wsmall * ( self.mapping * (m - self.mref) ) From 2a802c1aa3ae90b120cb828194c36dfd2d57ccbf Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 16:39:47 -0700 Subject: [PATCH 06/17] weights --> cell_weights, removed vol term from simple regularization --- SimPEG/Regularization.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 361a3206..5d40adf8 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -413,40 +413,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. + cell_weights = 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 + if isinstance(self.cell_weights,float): + self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights @property def Wsmall(self): """Regularization matrix Wsmall""" if getattr(self,'_Wsmall', None) is None: - self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5) + self._Wsmall = Utils.sdiag((self.alpha_s*self.cell_weights)**0.5) return self._Wsmall @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*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil + self._Wx = Utils.sdiag((self.alpha_x * (self.regmesh.aveCC2Fx*self.cell_weights))**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*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil + self._Wy = Utils.sdiag((self.alpha_y * (self.regmesh.aveCC2Fy*self.cell_weights))**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*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil + self._Wz = Utils.sdiag((self.alpha_z * (self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.regmesh.cellDiffzStencil return self._Wz @property @@ -779,13 +779,13 @@ class Sparse(Simple): curModel = None # use a model to compute the weights gamma = 1. norms = [0., 2., 2., 2.] - wght = 1. + cell_weights = 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 + if isinstance(self.cell_weights,float): + self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights @property def Wsmall(self): @@ -799,7 +799,7 @@ class Sparse(Simple): #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*self.wght)**0.5)*self.Rs + return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs @property @@ -814,7 +814,7 @@ class Sparse(Simple): self.rx = self.R( f_m , self.eps_q, self.norms[1]) self.Rx = Utils.sdiag( self.rx ) - 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 + 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 @property def Wy(self): @@ -828,7 +828,7 @@ class Sparse(Simple): self.ry = self.R( f_m , self.eps_q, self.norms[2]) self.Ry = Utils.sdiag( self.ry ) - 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 + 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 @property def Wz(self): @@ -842,7 +842,7 @@ class Sparse(Simple): 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.wght))**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.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil def R(self, f_m , eps, exponent): From 955bd54019b9e6fccf8f766ce83d4fdde9ac61f0 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 16:46:11 -0700 Subject: [PATCH 07/17] notation cleanup in Regularization --- SimPEG/Regularization.py | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 5d40adf8..16bb4b9f 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -1,4 +1,6 @@ -import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp +import Utils, Maps, Mesh +import numpy as np +import scipy.sparse as sp class RegularizationMesh(object): """ @@ -476,7 +478,7 @@ class Simple(BaseRegularization): return 0.5 * r.dot(r) @Utils.timeIt - def _evalSmallDeriv(self,m): + def _evalSmallDeriv(self, m): r = self.Wsmall * ( self.mapping * (m - self.mref) ) return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) @@ -514,7 +516,7 @@ class Simple(BaseRegularization): return phiSmooth @Utils.timeIt - def _evalSmoothxDeriv(self,m): + def _evalSmoothxDeriv(self, m): if self.mrefInSmooth == True: r = self.Wx * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wx * self.mapping.deriv(m - self.mref) ) @@ -523,7 +525,7 @@ class Simple(BaseRegularization): return r.T * ( self.Wx * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmoothyDeriv(self,m): + def _evalSmoothyDeriv(self, m): if self.mrefInSmooth == True: r = self.Wy * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wy * self.mapping.deriv(m - self.mref) ) @@ -532,7 +534,7 @@ class Simple(BaseRegularization): return r.T * ( self.Wy * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmoothzDeriv(self,m): + def _evalSmoothzDeriv(self, m): if self.mrefInSmooth == True: r = self.Wz * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wz * self.mapping.deriv(m - self.mref) ) @@ -541,7 +543,7 @@ class Simple(BaseRegularization): return r.T * ( self.Wz * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmoothDeriv(self,m): + def _evalSmoothDeriv(self, m): deriv = self._evalSmoothxDeriv(m) if self.regmesh.dim > 1: deriv += self._evalSmoothyDeriv(m) @@ -711,7 +713,7 @@ class Tikhonov(Simple): return phiSmooth2 @Utils.timeIt - def _evalSmoothxxDeriv(self,m): + def _evalSmoothxxDeriv(self, m): if self.mrefInSmooth == True: r = self.Wxx * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wxx * self.mapping.deriv(m - self.mref) ) @@ -720,7 +722,7 @@ class Tikhonov(Simple): return r.T * ( self.Wxx * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmoothyyDeriv(self,m): + def _evalSmoothyyDeriv(self, m): if self.mrefInSmooth == True: r = self.Wyy * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wyy * self.mapping.deriv(m - self.mref) ) @@ -729,7 +731,7 @@ class Tikhonov(Simple): return r.T * ( self.Wyy * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmoothzzDeriv(self,m): + def _evalSmoothzzDeriv(self, m): if self.mrefInSmooth == True: r = self.Wzz * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wzz * self.mapping.deriv(m - self.mref) ) @@ -738,7 +740,7 @@ class Tikhonov(Simple): return r.T * ( self.Wzz * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmooth2Deriv(self,m): + def _evalSmooth2Deriv(self, m): deriv = self._evalSmoothxxDeriv(m) if self.regmesh.dim > 1: deriv += self._evalSmoothyyDeriv(m) @@ -774,12 +776,12 @@ class Tikhonov(Simple): class Sparse(Simple): # set default values - eps_p = 1e-1 - eps_q = 1e-1 + eps_p = 1e-1 + eps_q = 1e-1 curModel = None # use a model to compute the weights - gamma = 1. - norms = [0., 2., 2., 2.] - cell_weights = 1. + gamma = 1. + norms = [0., 2., 2., 2.] + cell_weights = 1. def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) From 7964ebce50f86b87f11873cfb0e40632f4eb930f Mon Sep 17 00:00:00 2001 From: D Fournier Date: Tue, 10 May 2016 17:20:46 -0700 Subject: [PATCH 08/17] Update directive to None the Wsmooth after iteration. --- SimPEG/Directives.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 4fc1ffc5..008cc57b 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -303,26 +303,27 @@ class Update_IRLS(InversionDirective): # Set the weighting matrix to None so that it is recomputed next time # it is called in the inversion self.reg._W = None + self.reg._Wsmooth = None class Update_lin_PreCond(InversionDirective): """ Create a Jacobi preconditioner for the linear problem """ onlyOnStart=False - + def initialize(self): - + if getattr(self.opt, 'approxHinv', None) is None: # 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.curModel.size))**2. PC = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.) self.opt.approxHinv = PC - + def endIter(self): # Cool the threshold parameter if self.onlyOnStart==True: return - + if getattr(self.opt, 'approxHinv', None) is not None: # 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.curModel.size))**2. From 90a3030796234798c3755745810415f52ef5a77e Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 21:39:15 -0700 Subject: [PATCH 09/17] fixed 2 deriv --- SimPEG/Examples/Inversion_IRLS.py | 4 +- SimPEG/Regularization.py | 123 ++++++++++++++++++++++++++++-- 2 files changed, 119 insertions(+), 8 deletions(-) diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index 06ef4be4..c925cce2 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -84,12 +84,12 @@ def run(N=200, plotIt=True): #============================================================================== #reg.recModel = mrec - reg.wght = np.ones(mesh.nC) + # reg.cell_weight = np.ones(mesh.nC) reg.mref = np.zeros(mesh.nC) reg.eps_p = 5e-2 reg.eps_q = 1e-2 reg.norms = [0., 0., 2., 2.] - reg.wght = wr + reg.cell_weight = wr opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 16bb4b9f..5947ed03 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -410,11 +410,11 @@ class Simple(BaseRegularization): Simple regularization that does not include length scales in the derivatives. """ - mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options + mrefInSmooth = False #: include mref in the smoothness? alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "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_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") cell_weights = 1. def __init__(self, mesh, mapping=None, indActive=None, **kwargs): @@ -454,6 +454,8 @@ class Simple(BaseRegularization): @property def Wsmooth(self): """Full smoothness regularization matrix W""" + print 'wtf why are we using Wsmooth' + raise NotImplementedError if getattr(self, '_Wsmooth', None) is None: wlist = (self.Wx,) if self.regmesh.dim > 1: @@ -482,6 +484,13 @@ class Simple(BaseRegularization): r = self.Wsmall * ( self.mapping * (m - self.mref) ) return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) + @Utils.timeIt + def _evalSmall2Deriv(self, m, v = None): + rDeriv = self.Wsmall * ( self.mapping.deriv(m - self.mref) ) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothx(self, m): if self.mrefInSmooth == True: @@ -524,6 +533,17 @@ class Simple(BaseRegularization): r = self.Wx * ( self.mapping * m ) return r.T * ( self.Wx * self.mapping.deriv(m) ) + @Utils.timeIt + def _evalSmoothx2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wx * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wx * ( self.mapping.deriv(m) ) + + if v is not None: + return rDeriv.T * ( rDeriv * v ) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothyDeriv(self, m): if self.mrefInSmooth == True: @@ -533,6 +553,17 @@ class Simple(BaseRegularization): r = self.Wy * ( self.mapping * m ) return r.T * ( self.Wy * self.mapping.deriv(m) ) + @Utils.timeIt + def _evalSmoothy2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wy * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wy * ( self.mapping.deriv(m) ) + + if v is not None: + return rDeriv.T * ( rDeriv * v ) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothzDeriv(self, m): if self.mrefInSmooth == True: @@ -542,6 +573,17 @@ class Simple(BaseRegularization): r = self.Wz * ( self.mapping * m ) return r.T * ( self.Wz * self.mapping.deriv(m) ) + @Utils.timeIt + def _evalSmoothz2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wz * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wz * ( self.mapping.deriv(m) ) + + if v is not None: + return rDeriv.T * ( rDeriv * v ) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothDeriv(self, m): deriv = self._evalSmoothxDeriv(m) @@ -551,6 +593,15 @@ class Simple(BaseRegularization): deriv += self._evalSmoothzDeriv(m) return deriv + @Utils.timeIt + def _evalSmooth2Deriv(self, m, v=None): + deriv = self._evalSmoothx2Deriv(m, v) + if self.regmesh.dim > 1: + deriv += self._evalSmoothy2Deriv(m, v) + if self.regmesh.dim > 2: + deriv += self._evalSmoothz2Deriv(m, v) + return deriv + @Utils.timeIt def eval(self, m): @@ -574,6 +625,10 @@ class Simple(BaseRegularization): """ return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + @Utils.timeIt + def eval2Deriv(self, m, v=None): + return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + class Tikhonov(Simple): @@ -740,7 +795,37 @@ class Tikhonov(Simple): return r.T * ( self.Wzz * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmooth2Deriv(self, m): + def _evalSmoothxx2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wxx * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wxx * self.mapping.deriv(m) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + + @Utils.timeIt + def _evalSmoothyy2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wyy * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wyy * self.mapping.deriv(m) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + + @Utils.timeIt + def _evalSmoothzz2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wzz * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wzz * self.mapping.deriv(m) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + + @Utils.timeIt + def _evalSmoothDeriv2(self, m): deriv = self._evalSmoothxxDeriv(m) if self.regmesh.dim > 1: deriv += self._evalSmoothyyDeriv(m) @@ -748,6 +833,15 @@ class Tikhonov(Simple): deriv += self._evalSmoothzzDeriv(m) return deriv + @Utils.timeIt + def _evalSmooth2Deriv2(self, m, v=None): + deriv = self._evalSmoothxx2Deriv(m, v) + if self.regmesh.dim > 1: + deriv += self._evalSmoothyy2Deriv(m, v) + if self.regmesh.dim > 2: + deriv += self._evalSmoothzz2Deriv(m, v) + return deriv + @Utils.timeIt def eval(self, m): @@ -769,7 +863,24 @@ class Tikhonov(Simple): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmooth2Deriv(m) + return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m) + + def eval2Deriv(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})} + + """ + return self._evalSmall2Deriv(m) + self._evalSmooth2Deriv(m) + self._evalSmooth2Deriv2(m) From 3dd9ecc9cd7f95b596d9027595013f9ca4647035 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 22:10:19 -0700 Subject: [PATCH 10/17] fix tikhonov 2Deriv --- SimPEG/Regularization.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 5947ed03..8f1b2d2c 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -865,7 +865,7 @@ class Tikhonov(Simple): """ return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m) - def eval2Deriv(self, m): + def eval2Deriv(self, m, v=None): """ The regularization is: @@ -880,7 +880,7 @@ class Tikhonov(Simple): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - return self._evalSmall2Deriv(m) + self._evalSmooth2Deriv(m) + self._evalSmooth2Deriv2(m) + return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + self._evalSmooth2Deriv2(m, v) From e10d6878fb86b129953e68c3749307321a7d8eea Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 11 May 2016 07:58:06 -0700 Subject: [PATCH 11/17] Remove Wsmooth from def W and replace by parts --- SimPEG/Regularization.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 8f1b2d2c..3263ba9a 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -469,7 +469,11 @@ class Simple(BaseRegularization): def W(self): """Full regularization matrix W""" if getattr(self, '_W', None) is None: - wlist = (self.Wsmall, self.Wsmooth) + wlist = (self.Wsmall, self.Wx) + if self.regmesh.dim > 1: + wlist += (self.Wy,) + if self.regmesh.dim > 2: + wlist += (self.Wz,) self._W = sp.vstack(wlist) return self._W From cd2360b8153a9b4f343cfdb81d084e5f36c49132 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 11 May 2016 23:04:14 -0700 Subject: [PATCH 12/17] Stash the regularization between each beta --- SimPEG/Directives.py | 24 ++++++++++++-- SimPEG/Regularization.py | 67 ++++++++++++++++++++++------------------ 2 files changed, 58 insertions(+), 33 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 008cc57b..a9b74276 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -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): """ diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 3263ba9a..5dc5aba0 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -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): From 3cc46131a3138096441fb4c423c5f4f21543d2b4 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Thu, 12 May 2016 08:31:01 -0700 Subject: [PATCH 13/17] Temporary change ... comment out W and Wsmooth --- SimPEG/Directives.py | 6 ++--- SimPEG/Regularization.py | 52 ++++++++++++++++++++-------------------- 2 files changed, 28 insertions(+), 30 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index a9b74276..6b15be9f 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -278,7 +278,6 @@ class Update_IRLS(InversionDirective): 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 @@ -308,20 +307,19 @@ class Update_IRLS(InversionDirective): 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 - + print "New gamma ", np.linalg.norm(self.reg.gamma) + # 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 class Update_lin_PreCond(InversionDirective): """ diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 5dc5aba0..422c8f9c 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -451,32 +451,32 @@ class Simple(BaseRegularization): self._Wz = Utils.sdiag((self.alpha_z * (self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.regmesh.cellDiffzStencil return self._Wz - @property - def Wsmooth(self): - """Full smoothness regularization matrix W""" - print 'wtf why are we using Wsmooth' - raise NotImplementedError - 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 self._Wsmooth - - @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: - wlist += (self.Wy,) - if self.regmesh.dim > 2: - wlist += (self.Wz,) - self._W = sp.vstack(wlist) - return self._W +# @property +# def Wsmooth(self): +# """Full smoothness regularization matrix W""" +# print 'wtf why are we using Wsmooth' +# raise NotImplementedError +# 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 self._Wsmooth +# +# @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: +# wlist += (self.Wy,) +# if self.regmesh.dim > 2: +# wlist += (self.Wz,) +# self._W = sp.vstack(wlist) +# return self._W @Utils.timeIt From fd3bde787f3a9bc570ad8386e36eb16a7d3da1d3 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Thu, 12 May 2016 14:58:16 -0700 Subject: [PATCH 14/17] Propose change to the Projected_GNCG solver. Add inner GN iterations. Nice improvement to the convergence of IRLS --- SimPEG/Directives.py | 49 ++++++++++++++++++---- SimPEG/Examples/Inversion_IRLS.py | 20 ++++----- SimPEG/Optimization.py | 70 +++++++++++++++++++++++++++++++ 3 files changed, 118 insertions(+), 21 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 6b15be9f..f5c25249 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -144,6 +144,35 @@ class BetaSchedule(InversionDirective): if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter self.invProb.beta /= self.coolingFactor +#class BetaSchedule_PGN_CG(InversionDirective): +# """BetaSchedule""" +# +# coolingFactor = 5. +# coolingRate = 1 +# GN_step_last = None +# GN_step_c = None +# +# def endIter(self): +# +# """ Compute the change in GN step, and proceed with cooling if below tol""" +# if self.opt.iter == 1: +# self.GN_step_last = np.linalg.norm(self.opt.xc - self.opt.x_last) +# d_GN_step = 1. +# +# else: +# self.GN_step_c = np.linalg.norm(self.opt.xc - self.opt.x_last) +# d_GN_step = self.GN_step_c / self.GN_step_last +# +# # Re-initiate last GN step +# self.GN_step_last = self.GN_step_c +# +# print "GN_step_last: ", self.GN_step_last +# print "d_GN_step: ", d_GN_step +# +# if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: +# if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter +# self.invProb.beta /= self.coolingFactor + class TargetMisfit(InversionDirective): @property @@ -265,13 +294,16 @@ class Update_IRLS(InversionDirective): if getattr(self, 'phi_m_last', None) is not None: self.reg.curModel = self.invProb.curModel + self.reg._Wsmall = None + self.reg._Wx = None + self.reg._Wy = None + self.reg._Wz = None self.reg.gamma = 1. 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 - 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 @@ -295,12 +327,6 @@ class Update_IRLS(InversionDirective): # Get phi_m at the end of current iteration self.phi_m_last = self.invProb.phi_m_last - # Update the model used for the IRLS weights - self.reg.curModel = self.invProb.curModel - - # 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 @@ -308,13 +334,18 @@ class Update_IRLS(InversionDirective): self.reg._Wy = None self.reg._Wz = None + # Update the model used for the IRLS weights + self.reg.curModel = self.invProb.curModel + + # Temporarely set gamma to 1. to get raw phi_m + self.reg.gamma = 1. + # 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 - print "New gamma ", np.linalg.norm(self.reg.gamma) - + # Reset the regularization matrices again for new gamma self.reg._Wsmall = None self.reg._Wx = None diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index c925cce2..6551bf21 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -18,6 +18,8 @@ def run(N=200, plotIt=True): mesh = Mesh.TensorMesh([N]) m0 = np.ones(mesh.nC) * 1e-4 + mref = np.zeros(mesh.nC) + nk = 10 jk = np.linspace(1.,nk,nk) p = -2. @@ -51,12 +53,13 @@ def run(N=200, plotIt=True): wr = ( wr/np.max(wr) ) reg = Regularization.Simple(mesh) - reg.wght = wr + reg.mref = mref + reg.cell_weights = wr dmis = DataMisfit.l2_DataMisfit(survey) dmis.Wd = 1./wd - opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4) + opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, maxIterGN=1, tolCG = 1e-4) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) invProb.curModel = m0 @@ -76,22 +79,15 @@ def run(N=200, plotIt=True): phid = invProb.phi_d reg = Regularization.Sparse(mesh) + reg.mref = mref + reg.cell_weights = wr -#============================================================================== -# fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2)) -# dmdx = reg.mesh.cellDiffxStencil * mrec -# plt.plot(np.sort(dmdx)) -#============================================================================== - - #reg.recModel = mrec - # reg.cell_weight = np.ones(mesh.nC) reg.mref = np.zeros(mesh.nC) reg.eps_p = 5e-2 reg.eps_q = 1e-2 reg.norms = [0., 0., 2., 2.] - reg.cell_weight = wr - opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3) + opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) #betaest = Directives.BetaEstimate_ByEig() diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 77704733..95f53320 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -893,6 +893,10 @@ class ProjectedGNCG(BFGS, Minimize, Remember): lower = -np.inf upper = np.inf + # Variables to control inner GN iterations + rdm_tol = 1e-2 # Tolerance for largest change in step (Default 1%) + maxIterGN = 3 # Maximum number of GN inner iterations + def _startup(self, x0): # ensure bound vectors are the same size as the model if type(self.lower) is not np.ndarray: @@ -938,6 +942,72 @@ class ProjectedGNCG(BFGS, Minimize, Remember): def approxHinv(self, value): self._approxHinv = value + @Utils.timeIt + def minimize(self, evalFunction, x0): + """minimize(evalFunction, x0) + + Minimizes the function (evalFunction) starting at the location x0. + + :param def evalFunction: function handle that evaluates: f, g, H = F(x) + :param numpy.ndarray x0: starting location + :rtype: numpy.ndarray + :return: x, the last iterate of the optimization algorithm + + The GN newton steps are repeated until it the maxIterGN is reached or the + relative step change falls below some tolrerance. + + """ + self.evalFunction = evalFunction + self.startup(x0) + self.printInit() + + + while True: + self.doStartIteration() + self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) + self.printIter() + if self.stoppingCriteria(): break + + # Inner GN iterations, stop on maximum number of iterations + # or on tolerance for step length change + GN_count = 0 + dm0 = None # Initial GN step length + dmc = None # Current GN step length + rdm = 1. # Relative change in step length + + while rdm > self.rdm_tol and GN_count < self.maxIterGN: + + GN_count += 1 + self.searchDirection = self.findSearchDirection() + p = self.scaleSearchDirection(self.searchDirection) + xt, passLS = self.modifySearchDirection(p) + if not passLS: + xt, caught = self.modifySearchDirectionBreak(p) + if not caught: return self.xc + + if GN_count == 1: + dm0 = np.linalg.norm(self.xc - xt) + dmc = dm0 + + else: + dmc = np.linalg.norm(self.xc - xt) + + rdm = dmc / dm0 + self.xc = xt # Update current model + + # Form system for next iteration + self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) + + print "GN iter: %i,\t dm: %8.5e,\t rdm: %8.5e"% (GN_count, dmc, rdm) + + self.doEndIteration(xt) + if self.stopNextIteration: break + + self.printDone() + self.finish() + + return self.xc + @Utils.timeIt def findSearchDirection(self): From 022e1f7660dd0b20e44b00ebb919a0a2cf77b036 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 27 May 2016 13:11:31 -0700 Subject: [PATCH 15/17] Update IRLS directive to allow multiple GN iterations. Remove modifications to the ProjGN solver. Update IRLS example. --- SimPEG/Directives.py | 156 ++++++++++++++++-------------- SimPEG/Examples/Inversion_IRLS.py | 14 +-- SimPEG/Optimization.py | 72 +------------- SimPEG/Regularization.py | 38 ++++++-- 4 files changed, 121 insertions(+), 159 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index f5c25249..b84d78d1 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -144,34 +144,6 @@ class BetaSchedule(InversionDirective): if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter self.invProb.beta /= self.coolingFactor -#class BetaSchedule_PGN_CG(InversionDirective): -# """BetaSchedule""" -# -# coolingFactor = 5. -# coolingRate = 1 -# GN_step_last = None -# GN_step_c = None -# -# def endIter(self): -# -# """ Compute the change in GN step, and proceed with cooling if below tol""" -# if self.opt.iter == 1: -# self.GN_step_last = np.linalg.norm(self.opt.xc - self.opt.x_last) -# d_GN_step = 1. -# -# else: -# self.GN_step_c = np.linalg.norm(self.opt.xc - self.opt.x_last) -# d_GN_step = self.GN_step_c / self.GN_step_last -# -# # Re-initiate last GN step -# self.GN_step_last = self.GN_step_c -# -# print "GN_step_last: ", self.GN_step_last -# print "d_GN_step: ", d_GN_step -# -# if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: -# if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter -# self.invProb.beta /= self.coolingFactor class TargetMisfit(InversionDirective): @@ -286,10 +258,16 @@ class Update_IRLS(InversionDirective): gamma = None phi_m_last = None phi_d_last = None - - + f_old = None + f_min_change = 1e-1 + coolingRate = 3 + maxIRLSiter = 10 + + def initialize(self): - + + self.IRLSiter = 0 + # Scale the regularization for changes in norm if getattr(self, 'phi_m_last', None) is not None: @@ -310,48 +288,75 @@ class Update_IRLS(InversionDirective): self.reg._Wx = None self.reg._Wy = None self.reg._Wz = None - + if getattr(self, 'phi_d_last', None) is None: self.phi_d_last = self.invProb.phi_d + if getattr(self, 'f_last', None) is None: + self.f_old = self.invProb.evalFunction(self.reg.curModel, return_g=False, return_H=False) + print self.f_old def endIter(self): - # Cool the threshold parameter if required - if getattr(self, 'factor', None) is not None: - eps = self.reg.eps / self.factor + + + # Only update after GN iterations + if self.opt.iter % self.coolingRate == 0: + + self.IRLSiter += 1 - if getattr(self, 'eps_min', None) is not None: - self.reg.eps = np.max([self.eps_min,eps]) - else: - self.reg.eps = eps - - # Get phi_m at the end of current iteration - self.phi_m_last = self.invProb.phi_m_last - - # 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 - - # Update the model used for the IRLS weights - self.reg.curModel = self.invProb.curModel - - # Temporarely set gamma to 1. to get raw phi_m - self.reg.gamma = 1. - - # 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 - - # 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.f_change = np.abs(self.f_old - self.opt.f_last) / self.f_old + + print "Function decrease" + str(self.f_change) + + # Check for maximum number of IRLS cycles + if self.IRLSiter == self.maxIRLSiter: + self.opt.stopNextIteration = True + return + + # Check if the function has changed enough + if self.f_change < self.f_min_change: + self.opt.stopNextIteration = True + return + else: + self.f_old = self.opt.f_last + + # Cool the threshold parameter if required + if getattr(self, 'factor', None) is not None: + eps = self.reg.eps / self.factor + + if getattr(self, 'eps_min', None) is not None: + self.reg.eps = np.max([self.eps_min,eps]) + else: + self.reg.eps = eps + + # Get phi_m at the end of current iteration + self.phi_m_last = self.invProb.phi_m_last + + # 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 + + # Update the model used for the IRLS weights + self.reg.curModel = self.invProb.curModel + + # Temporarely set gamma to 1. to get raw phi_m + self.reg.gamma = 1. + + # 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 + + # Reset the regularization matrices again for new gamma + self.reg._Wsmall = None + self.reg._Wx = None + self.reg._Wy = None + self.reg._Wz = None + # Compute the change in class Update_lin_PreCond(InversionDirective): """ Create a Jacobi preconditioner for the linear problem @@ -411,11 +416,14 @@ class Scale_Beta(InversionDirective): update is done only if the misfit is outside some threshold bounds. """ tol = 0.05 - + coolingRate=5 + def endIter(self): - - # Check if misfit is within the tolerance, otherwise adjust beta - val = self.invProb.phi_d / (self.survey.nD*0.5) - - if np.abs(1.-val) > self.tol: - self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d + + # Only update after GN iterations + if self.opt.iter % self.coolingRate == 0: + # Check if misfit is within the tolerance, otherwise adjust beta + val = self.invProb.phi_d / (self.survey.nD*0.5) + + if np.abs(1.-val) > self.tol: + self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index 6551bf21..d1da048f 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -20,7 +20,7 @@ def run(N=200, plotIt=True): m0 = np.ones(mesh.nC) * 1e-4 mref = np.zeros(mesh.nC) - nk = 10 + nk = 15 jk = np.linspace(1.,nk,nk) p = -2. q = 1. @@ -59,7 +59,7 @@ def run(N=200, plotIt=True): dmis = DataMisfit.l2_DataMisfit(survey) dmis.Wd = 1./wd - opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, maxIterGN=1, tolCG = 1e-4) + opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) invProb.curModel = m0 @@ -83,18 +83,18 @@ def run(N=200, plotIt=True): reg.cell_weights = wr reg.mref = np.zeros(mesh.nC) - reg.eps_p = 5e-2 + reg.eps_p = 1e-3 reg.eps_q = 1e-2 reg.norms = [0., 0., 2., 2.] - opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) + opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) - #betaest = Directives.BetaEstimate_ByEig() + update_beta = Directives.Scale_Beta(tol = 0.05, coolingRate=5) target = Directives.TargetMisfit() - IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid ) + IRLS = Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid, coolingRate=5 ) - inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS]) + inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS,update_beta]) m0 = mrec diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 95f53320..05214dd5 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -893,10 +893,6 @@ class ProjectedGNCG(BFGS, Minimize, Remember): lower = -np.inf upper = np.inf - # Variables to control inner GN iterations - rdm_tol = 1e-2 # Tolerance for largest change in step (Default 1%) - maxIterGN = 3 # Maximum number of GN inner iterations - def _startup(self, x0): # ensure bound vectors are the same size as the model if type(self.lower) is not np.ndarray: @@ -942,72 +938,6 @@ class ProjectedGNCG(BFGS, Minimize, Remember): def approxHinv(self, value): self._approxHinv = value - @Utils.timeIt - def minimize(self, evalFunction, x0): - """minimize(evalFunction, x0) - - Minimizes the function (evalFunction) starting at the location x0. - - :param def evalFunction: function handle that evaluates: f, g, H = F(x) - :param numpy.ndarray x0: starting location - :rtype: numpy.ndarray - :return: x, the last iterate of the optimization algorithm - - The GN newton steps are repeated until it the maxIterGN is reached or the - relative step change falls below some tolrerance. - - """ - self.evalFunction = evalFunction - self.startup(x0) - self.printInit() - - - while True: - self.doStartIteration() - self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) - self.printIter() - if self.stoppingCriteria(): break - - # Inner GN iterations, stop on maximum number of iterations - # or on tolerance for step length change - GN_count = 0 - dm0 = None # Initial GN step length - dmc = None # Current GN step length - rdm = 1. # Relative change in step length - - while rdm > self.rdm_tol and GN_count < self.maxIterGN: - - GN_count += 1 - self.searchDirection = self.findSearchDirection() - p = self.scaleSearchDirection(self.searchDirection) - xt, passLS = self.modifySearchDirection(p) - if not passLS: - xt, caught = self.modifySearchDirectionBreak(p) - if not caught: return self.xc - - if GN_count == 1: - dm0 = np.linalg.norm(self.xc - xt) - dmc = dm0 - - else: - dmc = np.linalg.norm(self.xc - xt) - - rdm = dmc / dm0 - self.xc = xt # Update current model - - # Form system for next iteration - self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) - - print "GN iter: %i,\t dm: %8.5e,\t rdm: %8.5e"% (GN_count, dmc, rdm) - - self.doEndIteration(xt) - if self.stopNextIteration: break - - self.printDone() - self.finish() - - return self.xc - @Utils.timeIt def findSearchDirection(self): @@ -1078,4 +1008,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember): indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0)) delx[indx] = 0. - return delx + return delx \ No newline at end of file diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 422c8f9c..c3b1c9e5 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -890,14 +890,37 @@ class Tikhonov(Simple): class Sparse(Simple): - + """ + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R W(m-m_\\text{ref})} + + where the IRLS weight + + .. math:: + + R = \eta TO FINISH LATER!!! + + So the derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})} + + The IRLS weights are recomputed after each beta solves. + It is strongly recommended to do a few Gauss-Newton iterations + before updating. + """ + # set default values - eps_p = 1e-1 - eps_q = 1e-1 - curModel = None # use a model to compute the weights - gamma = 1. - norms = [0., 2., 2., 2.] - cell_weights = 1. + eps_p = 1e-1 # Threshold value for the model norm + eps_q = 1e-1 # Threshold value for the model gradient norm + curModel = None # Requires model to compute the weights + gamma = 1. # Model norm scaling to smooth out convergence + norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz) + cell_weights = 1. # Consider overwriting with sensitivity weights def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) @@ -971,6 +994,7 @@ class Sparse(Simple): def R(self, f_m , eps, exponent): + # Eta scaling is important for mix-norms...do not mess with it eta = (eps**(1.-exponent/2.))**0.5 r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.) From 3b4bec9c0b40b0ce949cea8422adef1f29b81bb7 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Sat, 28 May 2016 11:27:09 -0700 Subject: [PATCH 16/17] Refactor IRLS iterations, full solves from l2->lp Adapt Example --- SimPEG/Directives.py | 206 ++++++++++++++++++------------ SimPEG/Examples/Inversion_IRLS.py | 66 +++++----- SimPEG/Regularization.py | 1 + 3 files changed, 158 insertions(+), 115 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index b84d78d1..8727d925 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -144,7 +144,7 @@ class BetaSchedule(InversionDirective): if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter self.invProb.beta /= self.coolingFactor - + class TargetMisfit(InversionDirective): @property @@ -238,12 +238,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): # Save the file as a npz np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred) - -# class UpdateReferenceModel(Parameter): - -# mref0 = None - -# def nextIter(self): # mref = getattr(self, 'm_prev', None) # if mref is None: # if self.debug: print 'UpdateReferenceModel is using mref0' @@ -254,109 +248,165 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): class Update_IRLS(InversionDirective): eps_min = None + eps_p = None + eps_q = None + norms = [2.,2.,2.,2.] factor = None gamma = None phi_m_last = None phi_d_last = None f_old = None - f_min_change = 1e-1 - coolingRate = 3 - maxIRLSiter = 10 - + f_min_change = 1e-2 + beta_tol = 5e-2 + # Solving parameter for IRLS (mode:2) + IRLSiter = 0 + minGNiter = 6 + maxIRLSiter = 10 + iterStart = 0 + + # Beta schedule + coolingFactor = 2. + coolingRate = 1 + + mode = 1 + + @property + def target(self): + if getattr(self, '_target', None) is None: + self._target = self.survey.nD + return self._target + @target.setter + def target(self, val): + self._target = val + def initialize(self): - - self.IRLSiter = 0 - - # Scale the regularization for changes in norm - if getattr(self, 'phi_m_last', None) is not None: - self.reg.curModel = self.invProb.curModel - self.reg._Wsmall = None - self.reg._Wx = None - self.reg._Wy = None - self.reg._Wz = None - self.reg.gamma = 1. - phim_new = self.reg.eval(self.invProb.curModel) - self.gamma = self.phi_m_last / phim_new + if self.mode == 1: + self.reg.norms = [2., 2., 2., 2.] +# # Scale the regularization for changes in norm +# if getattr(self, 'phi_m_last', None) is not None: +# +# self.reg.curModel = self.invProb.curModel +# self.reg._Wsmall = None +# self.reg._Wx = None +# self.reg._Wy = None +# self.reg._Wz = None +# self.reg.gamma = 1. +# 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 +# # 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.curModel = self.invProb.curModel - self.reg.gamma = self.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 - - if getattr(self, 'phi_d_last', None) is None: - self.phi_d_last = self.invProb.phi_d - - if getattr(self, 'f_last', None) is None: - self.f_old = self.invProb.evalFunction(self.reg.curModel, return_g=False, return_H=False) - print self.f_old - def endIter(self): - - - # Only update after GN iterations - if self.opt.iter % self.coolingRate == 0: +# if getattr(self, 'phi_d_last', None) is None: +# self.phi_d_last = self.invProb.phi_d +# +# if getattr(self, 'f_last', None) is None: +# self.f_old = self.invProb.evalFunction(self.reg.curModel, return_g=False, return_H=False) +# print self.f_old + def endIter(self): + + # After reaching target misfit with l2-norm, switch to IRLS (mode:2) + if self.invProb.phi_d < self.target and self.mode == 1: + print "Convergence with smooth l2-norm regularization: Start IRLS steps..." + + self.mode = 2 + print self.eps_p, self.eps_q, self.norms + self.reg.eps_p = self.eps_p + self.reg.eps_q = self.eps_q + self.reg.norms = self.norms + self.coolingFactor = 1. + self.coolingRate = 1 + self.iterStart = self.opt.iter + self.phi_d_last = self.invProb.phi_d + self.phi_m_last = self.invProb.phi_m_last + + self.reg.l2model = self.invProb.curModel + self.reg.curModel = self.invProb.curModel + + if getattr(self, 'f_old', None) is None: + self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False) + + if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: + if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter + self.invProb.beta /= self.coolingFactor + + + # Only update after GN iterations + if (self.opt.iter-self.iterStart) % self.minGNiter == 0 and self.mode==2: + self.IRLSiter += 1 - self.f_change = np.abs(self.f_old - self.opt.f_last) / self.f_old - - print "Function decrease" + str(self.f_change) - + phim_new = self.reg.eval(self.invProb.curModel) + self.f_change = np.abs(self.f_old - phim_new) / self.f_old + + print "Regularization decrease: %6.3e" % (self.f_change) + # Check for maximum number of IRLS cycles if self.IRLSiter == self.maxIRLSiter: + print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter self.opt.stopNextIteration = True return - + # Check if the function has changed enough - if self.f_change < self.f_min_change: + if self.f_change < self.f_min_change and self.IRLSiter > 1: + print "Minimum decrease in regularization. End of IRLS" self.opt.stopNextIteration = True - return - else: - self.f_old = self.opt.f_last - + return + else: + self.f_old = phim_new + # Cool the threshold parameter if required if getattr(self, 'factor', None) is not None: eps = self.reg.eps / self.factor - + if getattr(self, 'eps_min', None) is not None: self.reg.eps = np.max([self.eps_min,eps]) else: self.reg.eps = eps - + # Get phi_m at the end of current iteration self.phi_m_last = self.invProb.phi_m_last - + # 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 - + # Update the model used for the IRLS weights self.reg.curModel = self.invProb.curModel - + # Temporarely set gamma to 1. to get raw phi_m self.reg.gamma = 1. - + # 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 - + # Reset the regularization matrices again for new gamma self.reg._Wsmall = None self.reg._Wx = None self.reg._Wy = None self.reg._Wz = None - # Compute the change in + # Check if misfit is within the tolerance, otherwise adjust beta + val = self.invProb.phi_d / (self.survey.nD*0.5) + + if np.abs(1.-val) > self.beta_tol: + self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d + class Update_lin_PreCond(InversionDirective): """ Create a Jacobi preconditioner for the linear problem @@ -409,21 +459,15 @@ class Update_Wj(InversionDirective): self.reg.wght = JtJdiag -class Scale_Beta(InversionDirective): - """ - Instead of a linear cooling schedule, beta is allowed to change based - on the ratio between the target misfit and the current data misfit. The - update is done only if the misfit is outside some threshold bounds. - """ - tol = 0.05 - coolingRate=5 - - def endIter(self): - - # Only update after GN iterations - if self.opt.iter % self.coolingRate == 0: - # Check if misfit is within the tolerance, otherwise adjust beta - val = self.invProb.phi_d / (self.survey.nD*0.5) - - if np.abs(1.-val) > self.tol: - self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d +#class Scale_Beta(InversionDirective): +# """ +# Instead of a linear cooling schedule, beta is allowed to change based +# on the ratio between the target misfit and the current data misfit. The +# update is done only if the misfit is outside some threshold bounds. +# """ +# tol = 0.05 +# coolingRate=5 +# +# def endIter(self): +# +# diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index d1da048f..5c1a8770 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -52,51 +52,49 @@ def run(N=200, plotIt=True): wr = np.sum(prob.G**2.,axis=0)**0.5 wr = ( wr/np.max(wr) ) - reg = Regularization.Simple(mesh) - reg.mref = mref - reg.cell_weights = wr - +# reg = Regularization.Simple(mesh) +# reg.mref = mref +# reg.cell_weights = wr +# dmis = DataMisfit.l2_DataMisfit(survey) dmis.Wd = 1./wd - - opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4) - invProb = InvProblem.BaseInvProblem(dmis, reg, opt) - invProb.curModel = m0 - - beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1) - target = Directives.TargetMisfit() - +# +# opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4) +# invProb = InvProblem.BaseInvProblem(dmis, reg, opt) +# invProb.curModel = m0 +# +# beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1) +# target = Directives.TargetMisfit() +# betaest = Directives.BetaEstimate_ByEig() - inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target]) - - - mrec = inv.run(m0) - ml2 = mrec - print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) - - # Switch regularization to sparse - phim = invProb.phi_m_last - phid = invProb.phi_d +# inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target]) +# +# +# mrec = inv.run(m0) +# ml2 = mrec +# print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) +# +# # Switch regularization to sparse +# phim = invProb.phi_m_last +# phid = invProb.phi_d reg = Regularization.Sparse(mesh) reg.mref = mref reg.cell_weights = wr reg.mref = np.zeros(mesh.nC) - reg.eps_p = 1e-3 - reg.eps_q = 1e-2 - reg.norms = [0., 0., 2., 2.] + eps_p = 1e-3 + eps_q = 1e-2 + norms = [0., 0., 2., 2.] opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) - invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) - beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) - update_beta = Directives.Scale_Beta(tol = 0.05, coolingRate=5) - target = Directives.TargetMisfit() - IRLS = Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid, coolingRate=5 ) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt) + #beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) + #update_beta = Directives.Scale_Beta(tol = 0.05, coolingRate=5) +# target = Directives.TargetMisfit() + IRLS = Directives.Update_IRLS( norms=norms, eps_p=eps_p, eps_q=eps_q) - inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS,update_beta]) - - m0 = mrec + inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest]) # Run inversion mrec = inv.run(m0) @@ -113,7 +111,7 @@ def run(N=200, plotIt=True): axes[0].set_title('Columns of matrix G') axes[1].plot(mesh.vectorCCx, mtrue, 'b-') - axes[1].plot(mesh.vectorCCx, ml2, 'r-') + axes[1].plot(mesh.vectorCCx, reg.l2model, 'r-') #axes[1].legend(('True Model', 'Recovered Model')) axes[1].set_ylim(-1.0,1.25) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index c3b1c9e5..3304022a 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -918,6 +918,7 @@ class Sparse(Simple): eps_p = 1e-1 # Threshold value for the model norm eps_q = 1e-1 # Threshold value for the model gradient norm curModel = None # Requires model to compute the weights + l2model = None gamma = 1. # Model norm scaling to smooth out convergence norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz) cell_weights = 1. # Consider overwriting with sensitivity weights From e476bf0059b3c6f8e75665bffa1e699b586919a1 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Sun, 29 May 2016 14:09:29 -0700 Subject: [PATCH 17/17] Cleanup Sparse Reg and Directives --- SimPEG/Directives.py | 71 +++++++------------------------ SimPEG/Examples/Inversion_IRLS.py | 16 +++---- 2 files changed, 23 insertions(+), 64 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 8727d925..98770666 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -258,66 +258,39 @@ class Update_IRLS(InversionDirective): f_old = None f_min_change = 1e-2 beta_tol = 5e-2 - + # Solving parameter for IRLS (mode:2) IRLSiter = 0 - minGNiter = 6 + minGNiter = 5 maxIRLSiter = 10 iterStart = 0 - + # Beta schedule coolingFactor = 2. coolingRate = 1 - + mode = 1 - + @property def target(self): if getattr(self, '_target', None) is None: - self._target = self.survey.nD + self._target = self.survey.nD*0.5 return self._target @target.setter def target(self, val): self._target = val - + def initialize(self): if self.mode == 1: self.reg.norms = [2., 2., 2., 2.] -# # Scale the regularization for changes in norm -# if getattr(self, 'phi_m_last', None) is not None: -# -# self.reg.curModel = self.invProb.curModel -# self.reg._Wsmall = None -# self.reg._Wx = None -# self.reg._Wy = None -# self.reg._Wz = None -# self.reg.gamma = 1. -# 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 -# # 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 -# if getattr(self, 'phi_d_last', None) is None: -# self.phi_d_last = self.invProb.phi_d -# -# if getattr(self, 'f_last', None) is None: -# self.f_old = self.invProb.evalFunction(self.reg.curModel, return_g=False, return_H=False) -# print self.f_old - def endIter(self): # After reaching target misfit with l2-norm, switch to IRLS (mode:2) if self.invProb.phi_d < self.target and self.mode == 1: print "Convergence with smooth l2-norm regularization: Start IRLS steps..." - + self.mode = 2 print self.eps_p, self.eps_q, self.norms self.reg.eps_p = self.eps_p @@ -328,17 +301,18 @@ class Update_IRLS(InversionDirective): self.iterStart = self.opt.iter self.phi_d_last = self.invProb.phi_d self.phi_m_last = self.invProb.phi_m_last - + self.reg.l2model = self.invProb.curModel self.reg.curModel = self.invProb.curModel - + if getattr(self, 'f_old', None) is None: self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False) - + + # Beta Schedule if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter self.invProb.beta /= self.coolingFactor - + # Only update after GN iterations if (self.opt.iter-self.iterStart) % self.minGNiter == 0 and self.mode==2: @@ -401,12 +375,12 @@ class Update_IRLS(InversionDirective): self.reg._Wy = None self.reg._Wz = None - # Check if misfit is within the tolerance, otherwise adjust beta + # Check if misfit is within the tolerance, otherwise scale beta val = self.invProb.phi_d / (self.survey.nD*0.5) - + if np.abs(1.-val) > self.beta_tol: self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d - + class Update_lin_PreCond(InversionDirective): """ Create a Jacobi preconditioner for the linear problem @@ -458,16 +432,3 @@ class Update_Wj(InversionDirective): JtJdiag = JtJdiag / max(JtJdiag) self.reg.wght = JtJdiag - -#class Scale_Beta(InversionDirective): -# """ -# Instead of a linear cooling schedule, beta is allowed to change based -# on the ratio between the target misfit and the current data misfit. The -# update is done only if the misfit is outside some threshold bounds. -# """ -# tol = 0.05 -# coolingRate=5 -# -# def endIter(self): -# -# diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index 5c1a8770..afd90525 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -1,7 +1,7 @@ from SimPEG import * -def run(N=200, plotIt=True): +def run(N=100, plotIt=True): """ Inversion: Linear Problem ========================= @@ -19,8 +19,8 @@ def run(N=200, plotIt=True): m0 = np.ones(mesh.nC) * 1e-4 mref = np.zeros(mesh.nC) - - nk = 15 + + nk = 10 jk = np.linspace(1.,nk,nk) p = -2. q = 1. @@ -83,18 +83,16 @@ def run(N=200, plotIt=True): reg.cell_weights = wr reg.mref = np.zeros(mesh.nC) - eps_p = 1e-3 - eps_q = 1e-2 + eps_p = 5e-2 + eps_q = 5e-2 norms = [0., 0., 2., 2.] opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) - #beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) - #update_beta = Directives.Scale_Beta(tol = 0.05, coolingRate=5) -# target = Directives.TargetMisfit() + update_Jacobi = Directives.Update_lin_PreCond() IRLS = Directives.Update_IRLS( norms=norms, eps_p=eps_p, eps_q=eps_q) - inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest]) + inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi]) # Run inversion mrec = inv.run(m0)