From 8f73b2e7be524a15411b4ab6765fae8b9fc5b81d Mon Sep 17 00:00:00 2001 From: D Fournier Date: Sun, 3 Apr 2016 17:18:39 -0700 Subject: [PATCH] Update directives Add IRLS example --- SimPEG/Directives.py | 101 +++++++++++++--------- SimPEG/Examples/Inversion_IRLS.py | 134 ++++++++++++++++++++---------- SimPEG/Regularization.py | 22 +++-- SimPEG/Survey.py | 2 +- 4 files changed, 160 insertions(+), 99 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index e5a63547..de7e0b6b 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -237,39 +237,41 @@ 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 SaveOutputDictEveryIteration(_SaveEveryIteration): - """SaveOutputDictEveryIteration - A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info). - """ - - def initialize(self): - print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName - - def endIter(self): - # Save the data. - ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) ) - phi_ms = 0.5*ms.dot(ms) - if self.reg.mrefInSmooth == True: - mref = self.reg.mref - else: - mref = 0 - mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) ) - phi_mx = 0.5 * mx.dot(mx) - if self.prob.mesh.dim==2: - my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) ) - phi_my = 0.5 * my.dot(my) - else: - phi_my = 'NaN' - if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType: - mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) ) - phi_mz = 0.5 * mz.dot(mz) - else: - phi_mz = 'NaN' - - - # Save the file as a npz - np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), 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 SaveOutputDictEveryIteration(_SaveEveryIteration): +# """SaveOutputDictEveryIteration +# A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info). +# """ +# +# def initialize(self): +# print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName +# +# def endIter(self): +# # Save the data. +# ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) ) +# phi_ms = 0.5*ms.dot(ms) +# if self.reg.mrefInSmooth == True: +# mref = self.reg.mref +# else: +# mref = 0 +# mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) ) +# phi_mx = 0.5 * mx.dot(mx) +# if self.prob.mesh.dim==2: +# my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) ) +# phi_my = 0.5 * my.dot(my) +# else: +# phi_my = 'NaN' +# if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType: +# mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) ) +# phi_mz = 0.5 * mz.dot(mz) +# else: +# phi_mz = 'NaN' +# +# +# # Save the file as a npz +# np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), 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): @@ -295,6 +297,8 @@ class update_IRLS(InversionDirective): # 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.gamma = 1. phim_new = self.reg.eval(self.invProb.curModel) self.gamma = self.phi_m_last / phim_new @@ -320,22 +324,39 @@ class update_IRLS(InversionDirective): # Update the model used for the IRLS weights self.reg.curModel = self.invProb.curModel - - # Update the pre-conditioner - diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.curModel.size))**2. - PC = Utils.sdiag(diagA**-1.) - self.opt.approxHinv = PC # Temporarely set gamma to 1. self.reg.gamma = 1. # Compute change in model objective function and update scaling phim_new = self.reg.eval(self.invProb.curModel) + self.reg.gamma = self.phi_m_last / phim_new - # TO DO: Re-scale beta if too much change in misfit - self.invProb.beta = self.invProb.beta * self.phi_d_last / self.invProb.phi_d + # TO DO: Check optimization class, data misfit not matching reality + #dpred = self.prob.fields(self.invProb.curModel) + #phid = self.invProb.dmisfit.eval(self.invProb.curModel) + #print self.survey.std[0] + #print phid + #print self.invProb.phi_d + #print self.invProb.phi_d_last + + self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d + +class update_lin_PreCond(InversionDirective): + + + def endIter(self): + # Cool the threshold parameter + + 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. + PC = Utils.sdiag(diagA**-1.) + self.opt.approxHinv = PC + print 'Updated pre-cond' + #============================================================================== # import pylab as plt # plt.figure() diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index a7f5bd54..53240139 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -1,7 +1,7 @@ from SimPEG import * -def run(N=100, plotIt=True): +def run(N=200, plotIt=True): """ Inversion: Linear Problem ========================= @@ -9,39 +9,21 @@ def run(N=100, plotIt=True): Here we go over the basics of creating a linear problem and inversion. """ - - class LinearSurvey(Survey.BaseSurvey): - def projectFields(self, u): - return u - - class LinearProblem(Problem.BaseProblem): - - surveyPair = LinearSurvey - - def __init__(self, mesh, G, **kwargs): - Problem.BaseProblem.__init__(self, mesh, **kwargs) - self.G = G - - def fields(self, m, u=None): - return self.G.dot(m) - - def Jvec(self, m, v, u=None): - return self.G.dot(v) - - def Jtvec(self, m, v, u=None): - return self.G.T.dot(v) - - + + np.random.seed(1) + std_noise = 1e-2 + mesh = Mesh.TensorMesh([N]) + + m0 = np.ones(mesh.nC) * 1e-4 + nk = 10 + jk = np.linspace(1.,nk,nk) + p = -2. + q = 1. - nk = 20 - jk = np.linspace(1.,20.,nk) - p = -0.25 - q = 0.25 - - g = lambda k: np.exp(p*jk[k]*mesh.vectorCCx)*np.cos(2*np.pi*q*jk[k]*mesh.vectorCCx) + g = lambda k: np.exp(p*jk[k]*mesh.vectorCCx)*np.cos(np.pi*q*jk[k]*mesh.vectorCCx) G = np.empty((nk, mesh.nC)) @@ -52,25 +34,43 @@ def run(N=100, plotIt=True): mtrue[mesh.vectorCCx > 0.3] = 1. mtrue[mesh.vectorCCx > 0.45] = -0.5 mtrue[mesh.vectorCCx > 0.6] = 0 + - prob = LinearProblem(mesh, G) - survey = LinearSurvey() + prob = Problem.LinearProblem(mesh, G) + survey = Survey.LinearSurvey() survey.pair(prob) - survey.makeSyntheticData(mtrue, std=0.01) - - M = prob.mesh - - reg = Regularization.Tikhonov(mesh) + survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk) + #survey.makeSyntheticData(mtrue, std=std_noise) + + wd = np.ones(nk) * std_noise + + #print survey.std[0] + #M = prob.mesh + # Distance weighting + wr = np.sum(prob.G**2.,axis=0)**0.5 + wr = ( wr/np.max(wr) )**0 + + reg = Regularization.Simple(mesh) + reg.wght = wr + dmis = DataMisfit.l2_DataMisfit(survey) - opt = Optimization.ProjectedGNCG(maxIter=20,lower=-1.,upper=1., maxIterCG= 20, tolCG = 1e-3) + dmis.Wd = 1./wd + + opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) - beta = Directives.BetaSchedule() + invProb.curModel = m0 + + beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1) + target = Directives.TargetMisfit() + betaest = Directives.BetaEstimate_ByEig() - inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest]) - m0 = np.zeros_like(survey.mtrue) + inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target]) + mrec = inv.run(m0) - + + print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) + if plotIt: import matplotlib.pyplot as plt @@ -78,12 +78,54 @@ def run(N=100, plotIt=True): for i in range(prob.G.shape[0]): axes[0].plot(prob.G[i,:]) axes[0].set_title('Columns of matrix G') - - axes[1].plot(M.vectorCCx, survey.mtrue, 'b-') - axes[1].plot(M.vectorCCx, mrec, 'r-') - axes[1].legend(('True Model', 'Recovered Model')) + + axes[1].plot(mesh.vectorCCx, mtrue, 'b-') + axes[1].plot(mesh.vectorCCx, mrec, 'r-') + #axes[1].legend(('True Model', 'Recovered Model')) + axes[1].set_ylim(-1.0,1.25) plt.show() + + # Switch regularization to sparse + phim = invProb.phi_m_last + phid = invProb.phi_d + + reg = Regularization.Sparse(mesh) + +#============================================================================== +# 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.wght = np.ones(mesh.nC) + reg.mref = np.zeros(mesh.nC) + reg.eps_p = 2e-3 + reg.eps_q = 2e-3 + reg.norms = [0., 0., 2., 2.] + reg.wght = wr + + opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterCG= 200, tolCG = 1e-3) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) + beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) + #betaest = Directives.BetaEstimate_ByEig() + target = Directives.TargetMisfit() + IRLS =Directives.update_IRLS( phi_m_last = phim, phi_d_last = phid ) + + inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS]) + m0 = mrec + + # Run inversion + mrec = inv.run(m0) + + print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) + + if plotIt: + axes[1].plot(mesh.vectorCCx, mrec, 'k-',lw = 2) + axes[1].legend(('True Model', 'Smooth l2-l2', + 'Sparse lp:' + str(reg.norms[0]) + ', lqx:' + str(reg.norms[1]) ), fontsize = 12) + return prob, survey, mesh, mrec if __name__ == '__main__': diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index ed75ada6..195e6203 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -643,13 +643,11 @@ class Simple(Tikhonov): class Sparse(Simple): # set default values - eps = 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., 1.] - qx = 2. - qy = 2. - qz = 2. + norms = [0., 2., 2., 2.] wght = 1. def __init__(self, mesh, mapping=None, indActive=None, **kwargs): @@ -666,7 +664,7 @@ class Sparse(Simple): else: f_m = self.curModel - self.reg.mref - self.rs = self.R(f_m , self.norms[0]) + 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 ) @@ -682,7 +680,7 @@ class Sparse(Simple): else: f_m = self.regmesh.cellDiffxStencil * self.curModel - self.rx = self.R( f_m , self.qx) + 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 @@ -696,7 +694,7 @@ class Sparse(Simple): else: f_m = self.regmesh.cellDiffyStencil * self.curModel - self.ry = self.R( f_m , self.qy) + 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 @@ -710,7 +708,7 @@ class Sparse(Simple): else: f_m = self.regmesh.cellDiffzStencil * self.curModel - self.rz = self.R( f_m , self.qz) + 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 @@ -735,9 +733,9 @@ class Sparse(Simple): #self._W = sp.vstack(wlist) return sp.vstack(wlist) - def R(self, f_m , exponent): + def R(self, f_m , eps, exponent): - eta = (self.eps**(1-exponent/2.))**0.5 - r = eta / (f_m**2.+self.eps**2.)**((1-exponent/2.)/2.) + eta = (eps**(1-exponent/2.))**0.5 + r = eta / (f_m**2.+ eps**2.)**((1-exponent/2.)/2.) return r diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 9f307c3f..eee2d538 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -372,7 +372,7 @@ class BaseSurvey(object): self.dtrue = self.dpred(m, u=u) noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape) self.dobs = self.dtrue+noise - self.std = self.dobs*0 + std + self.std = self.dobs*0. + std return self.dobs class LinearSurvey(BaseSurvey):