EM: TDEM: 1D: Inversion

Here we will create and run a TDEM 1D inversion.

Out:

SimPEG.Survey assigned new std of 5.00%
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***
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.81e+03  2.60e+03  0.00e+00  2.60e+03    3.13e+03      0
   1  2.81e+03  2.61e+02  7.46e-02  4.71e+02    4.03e+02      0
   2  5.63e+02  1.07e+02  1.08e-01  1.68e+02    2.85e+02      0   Skip BFGS
   3  5.63e+02  1.96e+01  1.68e-01  1.14e+02    1.28e+01      0   Skip BFGS
   4  1.13e+02  1.86e+01  1.69e-01  3.77e+01    7.46e+01      0   Skip BFGS
   5  1.13e+02  9.75e+00  1.97e-01  3.19e+01    2.95e+00      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 5.7587e+00 <= tolF*(1+|f0|) = 2.6055e+02
1 : |xc-x_last| = 2.0393e-01 <= tolX*(1+|x0|) = 3.0149e+00
0 : |proj(x-g)-x|    = 2.9460e+00 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 2.9460e+00 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =       5    <= iter          =      5
------------------------- DONE! -------------------------
/Users/lindseyjh/git/simpeg/simpeg/examples/09-tdem/plot_inversion_1D.py:72: UserWarning: Data has no positive values, and therefore cannot be log-scaled.
  ax[0].loglog(rx.times, survey.dtrue, 'b.-')
/Users/lindseyjh/git/simpeg/simpeg/examples/09-tdem/plot_inversion_1D.py:73: UserWarning: Data has no positive values, and therefore cannot be log-scaled.
  ax[0].loglog(rx.times, survey.dobs, 'r.-')
/Users/lindseyjh/git/simpeg/simpeg/examples/09-tdem/plot_inversion_1D.py:92: UserWarning: Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.
  plt.show()

import numpy as np
from SimPEG import (Mesh, Maps, SolverLU, DataMisfit, Regularization,
                    Optimization, InvProblem, Inversion, Directives, Utils)
import SimPEG.EM as EM
import matplotlib.pyplot as plt


def run(plotIt=True):

    cs, ncx, ncz, npad = 5., 25, 15, 15
    hx = [(cs, ncx),  (cs, npad, 1.3)]
    hz = [(cs, npad, -1.3), (cs, ncz), (cs, npad, 1.3)]
    mesh = Mesh.CylMesh([hx, 1, hz], '00C')

    active = mesh.vectorCCz < 0.
    layer = (mesh.vectorCCz < 0.) & (mesh.vectorCCz >= -100.)
    actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
    sig_half = 2e-3
    sig_air = 1e-8
    sig_layer = 1e-3
    sigma = np.ones(mesh.nCz)*sig_air
    sigma[active] = sig_half
    sigma[layer] = sig_layer
    mtrue = np.log(sigma[active])

    rxOffset = 1e-3
    rx = EM.TDEM.Rx.Point_dbdt(
        np.array([[rxOffset, 0., 30]]),
        np.logspace(-5, -3, 31),
        'z'
    )
    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 80]))
    survey = EM.TDEM.Survey([src])
    prb = EM.TDEM.Problem3D_e(mesh, sigmaMap=mapping)

    prb.Solver = SolverLU
    prb.timeSteps = [(1e-06, 20), (1e-05, 20), (0.0001, 20)]
    prb.pair(survey)

    # create observed data
    std = 0.05

    survey.dobs = survey.makeSyntheticData(mtrue, std)
    survey.std = std
    survey.eps = 1e-5*np.linalg.norm(survey.dobs)

    dmisfit = DataMisfit.l2_DataMisfit(survey)
    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
    reg = Regularization.Tikhonov(regMesh, alpha_s=1e-2, alpha_x=1.)
    opt = Optimization.InexactGaussNewton(maxIter=5, LSshorten=0.5)
    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)

    # Create an inversion object
    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
    betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest])
    m0 = np.log(np.ones(mtrue.size)*sig_half)
    prb.counter = opt.counter = Utils.Counter()
    opt.remember('xc')

    mopt = inv.run(m0)

    if plotIt:
        fig, ax = plt.subplots(1, 2, figsize=(10, 6))
        ax[0].loglog(rx.times, survey.dtrue, 'b.-')
        ax[0].loglog(rx.times, survey.dobs, 'r.-')
        ax[0].legend(('Noisefree', '$d^{obs}$'), fontsize=16)
        ax[0].set_xlabel('Time (s)', fontsize=14)
        ax[0].set_ylabel('$B_z$ (T)', fontsize=16)
        ax[0].set_xlabel('Time (s)', fontsize=14)
        ax[0].grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)

        plt.semilogx(sigma[active], mesh.vectorCCz[active])
        plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])
        ax[1].set_ylim(-600, 0)
        ax[1].set_xlim(1e-4, 1e-2)
        ax[1].set_xlabel('Conductivity (S/m)', fontsize=14)
        ax[1].set_ylabel('Depth (m)', fontsize=14)
        ax[1].grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
        plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'])


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