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):