Inversion: Linear: IRLSΒΆ

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

../../../_images/sphx_glr_plot_inversion_linear_irls_001.png

Out:

SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||

    SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
    ***Done using same Solver and solverOpts as the problem***
model has any nan: 0
=============================== Projected GNCG ===============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
   0  1.16e+04  5.36e+05  1.36e-07  5.36e+05    2.00e+01      0
   1  5.81e+03  2.97e+03  1.20e+00  9.92e+03    1.81e+01      0
   2  2.90e+03  1.22e+03  1.40e+00  5.29e+03    1.81e+01      0   Skip BFGS
   3  1.45e+03  4.69e+02  1.58e+00  2.76e+03    1.71e+01      0   Skip BFGS
   4  7.26e+02  1.74e+02  1.72e+00  1.42e+03    1.66e+01      0   Skip BFGS
   5  3.63e+02  6.31e+01  1.82e+00  7.24e+02    1.64e+01      0   Skip BFGS
   6  1.82e+02  2.34e+01  1.90e+00  3.68e+02    1.51e+01      0   Skip BFGS
   7  9.08e+01  1.05e+01  1.94e+00  1.87e+02    1.27e+01      0   Skip BFGS
Convergence with smooth l2-norm regularization: Start IRLS steps...
L[p qx qy qz]-norm : [ 0.  0.  2.  2.]
eps_p: 0.0167124628752 eps_q: 0.00735284148127
Regularization decrease: 8.956e-01
   8  9.08e+01  6.57e+00  1.97e+00  1.86e+02    1.92e+01      0   Skip BFGS
   9  9.08e+01  7.90e+00  1.11e+00  1.09e+02    1.92e+01      0
  10  9.08e+01  9.07e+00  9.42e-01  9.46e+01    1.84e+01      0
Regularization decrease: 4.147e-01
  11  9.08e+01  9.50e+00  9.42e-01  9.51e+01    1.93e+01      0
  12  9.08e+01  1.20e+01  6.41e-01  7.03e+01    1.83e+01      0
  13  9.08e+01  1.17e+01  6.33e-01  6.92e+01    1.66e+01      0
Regularization decrease: 2.710e-01
  14  7.81e+01  1.16e+01  6.33e-01  6.10e+01    1.97e+01      0
  15  7.81e+01  1.17e+01  5.27e-01  5.29e+01    1.70e+01      0
  16  7.81e+01  1.06e+01  5.27e-01  5.18e+01    1.81e+01      0
Regularization decrease: 1.449e-01
  17  7.35e+01  1.06e+01  5.27e-01  4.93e+01    1.61e+01      0
  18  7.35e+01  1.06e+01  4.56e-01  4.42e+01    1.85e+01      0
  19  7.35e+01  1.01e+01  4.57e-01  4.37e+01    1.71e+01      0
Regularization decrease: 8.661e-02
  20  7.35e+01  1.01e+01  4.57e-01  4.37e+01    1.90e+01      0
  21  7.35e+01  9.94e+00  4.00e-01  3.94e+01    1.84e+01      0
  22  7.35e+01  9.71e+00  3.99e-01  3.91e+01    1.68e+01      0
Regularization decrease: 1.267e-01
  23  7.35e+01  9.65e+00  3.99e-01  3.90e+01    1.76e+01      0
  24  7.35e+01  9.68e+00  3.78e-01  3.75e+01    1.91e+01      0
  25  7.35e+01  9.61e+00  3.75e-01  3.72e+01    1.77e+01      0
Regularization decrease: 1.562e-02
  26  7.77e+01  9.46e+00  3.75e-01  3.86e+01    1.74e+01      0
  27  7.77e+01  9.74e+00  3.45e-01  3.65e+01    1.63e+01      0
  28  7.77e+01  9.46e+00  3.43e-01  3.61e+01    1.44e+01      0
Regularization decrease: 5.361e-02
  29  7.77e+01  9.51e+00  3.43e-01  3.62e+01    1.29e+01      0
  30  7.77e+01  9.62e+00  3.23e-01  3.48e+01    1.82e+01      0
  31  7.77e+01  9.57e+00  3.22e-01  3.46e+01    1.15e+01      0
Regularization decrease: 5.457e-02
  32  7.77e+01  9.63e+00  3.22e-01  3.47e+01    8.97e+00      0
  33  7.77e+01  9.77e+00  3.17e-01  3.45e+01    1.49e+01      0
  34  7.77e+01  9.70e+00  3.18e-01  3.44e+01    1.32e+01      0
Regularization decrease: 3.906e-03
Minimum decrease in regularization. End of IRLS
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 5.3623e+04
1 : |xc-x_last| = 4.4467e-03 <= tolX*(1+|x0|) = 1.0010e-01
0 : |proj(x-g)-x|    = 1.3225e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.3225e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =     100    <= iter          =     35
------------------------- DONE! -------------------------
Final misfit:9.72478514183

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
from SimPEG import Maps


def run(N=100, plotIt=True):

    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.W = 1./wd

    betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e-2)

    # Creat reduced identity map
    idenMap = Maps.IdentityMap(nP=mesh.nC)

    reg = Regularization.Sparse(mesh, mapping=idenMap)
    reg.mref = mref
    reg.cell_weights = wr
    reg.norms = [0., 0., 2., 2.]
    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

    IRLS = Directives.Update_IRLS(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(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, IRLS.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()

Total running time of the script: ( 0 minutes 4.515 seconds)

Generated by Sphinx-Gallery