From 3b4bec9c0b40b0ce949cea8422adef1f29b81bb7 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Sat, 28 May 2016 11:27:09 -0700 Subject: [PATCH] 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