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 std of 5%
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***
Approximated diag(JtJ) with linear operator
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.15e+06  5.36e+05  1.36e-07  5.36e+05    2.00e+01      0
   1  5.74e+05  1.82e+05  8.95e-02  2.33e+05    1.94e+01      0
   2  2.87e+05  1.20e+05  1.67e-01  1.68e+05    1.92e+01      0   Skip BFGS
   3  1.43e+05  7.27e+04  2.83e-01  1.13e+05    1.91e+01      0   Skip BFGS
   4  7.17e+04  4.11e+04  4.38e-01  7.24e+04    1.89e+01      0   Skip BFGS
   5  3.59e+04  2.16e+04  6.27e-01  4.41e+04    1.86e+01      0   Skip BFGS
   6  1.79e+04  1.06e+04  8.41e-01  2.57e+04    1.85e+01      0   Skip BFGS
   7  8.96e+03  4.92e+03  1.06e+00  1.44e+04    1.83e+01      0   Skip BFGS
   8  4.48e+03  2.13e+03  1.28e+00  7.85e+03    1.82e+01      0   Skip BFGS
   9  2.24e+03  8.56e+02  1.47e+00  4.16e+03    1.77e+01      0   Skip BFGS
  10  1.12e+03  3.25e+02  1.63e+00  2.16e+03    1.70e+01      0   Skip BFGS
  11  5.60e+02  1.19e+02  1.76e+00  1.11e+03    1.67e+01      0   Skip BFGS
  12  2.80e+02  4.31e+01  1.85e+00  5.62e+02    1.61e+01      0   Skip BFGS
  13  1.40e+02  1.68e+01  1.92e+00  2.85e+02    1.42e+01      0   Skip BFGS
Reached starting chifact with l2-norm regularization: Start IRLS steps...
eps_p: 1.2115885600793632 eps_q: 1.2115885600793632
delta phim: 1.580e+07
  14  7.00e+01  8.44e+00  2.57e+00  1.89e+02    7.19e+00      0   Skip BFGS
Beta search step
  14  1.22e+02  8.44e+00  2.57e+00  3.24e+02    1.60e+01      0
Beta search step
  14  9.83e+01  8.44e+00  2.57e+00  2.62e+02    1.37e+01      0
delta phim: 4.203e-01
  15  9.83e+01  1.00e+01  2.70e+00  2.76e+02    1.78e+01      0   Skip BFGS
Beta search step
  15  7.82e+01  1.00e+01  2.70e+00  2.21e+02    1.63e+01      0
delta phim: 2.005e-01
  16  7.82e+01  9.39e+00  2.73e+00  2.23e+02    1.91e+01      0   Skip BFGS
Beta search step
  16  6.40e+01  9.39e+00  2.73e+00  1.85e+02    1.70e+01      0
delta phim: 2.136e-01
  17  6.40e+01  9.10e+00  2.71e+00  1.83e+02    1.54e+01      0   Skip BFGS
delta phim: 1.905e-01
  18  6.40e+01  1.01e+01  2.63e+00  1.79e+02    1.74e+01      0
Beta search step
  18  5.31e+01  1.01e+01  2.63e+00  1.50e+02    1.76e+01      0
delta phim: 1.640e-01
  19  5.31e+01  9.68e+00  2.51e+00  1.43e+02    1.54e+01      0   Skip BFGS
delta phim: 1.596e-01
  20  5.31e+01  1.04e+01  2.34e+00  1.35e+02    1.77e+01      0
Beta search step
  20  4.37e+01  1.04e+01  2.34e+00  1.13e+02    1.63e+01      0
delta phim: 1.351e-01
  21  4.37e+01  1.00e+01  2.17e+00  1.05e+02    1.56e+01      0   Skip BFGS
delta phim: 1.347e-01
  22  4.37e+01  1.02e+01  1.95e+00  9.55e+01    1.63e+01      0
delta phim: 1.114e-01
  23  4.37e+01  1.06e+01  1.75e+00  8.70e+01    1.67e+01      0
delta phim: 1.038e-01
  24  4.37e+01  1.08e+01  1.53e+00  7.79e+01    1.69e+01      0
Beta search step
  24  3.59e+01  1.08e+01  1.53e+00  6.58e+01    1.62e+01      0
delta phim: 8.640e-02
  25  3.59e+01  1.04e+01  1.36e+00  5.90e+01    1.56e+01      0   Skip BFGS
delta phim: 9.022e-02
  26  3.59e+01  1.01e+01  1.17e+00  5.23e+01    1.60e+01      0
delta phim: 6.574e-02
  27  3.59e+01  1.03e+01  1.01e+00  4.66e+01    1.52e+01      0
delta phim: 5.326e-02
  28  3.59e+01  1.01e+01  8.59e-01  4.10e+01    1.57e+01      0
delta phim: 3.247e-02
  29  3.59e+01  1.02e+01  7.29e-01  3.64e+01    1.43e+01      0
delta phim: 3.103e-02
  30  3.59e+01  1.00e+01  6.15e-01  3.21e+01    1.55e+01      0
delta phim: 1.467e-02
  31  3.59e+01  9.95e+00  5.15e-01  2.84e+01    1.44e+01      0
delta phim: 1.209e-02
  32  3.59e+01  9.72e+00  4.31e-01  2.52e+01    1.51e+01      0
delta phim: 2.704e-03
  33  3.59e+01  9.59e+00  3.59e-01  2.25e+01    1.44e+01      0
delta phim: 5.466e-03
  34  3.59e+01  9.44e+00  3.00e-01  2.02e+01    1.50e+01      0
delta phim: 2.864e-03
  35  3.59e+01  9.36e+00  2.50e-01  1.83e+01    1.37e+01      0
delta phim: 5.953e-03
  36  3.59e+01  9.15e+00  2.13e-01  1.68e+01    1.68e+01      0
delta phim: 1.450e-02
  37  3.59e+01  9.06e+00  1.80e-01  1.55e+01    1.35e+01      0
delta phim: 2.079e-02
  38  3.59e+01  9.03e+00  1.46e-01  1.43e+01    1.67e+01      0
Beta search step
  38  5.61e+01  9.03e+00  1.46e-01  1.72e+01    1.59e+01      0
delta phim: 3.068e-02
  39  5.61e+01  9.20e+00  1.19e-01  1.59e+01    1.57e+01      0
delta phim: 3.312e-02
  40  5.61e+01  9.16e+00  1.00e-01  1.48e+01    1.54e+01      0
delta phim: 1.273e-02
  41  5.61e+01  9.21e+00  8.20e-02  1.38e+01    1.75e+01      0
delta phim: 2.689e-02
  42  5.61e+01  9.11e+00  6.91e-02  1.30e+01    1.65e+01      0
delta phim: 5.760e-03
  43  5.61e+01  9.09e+00  5.88e-02  1.24e+01    1.60e+01      0
Beta search step
  43  8.74e+01  9.09e+00  5.88e-02  1.42e+01    1.84e+01      0
delta phim: 1.269e-02
  44  8.74e+01  9.19e+00  4.80e-02  1.34e+01    1.67e+01      0
delta phim: 2.146e-02
  45  8.74e+01  9.21e+00  3.91e-02  1.26e+01    1.45e+01      0
delta phim: 5.622e-02
  46  8.74e+01  9.27e+00  3.36e-02  1.22e+01    1.53e+01      0
delta phim: 1.443e-02
  47  8.74e+01  9.20e+00  2.98e-02  1.18e+01    1.90e+01      0
delta phim: 4.609e-02
  48  8.74e+01  9.26e+00  2.49e-02  1.14e+01    1.46e+01      0
delta phim: 1.249e-02
  49  8.74e+01  9.23e+00  2.13e-02  1.11e+01    1.55e+01      0
delta phim: 1.455e-02
  50  8.74e+01  9.25e+00  1.81e-02  1.08e+01    1.46e+01      0
delta phim: 1.147e-02
  51  8.74e+01  9.20e+00  1.57e-02  1.06e+01    1.63e+01      0
delta phim: 3.844e-02
  52  8.74e+01  9.17e+00  1.32e-02  1.03e+01    1.47e+01      0
delta phim: 5.251e-03
  53  8.74e+01  9.12e+00  1.11e-02  1.01e+01    1.61e+01      0
Reach maximum number of IRLS cycles: 40
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 5.3623e+04
1 : |xc-x_last| = 3.5634e-03 <= tolX*(1+|x0|) = 1.0010e-01
0 : |proj(x-g)-x|    = 1.6059e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.6059e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =     100    <= iter          =     54
------------------------- DONE! -------------------------
Final misfit:9.082738438198374

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.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=40, 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 13.138 seconds)

Gallery generated by Sphinx-Gallery