mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-05 10:17:03 +08:00
Refactor IRLS iterations, full solves from l2->lp
Adapt Example
This commit is contained in:
D Fournier
+125
-81
@@ -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):
|
||||
#
|
||||
#
|
||||
|
||||
@@ -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)
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
Reference in New Issue
Block a user