# Inversion: Linear: IRLSΒΆ

Here we go over the basics of creating a linear problem and inversion.

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 from __future__ import print_function import numpy as np import matplotlib.pyplot as plt from SimPEG import Mesh from SimPEG import Problem from SimPEG import Survey from SimPEG import DataMisfit from SimPEG import Directives from SimPEG import Optimization from SimPEG import Regularization from SimPEG import InvProblem from SimPEG import Inversion def run(N=100, plotIt=True): """ Inversion: Linear: IRLS ======================= Here we go over the basics of creating a linear problem and inversion. """ np.random.seed(1) std_noise = 1e-2 mesh = Mesh.TensorMesh([N]) m0 = np.ones(mesh.nC) * 1e-4 mref = np.zeros(mesh.nC) nk = 20 jk = np.linspace(1., 60., nk) p = -0.25 q = 0.25 def g(k): return ( np.exp(p*jk[k]*mesh.vectorCCx) * np.cos(np.pi*q*jk[k]*mesh.vectorCCx) ) G = np.empty((nk, mesh.nC)) for i in range(nk): G[i, :] = g(i) mtrue = np.zeros(mesh.nC) mtrue[mesh.vectorCCx > 0.3] = 1. mtrue[mesh.vectorCCx > 0.45] = -0.5 mtrue[mesh.vectorCCx > 0.6] = 0 prob = Problem.LinearProblem(mesh, G=G) survey = Survey.LinearSurvey() survey.pair(prob) survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk) wd = np.ones(nk) * std_noise # Distance weighting wr = np.sum(prob.G**2., axis=0)**0.5 wr = wr/np.max(wr) dmis = DataMisfit.l2_DataMisfit(survey) dmis.Wd = 1./wd betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e-2) reg = Regularization.Sparse(mesh) reg.mref = mref reg.cell_weights = wr reg.mref = np.zeros(mesh.nC) opt = Optimization.ProjectedGNCG( maxIter=100, lower=-2., upper=2., maxIterLS=20, maxIterCG=10, tolCG=1e-3 ) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) update_Jacobi = Directives.Update_lin_PreCond() # Set the IRLS directive, penalize the lowest 25 percentile of model values # Start with an l2-l2, then switch to lp-norms norms = [0., 0., 2., 2.] IRLS = Directives.Update_IRLS( norms=norms, prctile=25, maxIRLSiter=15, minGNiter=3 ) inv = Inversion.BaseInversion( invProb, directiveList=[IRLS, betaest, update_Jacobi] ) # Run inversion mrec = inv.run(m0) print("Final misfit:" + str(invProb.dmisfit.eval(mrec))) if plotIt: fig, axes = plt.subplots(1, 2, figsize=(12*1.2, 4*1.2)) 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(mesh.vectorCCx, mtrue, 'b-') axes[1].plot(mesh.vectorCCx, reg.l2model, 'r-') # axes[1].legend(('True Model', 'Recovered Model')) axes[1].set_ylim(-1.0, 1.25) axes[1].plot(mesh.vectorCCx, mrec, 'k-', lw=2) axes[1].legend( ( 'True Model', 'Smooth l2-l2', 'Sparse lp: {0}, lqx: {1}'.format(*reg.norms) ), fontsize=12 ) return prob, survey, mesh, mrec if __name__ == '__main__': run() plt.show()