.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/examples/10-pgi/plot_inv_0_PGI_Linear_1D.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_examples_10-pgi_plot_inv_0_PGI_Linear_1D.py: Petrophysically guided inversion (PGI): Linear example ====================================================== We do a comparison between the classic least-squares inversion and our formulation of a petrophysically constrained inversion. We explore it through the UBC linear example. .. GENERATED FROM PYTHON SOURCE LINES 12-14 Tikhonov Inversion# #################### .. GENERATED FROM PYTHON SOURCE LINES 14-90 .. code-block:: Python import discretize as Mesh import matplotlib.pyplot as plt import numpy as np from simpeg import ( data_misfit, directives, inverse_problem, inversion, maps, optimization, regularization, simulation, utils, ) # Random seed for reproductibility np.random.seed(1) # Mesh N = 100 mesh = Mesh.TensorMesh([N]) # Survey design parameters nk = 20 jk = np.linspace(1.0, 60.0, nk) p = -0.25 q = 0.25 # Physics def g(k): return np.exp(p * jk[k] * mesh.cell_centers_x) * np.cos( np.pi * q * jk[k] * mesh.cell_centers_x ) G = np.empty((nk, mesh.nC)) for i in range(nk): G[i, :] = g(i) # True model mtrue = np.zeros(mesh.nC) mtrue[mesh.cell_centers_x > 0.2] = 1.0 mtrue[mesh.cell_centers_x > 0.35] = 0.0 t = (mesh.cell_centers_x - 0.65) / 0.25 indx = np.abs(t) < 1 mtrue[indx] = -(((1 - t**2.0) ** 2.0)[indx]) mtrue = np.zeros(mesh.nC) mtrue[mesh.cell_centers_x > 0.3] = 1.0 mtrue[mesh.cell_centers_x > 0.45] = -0.5 mtrue[mesh.cell_centers_x > 0.6] = 0 # simpeg problem and survey prob = simulation.LinearSimulation(mesh, G=G, model_map=maps.IdentityMap()) std = 0.01 survey = prob.make_synthetic_data(mtrue, relative_error=std, add_noise=True) # Setup the inverse problem reg = regularization.WeightedLeastSquares(mesh, alpha_s=1.0, alpha_x=1.0) dmis = data_misfit.L2DataMisfit(data=survey, simulation=prob) opt = optimization.ProjectedGNCG(maxIter=10, maxIterCG=50, tolCG=1e-4) invProb = inverse_problem.BaseInvProblem(dmis, reg, opt) directiveslist = [ directives.BetaEstimate_ByEig(beta0_ratio=1e-5), directives.BetaSchedule(coolingFactor=10.0, coolingRate=2), directives.TargetMisfit(), ] inv = inversion.BaseInversion(invProb, directiveList=directiveslist) m0 = np.zeros_like(mtrue) mnormal = inv.run(m0) .. rst-class:: sphx-glr-script-out .. code-block:: none Running inversion with SimPEG v0.22.2 simpeg.InvProblem will set Regularization.reference_model to m0. simpeg.InvProblem will set Regularization.reference_model to m0. simpeg.InvProblem will set Regularization.reference_model to m0. simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv. ***Done using the default solver Pardiso and no solver_opts.*** 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.83e+01 2.00e+05 0.00e+00 2.00e+05 2.42e+06 0 1 1.83e+01 3.40e+01 4.59e+01 8.76e+02 7.54e-06 0 ------------------------- STOP! ------------------------- 0 : |fc-fOld| = 1.9912e+05 <= tolF*(1+|f0|) = 2.0000e+04 0 : |xc-x_last| = 4.0098e+00 <= tolX*(1+|x0|) = 1.0000e-01 1 : |proj(x-g)-x| = 7.5419e-06 <= tolG = 1.0000e-01 1 : |proj(x-g)-x| = 7.5419e-06 <= 1e3*eps = 1.0000e-02 0 : maxIter = 10 <= iter = 1 ------------------------- DONE! ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 91-93 Petrophysically constrained inversion # ######################################## .. GENERATED FROM PYTHON SOURCE LINES 93-168 .. code-block:: Python # fit a Gaussian Mixture Model with n components # on the true model to simulate the laboratory # petrophysical measurements n = 3 clf = utils.WeightedGaussianMixture( mesh=mesh, n_components=n, covariance_type="full", max_iter=100, n_init=3, reg_covar=5e-4, ) clf.fit(mtrue.reshape(-1, 1)) # Petrophyically constrained regularization reg = regularization.PGI( gmmref=clf, mesh=mesh, alpha_pgi=1.0, alpha_x=1.0, ) # Optimization opt = optimization.ProjectedGNCG(maxIter=20, maxIterCG=50, tolCG=1e-4) opt.remember("xc") # Setup new inverse problem invProb = inverse_problem.BaseInvProblem(dmis, reg, opt) # directives Alphas = directives.AlphasSmoothEstimate_ByEig(alpha0_ratio=10.0, verbose=True) beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-8) betaIt = directives.PGI_BetaAlphaSchedule( verbose=True, coolingFactor=2.0, warmingFactor=1.0, tolerance=0.1, update_rate=1, progress=0.2, ) targets = directives.MultiTargetMisfits(verbose=True) petrodir = directives.PGI_UpdateParameters() addmref = directives.PGI_AddMrefInSmooth(verbose=True) # Setup Inversion inv = inversion.BaseInversion( invProb, directiveList=[Alphas, beta, petrodir, targets, addmref, betaIt] ) # Initial model same as for WeightedLeastSquares mcluster = inv.run(m0) # Final Plot fig, axes = plt.subplots(1, 3, 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].hist(mtrue, bins=20, linewidth=3.0, density=True, color="k") axes[1].set_xlabel("Model value") axes[1].set_xlabel("Occurence") axes[1].hist(mnormal, bins=20, density=True, color="b") axes[1].hist(mcluster, bins=20, density=True, color="r") axes[1].legend(["Mtrue Hist.", "L2 Model Hist.", "PGI Model Hist."]) axes[2].plot(mesh.cell_centers_x, mtrue, color="black", linewidth=3) axes[2].plot(mesh.cell_centers_x, mnormal, color="blue") axes[2].plot(mesh.cell_centers_x, mcluster, "r-") axes[2].plot(mesh.cell_centers_x, invProb.reg.objfcts[0].reference_model, "r--") axes[2].legend(("True Model", "L2 Model", "PGI Model", "Learned Mref")) axes[2].set_ylim([-2, 2]) plt.show() .. image-sg:: /content/examples/10-pgi/images/sphx_glr_plot_inv_0_PGI_Linear_1D_001.png :alt: Columns of matrix G :srcset: /content/examples/10-pgi/images/sphx_glr_plot_inv_0_PGI_Linear_1D_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Running inversion with SimPEG v0.22.2 simpeg.InvProblem will set Regularization.reference_model to m0. simpeg.InvProblem will set Regularization.reference_model to m0. simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv. ***Done using the default solver Pardiso and no solver_opts.*** Alpha scales: [0.5281823872507709, 0.0] 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 3.16e-02 2.00e+05 0.00e+00 2.00e+05 2.42e+06 0 geophys. misfits: 2.0 (target 20.0 [True]) | smallness misfit: 4663.6 (target: 100.0 [False]) mref changed in 26 places Beta cooling evaluation: progress: [2.]; minimum progress targets: [160000.] Warming alpha_pgi to favor clustering: 10.23238713518336 1 3.16e-02 1.95e+00 5.13e+02 1.82e+01 1.92e+01 0 geophys. misfits: 3.2 (target 20.0 [True]) | smallness misfit: 910.7 (target: 100.0 [False]) mref changed in 0 places Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells Beta cooling evaluation: progress: [3.2]; minimum progress targets: [22.] Warming alpha_pgi to favor clustering: 63.141625636935046 2 3.16e-02 3.24e+00 6.55e+02 2.39e+01 4.36e+01 0 geophys. misfits: 4.7 (target 20.0 [True]) | smallness misfit: 478.4 (target: 100.0 [False]) mref changed in 0 places Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells Beta cooling evaluation: progress: [4.7]; minimum progress targets: [22.] Warming alpha_pgi to favor clustering: 270.2127421572115 3 3.16e-02 4.67e+00 1.42e+03 4.96e+01 2.69e+02 0 geophys. misfits: 7.4 (target 20.0 [True]) | smallness misfit: 394.0 (target: 100.0 [False]) mref changed in 0 places Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells Beta cooling evaluation: progress: [7.4]; minimum progress targets: [22.] Warming alpha_pgi to favor clustering: 725.5914794315696 4 3.16e-02 7.45e+00 3.00e+03 1.02e+02 1.45e+03 0 geophys. misfits: 15.2 (target 20.0 [True]) | smallness misfit: 341.1 (target: 100.0 [False]) mref changed in 0 places Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells Beta cooling evaluation: progress: [15.2]; minimum progress targets: [22.] Warming alpha_pgi to favor clustering: 953.3968727578674 5 3.16e-02 1.52e+01 3.40e+03 1.23e+02 3.53e+03 0 geophys. misfits: 20.3 (target 20.0 [False]) | smallness misfit: 321.9 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.3]; minimum progress targets: [22.] 6 3.16e-02 2.03e+01 3.22e+03 1.22e+02 4.49e+03 0 Skip BFGS geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.] 7 3.16e-02 2.02e+01 3.22e+03 1.22e+02 4.43e+03 0 geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.] 8 3.16e-02 2.02e+01 3.22e+03 1.22e+02 4.44e+03 0 geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.] 9 3.16e-02 2.02e+01 3.22e+03 1.22e+02 4.43e+03 0 geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.] 10 3.16e-02 2.02e+01 3.22e+03 1.22e+02 4.43e+03 0 geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.] 11 3.16e-02 2.02e+01 3.22e+03 1.22e+02 4.43e+03 0 geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.] 12 3.16e-02 2.02e+01 3.22e+03 1.22e+02 4.43e+03 0 geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.] 13 3.16e-02 2.02e+01 3.22e+03 1.22e+02 4.43e+03 0 geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.] 14 3.16e-02 2.07e+01 3.18e+03 1.21e+02 7.96e+00 0 Skip BFGS geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.] 15 3.16e-02 2.07e+01 3.18e+03 1.21e+02 6.78e+00 0 Skip BFGS geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.] 16 3.16e-02 2.07e+01 3.18e+03 1.21e+02 1.18e+00 0 geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.] 17 3.16e-02 2.07e+01 3.18e+03 1.21e+02 9.20e-01 0 Skip BFGS geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.] 18 3.16e-02 2.07e+01 3.18e+03 1.21e+02 6.82e-01 0 geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.] 19 3.16e-02 2.07e+01 3.18e+03 1.21e+02 5.17e-01 0 Skip BFGS geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False]) mref changed in 0 places Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.] 20 3.16e-02 2.07e+01 3.18e+03 1.21e+02 2.80e-01 0 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 3.6151e-08 <= tolF*(1+|f0|) = 2.0000e+04 1 : |xc-x_last| = 7.4356e-06 <= tolX*(1+|x0|) = 1.0000e-01 0 : |proj(x-g)-x| = 2.8026e-01 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 2.8026e-01 <= 1e3*eps = 1.0000e-02 1 : maxIter = 20 <= iter = 20 ------------------------- DONE! ------------------------- .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 12.770 seconds) **Estimated memory usage:** 166 MB .. _sphx_glr_download_content_examples_10-pgi_plot_inv_0_PGI_Linear_1D.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_inv_0_PGI_Linear_1D.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_inv_0_PGI_Linear_1D.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_inv_0_PGI_Linear_1D.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_