.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/user-guide/examples/01-maps/plot_sumMap.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_content_user-guide_examples_01-maps_plot_sumMap.py: Maps: ComboMaps =============== Invert synthetic magnetic data with variable background values and a single block anomaly buried at depth. We will use the Sum Map to invert for both the background values and an heterogeneous susceptibiilty model. .. code-block:: python :linenos: .. GENERATED FROM PYTHON SOURCE LINES 15-208 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /content/user-guide/examples/01-maps/images/sphx_glr_plot_sumMap_001.png :alt: plot sumMap :srcset: /content/user-guide/examples/01-maps/images/sphx_glr_plot_sumMap_001.png :class: sphx-glr-multi-img * .. image-sg:: /content/user-guide/examples/01-maps/images/sphx_glr_plot_sumMap_002.png :alt: plot sumMap :srcset: /content/user-guide/examples/01-maps/images/sphx_glr_plot_sumMap_002.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none Running inversion with SimPEG v0.25.2 ================================================= Projected GNCG ================================================= # beta phi_d phi_m f |proj(x-g)-x| LS iter_CG CG |Ax-b|/|b| CG |Ax-b| Comment ----------------------------------------------------------------------------------------------------------------- 0 5.33e+05 9.11e+06 4.27e-04 9.11e+06 0 inf inf 1 5.33e+05 1.30e+05 7.92e-02 1.72e+05 6.29e+01 0 3 7.68e-04 1.06e+06 2 2.67e+05 4.26e+04 2.45e-01 1.08e+05 3.53e+01 0 9 5.83e-04 2.61e+04 3 1.33e+05 1.81e+04 3.73e-01 6.78e+04 4.91e+01 0 10 1.04e-02 1.47e+04 4 6.67e+04 5.79e+03 4.83e-01 3.80e+04 5.55e+01 0 9 5.12e-04 5.53e+03 5 3.33e+04 4.35e+03 5.15e-01 2.15e+04 4.80e+01 1 10 1.21e-02 5.80e+03 6 1.67e+04 9.04e+02 6.03e-01 1.09e+04 3.57e+01 0 10 5.58e-04 6.26e+03 7 8.33e+03 8.98e+02 6.03e-01 5.93e+03 4.87e+01 4 10 8.93e-02 1.27e+04 8 4.17e+03 4.67e+02 6.46e-01 3.16e+03 5.35e+01 0 10 1.15e-01 1.34e+05 9 2.08e+03 4.66e+02 6.46e-01 1.81e+03 3.35e+01 7 10 1.24e-01 3.85e+04 10 1.04e+03 4.19e+02 6.64e-01 1.11e+03 3.52e+01 0 10 1.50e-01 5.43e+04 11 5.21e+02 4.17e+02 6.66e-01 7.64e+02 5.19e+01 2 10 6.32e-01 6.12e+04 12 2.60e+02 4.15e+02 6.68e-01 5.89e+02 6.08e+01 3 10 1.49e-01 9.27e+04 13 1.30e+02 4.15e+02 6.68e-01 5.02e+02 3.71e+01 13 10 4.68e-02 3.43e+04 14 6.51e+01 4.13e+02 6.68e-01 4.57e+02 3.71e+01 3 10 2.31e-02 1.69e+04 15 3.26e+01 4.04e+02 6.98e-01 4.27e+02 5.49e+01 0 10 2.27e-01 6.21e+04 16 1.63e+01 3.97e+02 6.98e-01 4.09e+02 3.90e+01 0 10 6.17e-03 7.88e+03 Reached starting chifact with l2-norm regularization: Start IRLS steps... irls_threshold 0.010169372574344441 irls_threshold 0.013260669096006662 17 1.63e+01 3.93e+02 9.90e-01 4.09e+02 3.87e+01 2 10 1.27e-02 1.39e+04 18 1.63e+01 3.93e+02 1.08e+00 4.10e+02 5.98e+01 1 10 2.19e-02 9.75e+03 19 1.63e+01 3.93e+02 1.14e+00 4.11e+02 4.71e+01 10 10 1.29e-01 5.23e+04 20 1.63e+01 3.93e+02 1.18e+00 4.12e+02 4.70e+01 6 10 3.06e-02 1.25e+04 21 1.63e+01 3.93e+02 1.21e+00 4.12e+02 4.95e+01 5 10 5.51e-01 2.02e+05 22 1.63e+01 3.91e+02 1.22e+00 4.11e+02 4.49e+01 3 10 4.11e-01 1.87e+05 23 1.63e+01 3.91e+02 1.21e+00 4.11e+02 5.79e+01 11 10 7.17e-02 2.13e+04 24 1.63e+01 3.91e+02 1.18e+00 4.10e+02 5.80e+01 5 10 8.83e-02 2.64e+04 25 1.63e+01 3.89e+02 1.17e+00 4.08e+02 4.69e+01 1 10 2.62e-01 2.80e+04 26 1.63e+01 3.79e+02 1.11e+00 3.97e+02 3.78e+01 0 10 1.47e-02 1.84e+04 27 1.63e+01 3.79e+02 1.06e+00 3.96e+02 4.50e+01 5 10 1.54e-01 1.85e+04 28 1.63e+01 3.79e+02 9.83e-01 3.95e+02 5.29e+01 2 10 1.62e-01 6.36e+03 29 1.63e+01 3.79e+02 9.03e-01 3.93e+02 3.71e+01 6 10 4.40e-02 5.43e+03 30 1.63e+01 3.78e+02 8.22e-01 3.92e+02 3.74e+01 2 10 9.31e-02 1.29e+04 31 1.63e+01 3.76e+02 7.55e-01 3.89e+02 6.11e+01 0 10 2.69e-02 8.54e+03 32 1.63e+01 3.75e+02 6.82e-01 3.86e+02 6.28e+01 0 10 1.73e-02 9.40e+03 33 1.63e+01 3.75e+02 6.18e-01 3.85e+02 3.78e+01 3 10 2.61e-02 1.03e+04 34 1.63e+01 3.74e+02 5.46e-01 3.83e+02 3.76e+01 2 10 2.87e-02 1.20e+04 35 1.63e+01 3.72e+02 4.90e-01 3.80e+02 6.12e+01 0 10 3.11e-02 1.02e+04 36 1.63e+01 3.72e+02 4.28e-01 3.79e+02 3.85e+01 0 10 2.17e-02 1.08e+04 Reach maximum number of IRLS cycles: 20 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 3.1406e-01 <= tolF*(1+|f0|) = 9.1101e+05 1 : |xc-x_last| = 2.9288e-03 <= tolX*(1+|x0|) = 1.0075e-01 0 : |proj(x-g)-x| = 3.8450e+01 <= tolG = 1.0000e-03 0 : |proj(x-g)-x| = 3.8450e+01 <= 1e3*eps = 1.0000e-03 0 : maxIter = 100 <= iter = 36 ------------------------- DONE! ------------------------- | .. code-block:: Python from discretize import TensorMesh from discretize.utils import active_from_xyz from simpeg import ( utils, maps, regularization, data_misfit, optimization, inverse_problem, directives, inversion, ) from simpeg.potential_fields import magnetics import numpy as np import matplotlib.pyplot as plt def run(plotIt=True): h0_amplitude, h0_inclination, h0_declination = (50000.0, 90.0, 0.0) # Create a mesh dx = 5.0 hxind = [(dx, 5, -1.3), (dx, 10), (dx, 5, 1.3)] hyind = [(dx, 5, -1.3), (dx, 10), (dx, 5, 1.3)] hzind = [(dx, 5, -1.3), (dx, 10)] mesh = TensorMesh([hxind, hyind, hzind], "CCC") # Lets create a simple Gaussian topo and set the active cells [xx, yy] = np.meshgrid(mesh.nodes_x, mesh.nodes_y) zz = -np.exp((xx**2 + yy**2) / 75**2) + mesh.nodes_z[-1] # We would usually load a topofile topo = np.c_[utils.mkvc(xx), utils.mkvc(yy), utils.mkvc(zz)] # Go from topo to array of indices of active cells actv = active_from_xyz(mesh, topo, "N") nC = int(actv.sum()) # Create and array of observation points xr = np.linspace(-20.0, 20.0, 20) yr = np.linspace(-20.0, 20.0, 20) X, Y = np.meshgrid(xr, yr) # Move the observation points 5m above the topo Z = -np.exp((X**2 + Y**2) / 75**2) + mesh.nodes_z[-1] + 5.0 # Create a MAGsurvey rxLoc = np.c_[utils.mkvc(X.T), utils.mkvc(Y.T), utils.mkvc(Z.T)] rxLoc = magnetics.Point(rxLoc) srcField = magnetics.UniformBackgroundField( receiver_list=[rxLoc], amplitude=h0_amplitude, inclination=h0_inclination, declination=h0_declination, ) survey = magnetics.Survey(srcField) # We can now create a susceptibility model and generate data model = np.zeros(mesh.nC) # Change values in half the domain model[mesh.gridCC[:, 0] < 0] = 0.01 # Add a block in half-space model = utils.model_builder.add_block( mesh.gridCC, model, np.r_[-10, -10, 20], np.r_[10, 10, 40], 0.05 ) model = utils.mkvc(model) model = model[actv] # Create active map to go from reduce set to full actvMap = maps.InjectActiveCells(mesh, actv, np.nan) # Create reduced identity map idenMap = maps.IdentityMap(nP=nC) # Create the forward model operator prob = magnetics.Simulation3DIntegral( mesh, survey=survey, chiMap=idenMap, active_cells=actv, store_sensitivities="forward_only", ) # Compute linear forward operator and compute some data data = prob.make_synthetic_data( model, relative_error=0.0, noise_floor=1, add_noise=True ) # Create a homogenous maps for the two domains domains = [mesh.gridCC[actv, 0] < 0, mesh.gridCC[actv, 0] >= 0] homogMap = maps.SurjectUnits(domains) # Create a wire map for a second model space, voxel based wires = maps.Wires(("homo", len(domains)), ("hetero", nC)) # Create Sum map sumMap = maps.SumMap([homogMap * wires.homo, wires.hetero]) # Create the forward model operator prob = magnetics.Simulation3DIntegral( mesh, survey=survey, chiMap=sumMap, active_cells=actv, store_sensitivities="ram" ) # Make sensitivity weighting # Take the cell number out of the scaling. # Want to keep high sens for large volumes wr = ( prob.getJtJdiag(np.ones(sumMap.shape[1])) / np.r_[homogMap.P.T * mesh.cell_volumes[actv], mesh.cell_volumes[actv]] ** 2.0 ) # Scale the model spaces independently wr[wires.homo.index] /= np.max((wires.homo * wr)) * utils.mkvc( homogMap.P.sum(axis=0).flatten() ) wr[wires.hetero.index] /= np.max(wires.hetero * wr) wr = wr**0.5 ## Create a regularization # For the homogeneous model regMesh = TensorMesh([len(domains)]) reg_m1 = regularization.Sparse(regMesh, mapping=wires.homo) reg_m1.set_weights(weights=wires.homo * wr) reg_m1.norms = [0, 2] reg_m1.reference_model = np.zeros(sumMap.shape[1]) # Regularization for the voxel model reg_m2 = regularization.Sparse( mesh, active_cells=actv, mapping=wires.hetero, gradient_type="components" ) reg_m2.set_weights(weights=wires.hetero * wr) reg_m2.norms = [0, 0, 0, 0] reg_m2.reference_model = np.zeros(sumMap.shape[1]) reg = reg_m1 + reg_m2 # Data misfit function dmis = data_misfit.L2DataMisfit(simulation=prob, data=data) # Add directives to the inversion opt = optimization.ProjectedGNCG( maxIter=100, lower=0.0, upper=1.0, maxIterLS=20, cg_maxiter=10, cg_rtol=1e-3, tolG=1e-3, eps=1e-6, ) invProb = inverse_problem.BaseInvProblem(dmis, reg, opt) betaest = directives.BetaEstimate_ByEig(beta0_ratio=1e-2) # Here is where the norms are applied # Use pick a threshold parameter empirically based on the distribution of # model parameters IRLS = directives.UpdateIRLS(f_min_change=1e-3) update_Jacobi = directives.UpdatePreconditioner() inv = inversion.BaseInversion(invProb, directiveList=[IRLS, betaest, update_Jacobi]) # Run the inversion m0 = np.ones(sumMap.shape[1]) * 1e-4 # Starting model prob.model = m0 mrecSum = inv.run(m0) if plotIt: mesh.plot_3d_slicer( actvMap * model, aspect="equal", zslice=30, pcolor_opts={"cmap": "inferno_r"}, transparent="slider", ) mesh.plot_3d_slicer( actvMap * sumMap * mrecSum, aspect="equal", zslice=30, pcolor_opts={"cmap": "inferno_r"}, transparent="slider", ) if __name__ == "__main__": run() plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 32.964 seconds) **Estimated memory usage:** 346 MB .. _sphx_glr_download_content_user-guide_examples_01-maps_plot_sumMap.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_sumMap.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_sumMap.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_sumMap.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_