Note
Go to the end to download the full example code.
Petrophysically guided inversion: Joint linear example with nonlinear relationships#
We do a comparison between the classic least-squares inversion and our formulation of a petrophysically guided inversion. We explore it through coupling two linear problems whose respective physical properties are linked by polynomial relationships that change between rock units.
Running inversion with SimPEG v0.25.0
Alpha scales: [np.float64(3.4259267988149684), np.float64(0.0), np.float64(3.4252313743201377e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09582987 0.90417013]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09582987 0.90417013]
================================================= Projected GNCG =================================================
# beta phi_d phi_m f |proj(x-g)-x| LS iter_CG CG |Ax-b|/|b| CG |Ax-b| Comment
-----------------------------------------------------------------------------------------------------------------
0 2.05e+01 3.00e+05 0.00e+00 3.00e+05 0 inf inf
1 2.05e+01 1.49e+03 1.70e+02 4.97e+03 1.40e+02 0 20 6.78e-04 6.21e+03
geophys. misfits: 8011.7 (target 30.0 [False]); 794.5 (target 30.0 [False]) | smallness misfit: 3997.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [8011.7 794.5]; minimum progress targets: [240000. 240000.]
2 2.05e+01 8.67e+01 4.02e+01 9.08e+02 1.40e+02 0 100 3.77e-03 2.34e+01 Skip BFGS
geophys. misfits: 512.3 (target 30.0 [False]); 41.6 (target 30.0 [False]) | smallness misfit: 1470.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [512.3 41.6]; minimum progress targets: [6409.3 635.6]
3 2.05e+01 8.64e+01 3.92e+01 8.89e+02 5.55e+01 0 100 2.98e-02 4.96e+00
geophys. misfits: 513.5 (target 30.0 [False]); 41.1 (target 30.0 [False]) | smallness misfit: 1300.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [513.5 41.1]; minimum progress targets: [409.9 33.3]
Decreasing beta to counter data misfit decrase plateau.
4 1.02e+01 5.19e+01 4.16e+01 4.77e+02 8.44e+01 0 100 4.60e-01 1.44e+02
geophys. misfits: 211.5 (target 30.0 [False]); 35.0 (target 30.0 [False]) | smallness misfit: 1330.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [211.5 35. ]; minimum progress targets: [410.8 32.9]
5 1.02e+01 5.17e+01 4.16e+01 4.77e+02 8.34e+01 0 100 1.08e+00 1.54e+02 Skip BFGS
geophys. misfits: 209.7 (target 30.0 [False]); 34.9 (target 30.0 [False]) | smallness misfit: 1322.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [209.7 34.9]; minimum progress targets: [169.2 30. ]
Decreasing beta to counter data misfit decrase plateau.
6 5.11e+00 3.69e+01 4.36e+01 2.60e+02 9.50e+01 0 100 7.95e-01 2.08e+02
geophys. misfits: 91.2 (target 30.0 [False]); 31.2 (target 30.0 [False]) | smallness misfit: 1438.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [91.2 31.2]; minimum progress targets: [167.8 30. ]
7 5.11e+00 3.68e+01 4.36e+01 2.60e+02 9.99e+01 0 100 9.91e-01 2.06e+02 Skip BFGS
geophys. misfits: 91.1 (target 30.0 [False]); 31.0 (target 30.0 [False]) | smallness misfit: 1444.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [91.1 31. ]; minimum progress targets: [72.9 30. ]
Decreasing beta to counter data misfit decrase plateau.
8 2.56e+00 2.97e+01 4.55e+01 1.46e+02 1.03e+02 0 100 9.04e-01 2.28e+02
geophys. misfits: 47.7 (target 30.0 [False]); 27.8 (target 30.0 [True]) | smallness misfit: 1751.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [47.7 27.8]; minimum progress targets: [72.9 30. ]
Updating scaling for data misfits by 1.079127780838916
New scales: [0.10263442 0.89736558]
9 2.56e+00 2.98e+01 4.55e+01 1.46e+02 1.01e+02 0 100 9.45e-01 2.23e+02
geophys. misfits: 44.1 (target 30.0 [False]); 28.2 (target 30.0 [True]) | smallness misfit: 1686.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [44.1 28.2]; minimum progress targets: [38.1 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.0649362896771621
New scales: [0.10857549 0.89142451]
10 1.28e+00 2.50e+01 4.81e+01 8.65e+01 1.04e+02 0 100 7.20e-01 1.83e+02
geophys. misfits: 25.5 (target 30.0 [True]); 24.9 (target 30.0 [True]) | smallness misfit: 2422.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [25.5 24.9]; minimum progress targets: [35.3 30. ]
Warming alpha_pgi to favor clustering: 1.1898192028056946
11 1.28e+00 2.51e+01 4.88e+01 8.75e+01 9.20e+01 0 100 1.17e+00 2.12e+02
geophys. misfits: 22.1 (target 30.0 [True]); 25.5 (target 30.0 [True]) | smallness misfit: 2097.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.1 25.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.5079395136510758
12 1.28e+00 2.59e+01 4.95e+01 8.92e+01 9.60e+01 0 100 4.46e+00 9.40e+02
geophys. misfits: 23.0 (target 30.0 [True]); 26.3 (target 30.0 [True]) | smallness misfit: 1812.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [23. 26.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.8461347994485606
13 1.28e+00 2.60e+01 5.07e+01 9.08e+01 1.05e+02 0 100 2.17e-01 2.04e+02
geophys. misfits: 21.7 (target 30.0 [True]); 26.5 (target 30.0 [True]) | smallness misfit: 1676.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [21.7 26.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.3218741298075396
14 1.28e+00 2.72e+01 5.13e+01 9.28e+01 1.01e+02 0 100 1.10e+00 2.22e+02
geophys. misfits: 23.5 (target 30.0 [True]); 27.7 (target 30.0 [True]) | smallness misfit: 1405.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [23.5 27.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.737172339635033
15 1.28e+00 2.78e+01 5.20e+01 9.43e+01 1.05e+02 0 100 1.64e+00 3.60e+02 Skip BFGS
geophys. misfits: 24.2 (target 30.0 [True]); 28.3 (target 30.0 [True]) | smallness misfit: 1277.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [24.2 28.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3.1514649870664284
16 1.28e+00 2.80e+01 5.30e+01 9.57e+01 1.08e+02 0 100 1.46e+00 5.25e+02
geophys. misfits: 24.5 (target 30.0 [True]); 28.4 (target 30.0 [True]) | smallness misfit: 1188.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [24.5 28.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3.5887384529036783
17 1.28e+00 2.60e+01 5.40e+01 9.51e+01 1.11e+02 0 100 8.69e+00 4.55e+03
geophys. misfits: 14.4 (target 30.0 [True]); 27.4 (target 30.0 [True]) | smallness misfit: 1143.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [14.4 27.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 5.688605867312924
18 1.28e+00 2.73e+01 5.74e+01 1.01e+02 9.03e+01 0 100 7.72e-01 3.52e+03
geophys. misfits: 15.6 (target 30.0 [True]); 28.7 (target 30.0 [True]) | smallness misfit: 895.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [15.6 28.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 8.455414319499214
19 1.28e+00 2.80e+01 6.15e+01 1.07e+02 9.29e+01 0 100 1.44e+00 5.07e+03
geophys. misfits: 15.7 (target 30.0 [True]); 29.5 (target 30.0 [True]) | smallness misfit: 760.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [15.7 29.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 12.374790674160009
20 1.28e+00 2.85e+01 6.68e+01 1.14e+02 9.16e+01 0 100 2.60e-01 1.32e+03
geophys. misfits: 15.8 (target 30.0 [True]); 30.0 (target 30.0 [True]) | smallness misfit: 650.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [15.8 30. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 17.93321648137953
21 1.28e+00 2.98e+01 7.28e+01 1.23e+02 8.95e+01 0 100 1.38e+00 1.82e+03
geophys. misfits: 18.5 (target 30.0 [True]); 31.2 (target 30.0 [False]) | smallness misfit: 554.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [18.5 31.2]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.6193175623673686
New scales: [0.06995506 0.93004494]
22 6.39e-01 2.70e+01 7.62e+01 7.57e+01 9.76e+01 0 100 3.09e-01 5.69e+02
geophys. misfits: 16.7 (target 30.0 [True]); 27.8 (target 30.0 [True]) | smallness misfit: 608.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [16.7 27.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 25.75444166769497
23 6.39e-01 2.69e+01 8.56e+01 8.16e+01 8.79e+01 0 100 2.42e+00 1.39e+03
geophys. misfits: 16.9 (target 30.0 [True]); 27.7 (target 30.0 [True]) | smallness misfit: 545.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [16.9 27.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 36.863232272907204
24 6.39e-01 2.64e+01 9.78e+01 8.89e+01 9.68e+01 0 100 2.82e+00 3.92e+03
geophys. misfits: 19.5 (target 30.0 [True]); 26.9 (target 30.0 [True]) | smallness misfit: 493.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.5 26.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 48.88623259680818
25 6.39e-01 2.61e+01 1.09e+02 9.55e+01 1.07e+02 0 100 9.01e-01 3.53e+03
geophys. misfits: 17.7 (target 30.0 [True]); 26.7 (target 30.0 [True]) | smallness misfit: 426.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [17.7 26.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 68.84232455912941
26 6.39e-01 2.60e+01 1.24e+02 1.05e+02 1.08e+02 0 100 2.01e+00 7.12e+03
geophys. misfits: 20.8 (target 30.0 [True]); 26.4 (target 30.0 [True]) | smallness misfit: 358.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.8 26.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 88.73248182308255
27 6.39e-01 2.74e+01 1.33e+02 1.13e+02 1.13e+02 0 100 8.59e-01 6.12e+03
geophys. misfits: 30.5 (target 30.0 [False]); 27.1 (target 30.0 [True]) | smallness misfit: 311.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [30.5 27.1]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.1054767765340727
New scales: [0.07676726 0.92323274]
28 3.20e-01 2.60e+01 1.34e+02 6.90e+01 9.55e+01 0 100 3.47e-01 2.12e+03
geophys. misfits: 17.3 (target 30.0 [True]); 26.7 (target 30.0 [True]) | smallness misfit: 296.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [17.3 26.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 126.59744987349254
29 3.20e-01 2.61e+01 1.58e+02 7.65e+01 1.02e+02 0 100 2.83e-01 6.02e+02
geophys. misfits: 20.3 (target 30.0 [True]); 26.6 (target 30.0 [True]) | smallness misfit: 278.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.3 26.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 164.86489446345593
30 3.20e-01 2.70e+01 1.78e+02 8.38e+01 9.59e+01 0 100 1.46e+00 8.91e+02
geophys. misfits: 27.2 (target 30.0 [True]); 27.0 (target 30.0 [True]) | smallness misfit: 260.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [27.2 27. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 182.5308232344503
31 3.20e-01 2.68e+01 1.87e+02 8.65e+01 9.38e+01 0 100 3.06e+00 2.74e+03
geophys. misfits: 22.3 (target 30.0 [True]); 27.2 (target 30.0 [True]) | smallness misfit: 247.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.3 27.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 223.51916594270446
32 3.20e-01 2.75e+01 2.05e+02 9.31e+01 1.02e+02 0 100 8.98e-01 2.47e+03
geophys. misfits: 23.4 (target 30.0 [True]); 27.9 (target 30.0 [True]) | smallness misfit: 230.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [23.4 27.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 263.5333482202469
33 3.20e-01 2.79e+01 2.23e+02 9.90e+01 1.03e+02 1 100 1.20e+00 2.97e+03
geophys. misfits: 25.7 (target 30.0 [True]); 28.0 (target 30.0 [True]) | smallness misfit: 217.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [25.7 28. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 294.8055749757426
34 3.20e-01 2.90e+01 2.32e+02 1.03e+02 1.04e+02 0 100 2.06e+00 5.58e+03
geophys. misfits: 32.0 (target 30.0 [False]); 28.8 (target 30.0 [True]) | smallness misfit: 204.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [32. 28.8]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.041766466772869
New scales: [0.07971796 0.92028204]
35 1.60e-01 2.69e+01 2.41e+02 6.54e+01 9.99e+01 0 100 3.84e-01 2.14e+03
geophys. misfits: 26.1 (target 30.0 [True]); 26.9 (target 30.0 [True]) | smallness misfit: 218.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.1 26.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 333.6964195863196
36 1.60e-01 2.72e+01 2.54e+02 6.77e+01 9.23e+01 0 100 3.47e-01 7.42e+02
geophys. misfits: 26.3 (target 30.0 [True]); 27.2 (target 30.0 [True]) | smallness misfit: 198.5 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [26.3 27.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 373.8959875052041
------------------------- STOP! -------------------------
1 : |fc-fOld| = 7.6914e-01 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 4.6693e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 9.2259e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 9.2259e+01 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 36
------------------------- DONE! -------------------------
Running inversion with SimPEG v0.25.0
Alpha scales: [np.float64(0.000356189101200682), np.float64(0.0), np.float64(3.910028918255035e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09582987 0.90417013]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09582987 0.90417013]
================================================= Projected GNCG =================================================
# beta phi_d phi_m f |proj(x-g)-x| LS iter_CG CG |Ax-b|/|b| CG |Ax-b| Comment
-----------------------------------------------------------------------------------------------------------------
0 1.94e+03 3.00e+05 0.00e+00 3.00e+05 0 inf inf
1 1.94e+03 6.60e+04 2.25e+01 1.10e+05 1.41e+02 0 15 3.05e-04 2.80e+03
geophys. misfits: 92582.7 (target 30.0 [False]); 63191.4 (target 30.0 [False]) | smallness misfit: 249.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [92582.7 63191.4]; minimum progress targets: [240000. 240000.]
2 1.94e+03 8.44e+01 5.55e-01 1.16e+03 1.35e+02 0 100 5.45e-03 1.47e+02 Skip BFGS
geophys. misfits: 577.0 (target 30.0 [False]); 32.1 (target 30.0 [False]) | smallness misfit: 118.0 (target: 200.0 [True])
Beta cooling evaluation: progress: [577. 32.1]; minimum progress targets: [74066.2 50553.1]
3 1.94e+03 3.88e+01 1.36e-01 3.02e+02 1.07e+02 0 100 7.84e-02 1.66e+02 Skip BFGS
geophys. misfits: 123.0 (target 30.0 [False]); 29.9 (target 30.0 [True]) | smallness misfit: 65.4 (target: 200.0 [True])
Beta cooling evaluation: progress: [123. 29.9]; minimum progress targets: [461.6 30. ]
Updating scaling for data misfits by 1.0025060216016555
New scales: [0.09604695 0.90395305]
4 1.94e+03 3.83e+01 1.27e-01 2.84e+02 6.31e+01 0 100 2.04e-01 2.93e+02
geophys. misfits: 118.7 (target 30.0 [False]); 29.8 (target 30.0 [True]) | smallness misfit: 46.2 (target: 200.0 [True])
Beta cooling evaluation: progress: [118.7 29.8]; minimum progress targets: [98.4 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.0083328033463455
New scales: [0.09676985 0.90323015]
5 9.68e+02 2.97e+01 1.33e-01 1.58e+02 9.89e+01 0 100 5.59e-02 2.19e+01
geophys. misfits: 46.9 (target 30.0 [False]); 27.8 (target 30.0 [True]) | smallness misfit: 49.0 (target: 200.0 [True])
Beta cooling evaluation: progress: [46.9 27.8]; minimum progress targets: [94.9 30. ]
Updating scaling for data misfits by 1.0775234008566383
New scales: [0.10349536 0.89650464]
6 9.68e+02 2.94e+01 1.33e-01 1.58e+02 2.74e+01 0 100 5.79e-01 1.59e+01 Skip BFGS
geophys. misfits: 43.0 (target 30.0 [False]); 27.8 (target 30.0 [True]) | smallness misfit: 49.2 (target: 200.0 [True])
Beta cooling evaluation: progress: [43. 27.8]; minimum progress targets: [37.5 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.07771424005532
New scales: [0.11064847 0.88935153]
7 4.84e+02 2.60e+01 1.38e-01 9.28e+01 8.25e+01 0 100 8.22e-01 1.04e+02
geophys. misfits: 20.4 (target 30.0 [True]); 26.7 (target 30.0 [True]) | smallness misfit: 51.4 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [20.4 26.7]; minimum progress targets: [34.4 30. ]
Warming alpha_pgi to favor clustering: 1.2955903477503043
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.1864e+00 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 1.9601e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 8.2516e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 8.2516e+01 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 7
------------------------- DONE! -------------------------
Running inversion with SimPEG v0.25.0
Alpha scales: [np.float64(3.5070151654117766e-05), np.float64(0.0), np.float64(3.5058668796646294e-05), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09582987 0.90417013]
/home/vsts/work/1/s/simpeg/directives/_directives.py:334: UserWarning:
There is no PGI regularization. Smallness target is turned off (TriggerSmall flag)
Initial data misfit scales: [0.09582987 0.90417013]
================================================= Projected GNCG =================================================
# beta phi_d phi_m f |proj(x-g)-x| LS iter_CG CG |Ax-b|/|b| CG |Ax-b| Comment
-----------------------------------------------------------------------------------------------------------------
0 1.04e+06 3.00e+05 0.00e+00 3.00e+05 0 inf inf
1 1.04e+06 3.88e+04 4.27e-02 8.33e+04 1.41e+02 0 33 6.98e-04 6.39e+03
geophys. misfits: 57761.8 (target 30.0 [False]); 36826.8 (target 30.0 [False])
2 2.08e+05 4.39e+03 1.07e-01 2.67e+04 1.38e+02 0 100 6.53e-04 1.14e+01 Skip BFGS
geophys. misfits: 8430.6 (target 30.0 [False]); 3966.7 (target 30.0 [False])
3 4.17e+04 3.06e+02 1.42e-01 6.22e+03 1.30e+02 0 100 8.09e-02 4.14e+02 Skip BFGS
geophys. misfits: 562.7 (target 30.0 [False]); 279.2 (target 30.0 [False])
4 8.33e+03 4.66e+01 1.53e-01 1.32e+03 1.04e+02 0 100 3.98e-01 4.92e+02 Skip BFGS
geophys. misfits: 39.0 (target 30.0 [False]); 47.4 (target 30.0 [False])
5 1.67e+03 2.75e+01 1.57e-01 2.89e+02 8.68e+01 0 100 3.63e-02 1.99e+01 Skip BFGS
geophys. misfits: 12.1 (target 30.0 [True]); 29.1 (target 30.0 [True])
All targets have been reached
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.2242e+01 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 5.1761e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 8.6763e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 8.6763e+01 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 5
------------------------- DONE! -------------------------
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:302: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "b.-" (-> marker='.'). The keyword argument will take precedence.
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:309: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "r.-" (-> marker='.'). The keyword argument will take precedence.
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:353: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "b.-" (-> marker='.'). The keyword argument will take precedence.
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:360: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "r.-" (-> marker='.'). The keyword argument will take precedence.
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:367: UserWarning:
The following kwargs were not used by contour: 'label'
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:412: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "b.-" (-> marker='.'). The keyword argument will take precedence.
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:419: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "r.-" (-> marker='.'). The keyword argument will take precedence.
import discretize as Mesh
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
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 = 30
jk = np.linspace(1.0, 59.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)
m0 = np.zeros(mesh.nC)
m0[20:41] = np.linspace(0.0, 1.0, 21)
m0[41:57] = np.linspace(-1, 0.0, 16)
poly0 = maps.PolynomialPetroClusterMap(coeffyx=np.r_[0.0, -4.0, 4.0])
poly1 = maps.PolynomialPetroClusterMap(coeffyx=np.r_[-0.0, 3.0, 6.0, 6.0])
poly0_inverse = maps.PolynomialPetroClusterMap(coeffyx=-np.r_[0.0, -4.0, 4.0])
poly1_inverse = maps.PolynomialPetroClusterMap(coeffyx=-np.r_[0.0, 3.0, 6.0, 6.0])
cluster_mapping = [maps.IdentityMap(), poly0_inverse, poly1_inverse]
m1 = np.zeros(100)
m1[20:41] = 1.0 + (poly0 * np.vstack([m0[20:41], m1[20:41]]).T)[:, 1]
m1[41:57] = -1.0 + (poly1 * np.vstack([m0[41:57], m1[41:57]]).T)[:, 1]
model2d = np.vstack([m0, m1]).T
m = utils.mkvc(model2d)
clfmapping = utils.GaussianMixtureWithNonlinearRelationships(
mesh=mesh,
n_components=3,
covariance_type="full",
tol=1e-8,
reg_covar=1e-3,
max_iter=1000,
n_init=100,
init_params="kmeans",
random_state=None,
warm_start=False,
means_init=np.array(
[
[0, 0],
[m0[20:41].mean(), m1[20:41].mean()],
[m0[41:57].mean(), m1[41:57].mean()],
]
),
verbose=0,
verbose_interval=10,
cluster_mapping=cluster_mapping,
)
clfmapping = clfmapping.fit(model2d)
clfnomapping = utils.WeightedGaussianMixture(
mesh=mesh,
n_components=3,
covariance_type="full",
tol=1e-8,
reg_covar=1e-3,
max_iter=1000,
n_init=100,
init_params="kmeans",
random_state=None,
warm_start=False,
verbose=0,
verbose_interval=10,
)
clfnomapping = clfnomapping.fit(model2d)
wires = maps.Wires(("m1", mesh.nC), ("m2", mesh.nC))
relatrive_error = 0.01
noise_floor = 0.0
prob1 = simulation.LinearSimulation(mesh, G=G, model_map=wires.m1)
survey1 = prob1.make_synthetic_data(
m, relative_error=relatrive_error, noise_floor=noise_floor, add_noise=True
)
prob2 = simulation.LinearSimulation(mesh, G=G, model_map=wires.m2)
survey2 = prob2.make_synthetic_data(
m, relative_error=relatrive_error, noise_floor=noise_floor, add_noise=True
)
dmis1 = data_misfit.L2DataMisfit(simulation=prob1, data=survey1)
dmis2 = data_misfit.L2DataMisfit(simulation=prob2, data=survey2)
dmis = dmis1 + dmis2
minit = np.zeros_like(m)
# Distance weighting
wr1 = np.sum(prob1.G**2.0, axis=0) ** 0.5 / mesh.cell_volumes
wr1 = wr1 / np.max(wr1)
wr2 = np.sum(prob2.G**2.0, axis=0) ** 0.5 / mesh.cell_volumes
wr2 = wr2 / np.max(wr2)
reg_simple = regularization.PGI(
mesh=mesh,
gmmref=clfmapping,
gmm=clfmapping,
approx_gradient=True,
wiresmap=wires,
non_linear_relationships=True,
weights_list=[wr1, wr2],
)
opt = optimization.ProjectedGNCG(
maxIter=50,
tolX=1e-6,
cg_maxiter=100,
cg_rtol=1e-3,
lower=-10,
upper=10,
)
invProb = inverse_problem.BaseInvProblem(dmis, reg_simple, opt)
# directives
scales = directives.ScalingMultipleDataMisfits_ByEig(
chi0_ratio=np.r_[1.0, 1.0], verbose=True, n_pw_iter=10
)
scaling_schedule = directives.JointScalingSchedule(verbose=True)
alpha0_ratio = np.r_[1e6, 1e4, 1, 1]
alphas = directives.AlphasSmoothEstimate_ByEig(
alpha0_ratio=alpha0_ratio, n_pw_iter=10, verbose=True
)
beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-5, n_pw_iter=10)
betaIt = directives.PGI_BetaAlphaSchedule(
verbose=True,
coolingFactor=2.0,
progress=0.2,
)
targets = directives.MultiTargetMisfits(verbose=True)
petrodir = directives.PGI_UpdateParameters(update_gmm=False)
# Setup Inversion
inv = inversion.BaseInversion(
invProb,
directiveList=[alphas, scales, beta, petrodir, targets, betaIt, scaling_schedule],
)
mcluster_map = inv.run(minit)
# Inversion with no nonlinear mapping
reg_simple_no_map = regularization.PGI(
mesh=mesh,
gmmref=clfnomapping,
gmm=clfnomapping,
approx_gradient=True,
wiresmap=wires,
non_linear_relationships=False,
weights_list=[wr1, wr2],
)
opt = optimization.ProjectedGNCG(
maxIter=50,
tolX=1e-6,
cg_maxiter=100,
cg_rtol=1e-3,
lower=-10,
upper=10,
)
invProb = inverse_problem.BaseInvProblem(dmis, reg_simple_no_map, opt)
# directives
scales = directives.ScalingMultipleDataMisfits_ByEig(
chi0_ratio=np.r_[1.0, 1.0], verbose=True, n_pw_iter=10
)
scaling_schedule = directives.JointScalingSchedule(verbose=True)
alpha0_ratio = np.r_[100.0 * np.ones(2), 1, 1]
alphas = directives.AlphasSmoothEstimate_ByEig(
alpha0_ratio=alpha0_ratio, n_pw_iter=10, verbose=True
)
beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-5, n_pw_iter=10)
betaIt = directives.PGI_BetaAlphaSchedule(
verbose=True,
coolingFactor=2.0,
progress=0.2,
)
targets = directives.MultiTargetMisfits(
chiSmall=1.0, TriggerSmall=True, TriggerTheta=False, verbose=True
)
petrodir = directives.PGI_UpdateParameters(update_gmm=False)
# Setup Inversion
inv = inversion.BaseInversion(
invProb,
directiveList=[alphas, scales, beta, petrodir, targets, betaIt, scaling_schedule],
)
mcluster_no_map = inv.run(minit)
# WeightedLeastSquares Inversion
reg1 = regularization.WeightedLeastSquares(
mesh, alpha_s=1.0, alpha_x=1.0, mapping=wires.m1, weights={"cell_weights": wr1}
)
reg2 = regularization.WeightedLeastSquares(
mesh, alpha_s=1.0, alpha_x=1.0, mapping=wires.m2, weights={"cell_weights": wr2}
)
reg = reg1 + reg2
opt = optimization.ProjectedGNCG(
maxIter=50,
tolX=1e-6,
cg_maxiter=100,
cg_rtol=1e-3,
lower=-10,
upper=10,
)
invProb = inverse_problem.BaseInvProblem(dmis, reg, opt)
# directives
alpha0_ratio = np.r_[1, 1, 1, 1]
alphas = directives.AlphasSmoothEstimate_ByEig(
alpha0_ratio=alpha0_ratio, n_pw_iter=10, verbose=True
)
scales = directives.ScalingMultipleDataMisfits_ByEig(
chi0_ratio=np.r_[1.0, 1.0], verbose=True, n_pw_iter=10
)
scaling_schedule = directives.JointScalingSchedule(verbose=True)
beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-5, n_pw_iter=10)
beta_schedule = directives.BetaSchedule(coolingFactor=5.0, coolingRate=1)
targets = directives.MultiTargetMisfits(
TriggerSmall=False,
verbose=True,
)
# Setup Inversion
inv = inversion.BaseInversion(
invProb,
directiveList=[alphas, scales, beta, targets, beta_schedule, scaling_schedule],
)
mtik = inv.run(minit)
# Final Plot
fig, axes = plt.subplots(3, 4, figsize=(25, 15))
axes = axes.reshape(12)
left, width = 0.25, 0.5
bottom, height = 0.25, 0.5
right = left + width
top = bottom + height
axes[0].set_axis_off()
axes[0].text(
0.5 * (left + right),
0.5 * (bottom + top),
("Using true nonlinear\npetrophysical relationships"),
horizontalalignment="center",
verticalalignment="center",
fontsize=20,
color="black",
transform=axes[0].transAxes,
)
axes[1].plot(mesh.cell_centers_x, wires.m1 * mcluster_map, "b.-", ms=5, marker="v")
axes[1].plot(mesh.cell_centers_x, wires.m1 * m, "k--")
axes[1].set_title("Problem 1")
axes[1].legend(["Recovered Model", "True Model"], loc=1)
axes[1].set_xlabel("X")
axes[1].set_ylabel("Property 1")
axes[2].plot(mesh.cell_centers_x, wires.m2 * mcluster_map, "r.-", ms=5, marker="v")
axes[2].plot(mesh.cell_centers_x, wires.m2 * m, "k--")
axes[2].set_title("Problem 2")
axes[2].legend(["Recovered Model", "True Model"], loc=1)
axes[2].set_xlabel("X")
axes[2].set_ylabel("Property 2")
x, y = np.mgrid[-1:1:0.01, -4:2:0.01]
pos = np.empty(x.shape + (2,))
pos[:, :, 0] = x
pos[:, :, 1] = y
CS = axes[3].contour(
x,
y,
np.exp(clfmapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
100,
alpha=0.25,
cmap="viridis",
)
cs_proxy = mlines.Line2D([], [], label="True Petrophysical Distribution")
ps = axes[3].scatter(
wires.m1 * mcluster_map,
wires.m2 * mcluster_map,
marker="v",
label="Recovered model crossplot",
)
axes[3].set_title("Petrophysical Distribution")
axes[3].legend(handles=[cs_proxy, ps])
axes[3].set_xlabel("Property 1")
axes[3].set_ylabel("Property 2")
axes[4].set_axis_off()
axes[4].text(
0.5 * (left + right),
0.5 * (bottom + top),
("Using a pure\nGaussian distribution"),
horizontalalignment="center",
verticalalignment="center",
fontsize=20,
color="black",
transform=axes[4].transAxes,
)
axes[5].plot(mesh.cell_centers_x, wires.m1 * mcluster_no_map, "b.-", ms=5, marker="v")
axes[5].plot(mesh.cell_centers_x, wires.m1 * m, "k--")
axes[5].set_title("Problem 1")
axes[5].legend(["Recovered Model", "True Model"], loc=1)
axes[5].set_xlabel("X")
axes[5].set_ylabel("Property 1")
axes[6].plot(mesh.cell_centers_x, wires.m2 * mcluster_no_map, "r.-", ms=5, marker="v")
axes[6].plot(mesh.cell_centers_x, wires.m2 * m, "k--")
axes[6].set_title("Problem 2")
axes[6].legend(["Recovered Model", "True Model"], loc=1)
axes[6].set_xlabel("X")
axes[6].set_ylabel("Property 2")
CSF = axes[7].contour(
x,
y,
np.exp(clfmapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
100,
alpha=0.5,
label="True Petro. Distribution",
)
CS = axes[7].contour(
x,
y,
np.exp(clfnomapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
500,
cmap="viridis",
linestyles="--",
)
axes[7].scatter(
wires.m1 * mcluster_no_map,
wires.m2 * mcluster_no_map,
marker="v",
label="Recovered model crossplot",
)
cs_modeled_proxy = mlines.Line2D(
[], [], linestyle="--", label="Modeled Petro. Distribution"
)
axes[7].set_title("Petrophysical Distribution")
axes[7].legend(handles=[cs_proxy, cs_modeled_proxy, ps])
axes[7].set_xlabel("Property 1")
axes[7].set_ylabel("Property 2")
# Tikonov
axes[8].set_axis_off()
axes[8].text(
0.5 * (left + right),
0.5 * (bottom + top),
("Least-Squares\n~Using a single cluster"),
horizontalalignment="center",
verticalalignment="center",
fontsize=20,
color="black",
transform=axes[8].transAxes,
)
axes[9].plot(mesh.cell_centers_x, wires.m1 * mtik, "b.-", ms=5, marker="v")
axes[9].plot(mesh.cell_centers_x, wires.m1 * m, "k--")
axes[9].set_title("Problem 1")
axes[9].legend(["Recovered Model", "True Model"], loc=1)
axes[9].set_xlabel("X")
axes[9].set_ylabel("Property 1")
axes[10].plot(mesh.cell_centers_x, wires.m2 * mtik, "r.-", ms=5, marker="v")
axes[10].plot(mesh.cell_centers_x, wires.m2 * m, "k--")
axes[10].set_title("Problem 2")
axes[10].legend(["Recovered Model", "True Model"], loc=1)
axes[10].set_xlabel("X")
axes[10].set_ylabel("Property 2")
CS = axes[11].contour(
x,
y,
np.exp(clfmapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
100,
alpha=0.25,
cmap="viridis",
)
axes[11].scatter(wires.m1 * mtik, wires.m2 * mtik, marker="v")
axes[11].set_title("Petro Distribution")
axes[11].legend(handles=[cs_proxy, ps])
axes[11].set_xlabel("Property 1")
axes[11].set_ylabel("Property 2")
plt.subplots_adjust(wspace=0.3, hspace=0.3, top=0.85)
plt.show()
Total running time of the script: (0 minutes 49.025 seconds)
Estimated memory usage: 319 MB