Inversion: Linear ProblemΒΆ

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

../../../_images/sphx_glr_plot_inversion_linear_001.png

Out:

Efficiency Warning: Interpolation will be slow, use setup.py!

            python setup.py build_ext --inplace

SimPEG.InvProblem will set Regularization.mref to m0.
SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                    ***Done using same Solver and solverOpts as the problem***
SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||
model has any nan: 0
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
   0  2.07e+01  7.99e+04  0.00e+00  7.99e+04    3.59e+05      0
   1  2.07e+01  3.36e+04  9.61e-01  3.36e+04    1.59e+05      0
   2  2.07e+01  2.13e+04  2.83e+00  2.13e+04    8.15e+04      0
   3  2.07e+01  7.48e+03  1.24e+01  7.74e+03    7.59e+04      0   Skip BFGS
   4  2.07e+01  5.85e+03  1.23e+01  6.10e+03    7.69e+04      0
   5  2.07e+01  5.41e+03  1.46e+01  5.72e+03    7.42e+04      0
   6  2.07e+01  5.17e+03  1.49e+01  5.47e+03    8.20e+04      0
   7  2.07e+01  5.00e+03  1.53e+01  5.31e+03    8.04e+04      0
   8  2.07e+01  4.78e+03  1.53e+01  5.10e+03    7.04e+04      0
   9  2.07e+01  4.29e+03  1.72e+01  4.65e+03    8.46e+04      0
  10  2.07e+01  3.98e+03  1.76e+01  4.34e+03    7.22e+04      0
  11  2.07e+01  3.67e+03  1.90e+01  4.06e+03    7.78e+04      0
  12  2.07e+01  3.61e+03  1.90e+01  4.01e+03    8.30e+04      0
  13  2.07e+01  3.57e+03  1.94e+01  3.97e+03    7.71e+04      0
  14  2.07e+01  3.06e+02  2.14e+01  7.50e+02    4.88e+04      0   Skip BFGS
  15  2.07e+01  2.00e+02  2.10e+01  6.34e+02    7.19e+03      0
  16  2.07e+01  1.65e+02  2.11e+01  6.02e+02    6.74e+03      0
  17  2.07e+01  1.58e+02  2.07e+01  5.86e+02    1.02e+04      0   Skip BFGS
  18  2.07e+01  1.42e+02  2.10e+01  5.77e+02    5.33e+03      0
  19  2.07e+01  1.19e+02  2.12e+01  5.58e+02    4.27e+03      0   Skip BFGS
  20  2.07e+01  1.06e+02  2.15e+01  5.52e+02    4.05e+03      0
  21  2.07e+01  1.03e+02  2.16e+01  5.51e+02    4.43e+03      0   Skip BFGS
  22  2.07e+01  1.02e+02  2.16e+01  5.50e+02    4.10e+03      0
  23  2.07e+01  1.00e+02  2.17e+01  5.50e+02    4.02e+03      0   Skip BFGS
  24  2.07e+01  1.02e+02  2.16e+01  5.50e+02    4.04e+03      0
  25  2.07e+01  9.82e+01  2.16e+01  5.46e+02    2.44e+03      0   Skip BFGS
  26  2.07e+01  9.42e+01  2.17e+01  5.43e+02    3.97e+03      0
  27  2.07e+01  6.40e+01  2.21e+01  5.22e+02    5.13e+03      0   Skip BFGS
  28  2.07e+01  6.21e+01  2.21e+01  5.20e+02    4.99e+03      0
  29  2.07e+01  6.18e+01  2.21e+01  5.20e+02    4.95e+03      0   Skip BFGS
  30  2.07e+01  6.20e+01  2.21e+01  5.20e+02    4.98e+03      0
  31  2.07e+01  4.82e+01  2.24e+01  5.12e+02    3.50e+03      0   Skip BFGS
  32  2.07e+01  4.81e+01  2.22e+01  5.07e+02    3.08e+03      0
  33  2.07e+01  4.06e+01  2.21e+01  4.98e+02    2.24e+03      0   Skip BFGS
  34  2.07e+01  3.80e+01  2.21e+01  4.96e+02    1.45e+03      0
  35  2.07e+01  3.26e+01  2.22e+01  4.92e+02    2.43e+03      0
  36  2.07e+01  3.12e+01  2.22e+01  4.91e+02    1.83e+03      0   Skip BFGS
  37  2.07e+01  3.17e+01  2.22e+01  4.91e+02    1.82e+03      0
  38  2.07e+01  3.25e+01  2.21e+01  4.91e+02    1.92e+03      0
  39  2.07e+01  3.23e+01  2.21e+01  4.91e+02    1.76e+03      0
  40  2.07e+01  3.14e+01  2.22e+01  4.91e+02    1.64e+03      0
  41  2.07e+01  3.14e+01  2.22e+01  4.91e+02    1.60e+03      0   Skip BFGS
  42  2.07e+01  3.11e+01  2.22e+01  4.91e+02    1.79e+03      0
  43  2.07e+01  3.18e+01  2.22e+01  4.91e+02    2.03e+03      0   Skip BFGS
  44  2.07e+01  3.12e+01  2.22e+01  4.91e+02    2.20e+03      0
  45  2.07e+01  3.12e+01  2.22e+01  4.91e+02    2.00e+03      0
  46  2.07e+01  2.60e+01  2.24e+01  4.89e+02    2.45e+03      0   Skip BFGS
  47  2.07e+01  2.81e+01  2.22e+01  4.89e+02    2.37e+03      0
  48  2.07e+01  2.82e+01  2.22e+01  4.88e+02    1.57e+03      0   Skip BFGS
  49  2.07e+01  2.74e+01  2.23e+01  4.88e+02    1.61e+03      0
  50  2.07e+01  2.57e+01  2.22e+01  4.86e+02    5.35e+02      0   Skip BFGS
  51  2.07e+01  2.55e+01  2.22e+01  4.86e+02    5.51e+02      0
  52  2.07e+01  2.56e+01  2.22e+01  4.86e+02    4.65e+02      0
  53  2.07e+01  2.55e+01  2.22e+01  4.86e+02    2.44e+02      0   Skip BFGS
  54  2.07e+01  2.55e+01  2.22e+01  4.86e+02    2.59e+02      0
  55  2.07e+01  2.55e+01  2.22e+01  4.86e+02    2.73e+02      0   Skip BFGS
  56  2.07e+01  2.55e+01  2.22e+01  4.86e+02    2.68e+02      0
  57  2.07e+01  2.56e+01  2.22e+01  4.86e+02    2.68e+02      0   Skip BFGS
  58  2.07e+01  2.56e+01  2.22e+01  4.86e+02    2.58e+02      0
  59  2.07e+01  2.55e+01  2.22e+01  4.86e+02    6.63e+02      0   Skip BFGS
  60  2.07e+01  2.54e+01  2.22e+01  4.86e+02    6.61e+02      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.7534e-03 <= tolF*(1+|f0|) = 7.9851e+03
1 : |xc-x_last| = 1.4427e-03 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x|    = 6.6098e+02 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 6.6098e+02 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      60    <= iter          =     60
------------------------- DONE! -------------------------

from __future__ import print_function
import numpy as np
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
import matplotlib.pyplot as plt


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

    np.random.seed(1)

    mesh = Mesh.TensorMesh([N])

    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.makeSyntheticData(mtrue, std=0.01)

    M = prob.mesh

    reg = Regularization.Tikhonov(mesh, alpha_s=1., alpha_x=1.)
    dmis = DataMisfit.l2_DataMisfit(survey)
    opt = Optimization.InexactGaussNewton(maxIter=60)
    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
    directives = [
        Directives.BetaEstimate_ByEig(beta0_ratio=1e-2),
        Directives.TargetMisfit()
    ]
    inv = Inversion.BaseInversion(invProb, directiveList=directives)
    m0 = np.zeros_like(survey.mtrue)

    mrec = inv.run(m0)

    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(M.vectorCCx, survey.mtrue, 'b-')
        axes[1].plot(M.vectorCCx, mrec, 'r-')
        axes[1].legend(('True Model', 'Recovered Model'))
        axes[1].set_ylim([-2, 2])

    return prob, survey, mesh, mrec

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

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

Generated by Sphinx-Gallery