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+06  5.36e+05  1.36e-07  5.36e+05    2.00e+01      0
   1  2.91e+05  1.83e+05  8.84e-02  2.09e+05    1.94e+01      0
   2  7.27e+04  7.34e+04  2.80e-01  9.38e+04    1.91e+01      0   Skip BFGS
   3  1.82e+04  2.19e+04  6.23e-01  3.33e+04    1.87e+01      0   Skip BFGS
   4  4.56e+03  5.00e+03  1.06e+00  9.83e+03    1.83e+01      0   Skip BFGS
   5  1.17e+03  8.78e+02  1.47e+00  2.59e+03    1.79e+01      0   Skip BFGS
   6  3.37e+02  1.30e+02  1.75e+00  7.19e+02    1.69e+01      0   Skip BFGS
   7  1.63e+02  2.13e+01  1.90e+00  3.32e+02    1.47e+01      0   Skip BFGS
Reached starting chifact with l2-norm regularization: Start IRLS steps...
eps_p: 0.1 eps_q: 0.1
Eps_p: 0.08333333333333334
Eps_q: 0.08333333333333334
delta phim: 1.568e+07
   8  1.63e+02  9.49e+00  1.06e+00  1.83e+02    1.89e+01      0   Skip BFGS
Beta search step
   9  1.13e+02  9.49e+00  1.06e+00  1.30e+02    1.86e+01      0
Eps_p: 0.06944444444444446
Eps_q: 0.06944444444444446
delta phim: 3.082e+00
  10  1.13e+02  9.65e+00  7.11e-01  9.00e+01    2.04e+01      0   Skip BFGS
Eps_p: 0.057870370370370385
Eps_q: 0.057870370370370385
delta phim: 1.831e-01
  11  1.13e+02  1.03e+01  5.76e-01  7.54e+01    1.60e+01      0   Skip BFGS
Eps_p: 0.04822530864197532
Eps_q: 0.04822530864197532
delta phim: 2.937e-02
  12  1.13e+02  1.00e+01  4.75e-01  6.38e+01    1.85e+01      0
Eps_p: 0.0401877572016461
Eps_q: 0.0401877572016461
delta phim: 3.171e-02
  13  1.13e+02  9.59e+00  3.77e-01  5.22e+01    1.60e+01      0
Eps_p: 0.03348979766803842
Eps_q: 0.03348979766803842
delta phim: 2.629e-02
  14  1.13e+02  9.12e+00  2.96e-01  4.25e+01    1.69e+01      0
Eps_p: 0.02790816472336535
Eps_q: 0.02790816472336535
delta phim: 6.605e-02
  15  1.13e+02  9.08e+00  2.44e-01  3.67e+01    1.51e+01      0   Skip BFGS
Beta search step
  16  1.76e+02  9.08e+00  2.44e-01  5.21e+01    1.78e+01      0
Eps_p: 0.023256803936137792
Eps_q: 0.023256803936137792
delta phim: 1.654e-02
  17  1.76e+02  9.57e+00  2.01e-01  4.50e+01    1.38e+01      0
Eps_p: 0.019380669946781493
Eps_q: 0.019380669946781493
delta phim: 1.558e-02
  18  1.76e+02  9.00e+00  1.68e-01  3.86e+01    1.78e+01      0
Eps_p: 0.016150558288984578
Eps_q: 0.016150558288984578
delta phim: 6.584e-03
  19  1.76e+02  9.30e+00  1.39e-01  3.37e+01    1.63e+01      0
Beta search step
  20  2.75e+02  9.30e+00  1.39e-01  4.75e+01    1.95e+01      0
Eps_p: 0.013458798574153816
Eps_q: 0.013458798574153816
delta phim: 1.539e-02
  21  2.75e+02  9.18e+00  1.16e-01  4.11e+01    1.63e+01      0
Eps_p: 0.011215665478461513
Eps_q: 0.011215665478461513
delta phim: 2.932e-02
  22  2.75e+02  9.25e+00  9.88e-02  3.64e+01    1.76e+01      0
Beta search step
  23  4.29e+02  9.25e+00  9.88e-02  5.16e+01    1.97e+01      0
Eps_p: 0.009346387898717928
Eps_q: 0.009346387898717928
delta phim: 1.760e-02
  24  4.29e+02  9.16e+00  8.22e-02  4.44e+01    1.93e+01      0
Eps_p: 0.00778865658226494
Eps_q: 0.00778865658226494
delta phim: 5.833e-02
  25  4.29e+02  9.08e+00  6.62e-02  3.75e+01    1.84e+01      0
Eps_p: 0.00649054715188745
Eps_q: 0.00649054715188745
delta phim: 7.242e-02
  26  4.29e+02  9.06e+00  5.24e-02  3.15e+01    1.63e+01      0
Eps_p: 0.005408789293239542
Eps_q: 0.005408789293239542
delta phim: 5.476e-02
  27  4.29e+02  9.03e+00  4.27e-02  2.73e+01    1.77e+01      0
Eps_p: 0.004507324411032952
Eps_q: 0.004507324411032952
delta phim: 2.722e-02
  28  4.29e+02  9.10e+00  3.55e-02  2.43e+01    1.63e+01      0
Eps_p: 0.0037561036758607933
Eps_q: 0.0037561036758607933
delta phim: 6.178e-04
  29  4.29e+02  9.11e+00  2.96e-02  2.18e+01    1.69e+01      0
Eps_p: 0.003130086396550661
Eps_q: 0.003130086396550661
delta phim: 3.599e-04
  30  4.29e+02  9.16e+00  2.41e-02  1.95e+01    1.79e+01      0
Eps_p: 0.0026084053304588845
Eps_q: 0.0026084053304588845
delta phim: 2.421e-02
  31  4.29e+02  9.15e+00  2.00e-02  1.77e+01    1.71e+01      0
Reach maximum number of IRLS cycles: 20
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 5.3623e+04
1 : |xc-x_last| = 3.4583e-02 <= tolX*(1+|x0|) = 1.0010e-01
0 : |proj(x-g)-x|    = 1.7130e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.7130e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =     100    <= iter          =     32
------------------------- DONE! -------------------------
Final misfit:9.178804142181734

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.getJ(m0)**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=1e0)

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

    reg = Regularization.Sparse(mesh, mapping=idenMap)
    reg.mref = mref
    reg.cell_weights = wr
    reg.eps_p, reg.eps_q = 1e-1, 1e-1
    reg.norms = np.c_[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.UpdatePreconditioner()

    # 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(
        maxIRLSiter=20, minGNiter=1, f_min_change=1e-4)
    saveDict = Directives.SaveOutputEveryIteration(save_txt=False)
    inv = Inversion.BaseInversion(
        invProb,
        directiveList=[IRLS, betaest, update_Jacobi, saveDict]
    )

    # Run inversion
    mrec = inv.run(m0)

    print("Final misfit:" + str(invProb.dmisfit(mrec)))

    if plotIt:
        fig, axes = plt.subplots(2, 2, figsize=(12*1.2, 8*1.2))
        for i in range(prob.G.shape[0]):
            axes[0, 0].plot(prob.G[i, :])
        axes[0, 0].set_title('Columns of matrix G')

        axes[0, 1].plot(mesh.vectorCCx, mtrue, 'b-')
        axes[0, 1].plot(mesh.vectorCCx, invProb.l2model, 'r-')
        # axes[0, 1].legend(('True Model', 'Recovered Model'))
        axes[0, 1].set_ylim(-1.0, 1.25)

        axes[0, 1].plot(mesh.vectorCCx, mrec, 'k-', lw=2)
        axes[0, 1].legend(
            (
                'True Model',
                'Smooth l2-l2',
                'Sparse norms: {0}'.format(*reg.norms)
            ),
            fontsize=12
        )

        axes[1, 1].plot(saveDict.phi_d, 'k', lw=2)

        twin = axes[1, 1].twinx()
        twin.plot(saveDict.phi_m, 'k--', lw=2)
        axes[1, 1].plot(
            np.r_[IRLS.iterStart, IRLS.iterStart],
            np.r_[0, np.max(saveDict.phi_d)], 'k:'
        )
        axes[1, 1].text(
            IRLS.iterStart, 0.,
            'IRLS Start', va='bottom', ha='center',
            rotation='vertical', size=12,
            bbox={'facecolor': 'white'}
        )

        axes[1, 1].set_ylabel('$\phi_d$', size=16, rotation=0)
        axes[1, 1].set_xlabel('Iterations', size=14)
        axes[1, 0].axis('off')
        twin.set_ylabel('$\phi_m$', size=16, rotation=0)

    return prob, survey, mesh, mrec


if __name__ == '__main__':
    run()
    plt.show()

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

Gallery generated by Sphinx-Gallery