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.24.0
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.***
Alpha scales: [np.float64(3.7245395620699626), np.float64(0.0), np.float64(3.5181247488491804e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09910251 0.90089749]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09910251 0.90089749]
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.82e+01 3.00e+05 0.00e+00 3.00e+05 1.41e+02 0
geophys. misfits: 999.7 (target 30.0 [False]); 66.7 (target 30.0 [False]) | smallness misfit: 2936.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [999.7 66.7]; minimum progress targets: [240000. 240000.]
1 1.82e+01 1.59e+02 4.38e+01 9.57e+02 8.92e+01 0
geophys. misfits: 487.8 (target 30.0 [False]); 21.4 (target 30.0 [True]) | smallness misfit: 1266.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [487.8 21.4]; minimum progress targets: [799.8 53.4]
Updating scaling for data misfits by 1.4001094228596667
New scales: [0.13346233 0.86653767]
2 1.82e+01 8.37e+01 4.24e+01 8.55e+02 8.79e+01 0 Skip BFGS
geophys. misfits: 320.4 (target 30.0 [False]); 21.4 (target 30.0 [True]) | smallness misfit: 1183.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [320.4 21.4]; minimum progress targets: [390.3 30. ]
Updating scaling for data misfits by 1.4016226429816727
New scales: [0.17754704 0.82245296]
3 1.82e+01 7.45e+01 4.34e+01 8.65e+02 7.04e+01 0 Skip BFGS
geophys. misfits: 216.0 (target 30.0 [False]); 21.4 (target 30.0 [True]) | smallness misfit: 1128.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [216. 21.4]; minimum progress targets: [256.3 30. ]
Updating scaling for data misfits by 1.3988861033803235
New scales: [0.23194173 0.76805827]
4 1.82e+01 6.66e+01 4.43e+01 8.73e+02 8.47e+01 0 Skip BFGS
geophys. misfits: 150.8 (target 30.0 [False]); 21.6 (target 30.0 [True]) | smallness misfit: 1089.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [150.8 21.6]; minimum progress targets: [172.8 30. ]
Updating scaling for data misfits by 1.390564146421494
New scales: [0.29573938 0.70426062]
5 1.82e+01 5.98e+01 4.50e+01 8.79e+02 8.02e+01 0 Skip BFGS
geophys. misfits: 110.0 (target 30.0 [False]); 21.8 (target 30.0 [True]) | smallness misfit: 1059.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [110. 21.8]; minimum progress targets: [120.6 30. ]
Updating scaling for data misfits by 1.3769468798169848
New scales: [0.36637469 0.63362531]
6 1.82e+01 5.41e+01 4.56e+01 8.84e+02 6.88e+01 0 Skip BFGS
geophys. misfits: 84.5 (target 30.0 [False]); 22.1 (target 30.0 [True]) | smallness misfit: 1035.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [84.5 22.1]; minimum progress targets: [88. 30.]
Updating scaling for data misfits by 1.3558145684402008
New scales: [0.43944894 0.56055106]
7 1.82e+01 4.95e+01 4.61e+01 8.87e+02 6.52e+01 0 Skip BFGS
geophys. misfits: 68.4 (target 30.0 [False]); 22.6 (target 30.0 [True]) | smallness misfit: 1016.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [68.4 22.6]; minimum progress targets: [67.6 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.3254531235868914
New scales: [0.50958772 0.49041228]
8 9.10e+00 4.59e+01 4.64e+01 4.68e+02 8.52e+01 0 Skip BFGS
geophys. misfits: 31.0 (target 30.0 [False]); 20.8 (target 30.0 [True]) | smallness misfit: 1014.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [31. 20.8]; minimum progress targets: [54.7 30. ]
Updating scaling for data misfits by 1.443718196889594
New scales: [0.60002696 0.39997304]
9 9.10e+00 2.69e+01 4.78e+01 4.61e+02 5.01e+01 0
geophys. misfits: 28.7 (target 30.0 [True]); 21.5 (target 30.0 [True]) | smallness misfit: 992.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [28.7 21.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.220955943392176
10 9.10e+00 2.58e+01 4.84e+01 4.66e+02 3.02e+01 0 Skip BFGS
geophys. misfits: 28.2 (target 30.0 [True]); 21.9 (target 30.0 [True]) | smallness misfit: 963.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [28.2 21.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.4866590916037858
11 9.10e+00 2.57e+01 4.90e+01 4.71e+02 4.89e+01 0
geophys. misfits: 27.7 (target 30.0 [True]); 22.5 (target 30.0 [True]) | smallness misfit: 926.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [27.7 22.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.7942062206553802
12 9.10e+00 2.56e+01 4.96e+01 4.77e+02 3.28e+01 0
geophys. misfits: 27.2 (target 30.0 [True]); 23.6 (target 30.0 [True]) | smallness misfit: 890.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [27.2 23.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.131996055411159
13 9.10e+00 2.57e+01 5.03e+01 4.83e+02 3.20e+01 0
geophys. misfits: 27.0 (target 30.0 [True]); 24.3 (target 30.0 [True]) | smallness misfit: 851.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [27. 24.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.497671087995785
14 9.10e+00 2.59e+01 5.09e+01 4.89e+02 4.46e+01 0
geophys. misfits: 26.7 (target 30.0 [True]); 25.6 (target 30.0 [True]) | smallness misfit: 813.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.7 25.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.8678124640105236
15 9.10e+00 2.63e+01 5.15e+01 4.95e+02 7.46e+01 0
geophys. misfits: 26.4 (target 30.0 [True]); 27.2 (target 30.0 [True]) | smallness misfit: 778.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.4 27.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3.2115742207052396
16 9.10e+00 2.67e+01 5.21e+01 5.00e+02 5.48e+01 0
geophys. misfits: 26.4 (target 30.0 [True]); 28.2 (target 30.0 [True]) | smallness misfit: 746.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.4 28.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3.53232874149047
17 9.10e+00 2.71e+01 5.26e+01 5.05e+02 6.00e+01 0 Skip BFGS
geophys. misfits: 26.4 (target 30.0 [True]); 29.5 (target 30.0 [True]) | smallness misfit: 719.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.4 29.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3.8034716900889873
18 9.10e+00 2.76e+01 5.29e+01 5.09e+02 5.61e+01 0 Skip BFGS
geophys. misfits: 26.5 (target 30.0 [True]); 30.6 (target 30.0 [False]) | smallness misfit: 699.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.5 30.6]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.1316094788539754
New scales: [0.57002092 0.42997908]
19 4.55e+00 2.83e+01 5.29e+01 2.69e+02 9.29e+01 0 Skip BFGS
geophys. misfits: 21.0 (target 30.0 [True]); 23.2 (target 30.0 [True]) | smallness misfit: 745.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [21. 23.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 5.174039965987582
20 4.55e+00 2.20e+01 5.61e+01 2.77e+02 7.94e+01 0
geophys. misfits: 20.7 (target 30.0 [True]); 25.2 (target 30.0 [True]) | smallness misfit: 674.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.7 25.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 6.825530703675283
21 4.55e+00 2.26e+01 5.83e+01 2.88e+02 7.21e+01 0
geophys. misfits: 20.9 (target 30.0 [True]); 27.1 (target 30.0 [True]) | smallness misfit: 594.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.9 27.1]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 8.680776045318122
22 4.55e+00 2.36e+01 6.05e+01 2.99e+02 8.22e+01 0
geophys. misfits: 21.0 (target 30.0 [True]); 29.4 (target 30.0 [True]) | smallness misfit: 536.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [21. 29.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 10.627376831235125
23 4.55e+00 2.46e+01 6.25e+01 3.09e+02 8.95e+01 0 Skip BFGS
geophys. misfits: 21.3 (target 30.0 [True]); 32.4 (target 30.0 [False]) | smallness misfit: 480.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [21.3 32.4]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.4088754406961084
New scales: [0.48479086 0.51520914]
24 2.27e+00 2.70e+01 6.20e+01 1.68e+02 9.96e+01 0
geophys. misfits: 19.6 (target 30.0 [True]); 25.3 (target 30.0 [True]) | smallness misfit: 541.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.6 25.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 14.459565932208442
25 2.27e+00 2.25e+01 6.79e+01 1.77e+02 8.17e+01 0
geophys. misfits: 19.6 (target 30.0 [True]); 23.3 (target 30.0 [True]) | smallness misfit: 473.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.6 23.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 20.36894416568931
26 2.27e+00 2.15e+01 7.40e+01 1.90e+02 9.92e+01 0
geophys. misfits: 19.3 (target 30.0 [True]); 22.6 (target 30.0 [True]) | smallness misfit: 440.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.3 22.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 29.389417331610204
27 2.27e+00 2.10e+01 8.27e+01 2.09e+02 1.06e+02 0
geophys. misfits: 19.7 (target 30.0 [True]); 21.2 (target 30.0 [True]) | smallness misfit: 378.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.7 21.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 43.171478718888295
28 2.27e+00 2.05e+01 9.36e+01 2.33e+02 1.13e+02 0
geophys. misfits: 19.8 (target 30.0 [True]); 21.0 (target 30.0 [True]) | smallness misfit: 350.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.8 21. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 63.48703479728032
29 2.27e+00 2.04e+01 1.09e+02 2.68e+02 1.19e+02 0
geophys. misfits: 21.2 (target 30.0 [True]); 21.9 (target 30.0 [True]) | smallness misfit: 295.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [21.2 21.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 88.43815116887382
30 2.27e+00 2.16e+01 1.23e+02 3.02e+02 1.21e+02 0
geophys. misfits: 22.0 (target 30.0 [True]); 22.2 (target 30.0 [True]) | smallness misfit: 285.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [22. 22.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 119.79764551216057
31 2.27e+00 2.22e+01 1.41e+02 3.43e+02 1.24e+02 1
geophys. misfits: 25.0 (target 30.0 [True]); 24.1 (target 30.0 [True]) | smallness misfit: 245.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [25. 24.1]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 146.49939524369262
32 2.27e+00 2.45e+01 1.53e+02 3.73e+02 1.28e+02 0
geophys. misfits: 25.0 (target 30.0 [True]); 30.0 (target 30.0 [True]) | smallness misfit: 219.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [25. 30.]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 161.05823908598995
33 2.27e+00 2.76e+01 1.55e+02 3.80e+02 1.26e+02 0
geophys. misfits: 28.4 (target 30.0 [True]); 36.7 (target 30.0 [False]) | smallness misfit: 195.9 (target: 200.0 [True])
Beta cooling evaluation: progress: [28.4 36.7]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.055790044178622
New scales: [0.47124559 0.52875441]
34 1.14e+00 3.28e+01 1.45e+02 1.98e+02 1.16e+02 1
geophys. misfits: 23.1 (target 30.0 [True]); 28.2 (target 30.0 [True]) | smallness misfit: 207.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [23.1 28.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 190.13082833314132
35 1.14e+00 2.58e+01 1.61e+02 2.09e+02 1.21e+02 0
geophys. misfits: 21.7 (target 30.0 [True]); 25.3 (target 30.0 [True]) | smallness misfit: 183.0 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [21.7 25.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 243.84607727100806
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 7.2887e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 1.2057e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 1.2057e+02 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 36
------------------------- DONE! -------------------------
Running inversion with SimPEG v0.24.0
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.***
Alpha scales: [np.float64(0.000347239924036788), np.float64(0.0), np.float64(3.472384035579762e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09910251 0.90089749]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09910251 0.90089749]
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.96e+03 3.00e+05 0.00e+00 3.00e+05 1.41e+02 0
geophys. misfits: 86899.5 (target 30.0 [False]); 63578.8 (target 30.0 [False]) | smallness misfit: 282.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [86899.5 63578.8]; minimum progress targets: [240000. 240000.]
1 1.96e+03 6.59e+04 6.76e-01 6.72e+04 9.16e+01 0
geophys. misfits: 619.3 (target 30.0 [False]); 31.8 (target 30.0 [False]) | smallness misfit: 88.2 (target: 200.0 [True])
Beta cooling evaluation: progress: [619.3 31.8]; minimum progress targets: [69519.6 50863. ]
2 1.96e+03 9.00e+01 2.41e-01 5.62e+02 7.28e+01 0 Skip BFGS
geophys. misfits: 124.9 (target 30.0 [False]); 33.2 (target 30.0 [False]) | smallness misfit: 54.3 (target: 200.0 [True])
Beta cooling evaluation: progress: [124.9 33.2]; minimum progress targets: [495.5 30. ]
3 1.96e+03 4.23e+01 1.54e-01 3.43e+02 7.65e+01 0 Skip BFGS
geophys. misfits: 127.5 (target 30.0 [False]); 30.7 (target 30.0 [False]) | smallness misfit: 49.6 (target: 200.0 [True])
Beta cooling evaluation: progress: [127.5 30.7]; minimum progress targets: [99.9 30. ]
Decreasing beta to counter data misfit decrase plateau.
4 9.79e+02 4.03e+01 1.31e-01 1.69e+02 9.40e+01 0
geophys. misfits: 57.8 (target 30.0 [False]); 27.3 (target 30.0 [True]) | smallness misfit: 48.0 (target: 200.0 [True])
Beta cooling evaluation: progress: [57.8 27.3]; minimum progress targets: [102. 30.]
Updating scaling for data misfits by 1.0982435093910916
New scales: [0.10778924 0.89221076]
5 9.79e+02 3.06e+01 1.31e-01 1.59e+02 3.70e+01 0
geophys. misfits: 52.7 (target 30.0 [False]); 27.5 (target 30.0 [True]) | smallness misfit: 48.2 (target: 200.0 [True])
Beta cooling evaluation: progress: [52.7 27.5]; minimum progress targets: [46.2 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.0919579202000722
New scales: [0.1165461 0.8834539]
6 4.89e+02 3.04e+01 1.31e-01 9.46e+01 8.42e+01 0 Skip BFGS
geophys. misfits: 29.8 (target 30.0 [True]); 25.2 (target 30.0 [True]) | smallness misfit: 51.0 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [29.8 25.2]; minimum progress targets: [42.1 30. ]
Warming alpha_pgi to favor clustering: 1.0990924372083506
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 2.0140e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 8.4173e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 8.4173e+01 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 7
------------------------- DONE! -------------------------
Running inversion with SimPEG v0.24.0
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: [np.float64(3.372927745135317e-05), np.float64(0.0), np.float64(3.3723126275065764e-05), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09910251 0.90089749]
/home/vsts/work/1/s/simpeg/directives/directives.py:339: UserWarning:
There is no PGI regularization. Smallness target is turned off (TriggerSmall flag)
Initial data misfit scales: [0.09910251 0.90089749]
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.07e+06 3.00e+05 0.00e+00 3.00e+05 1.41e+02 0
geophys. misfits: 56133.6 (target 30.0 [False]); 37206.4 (target 30.0 [False])
1 2.15e+05 3.91e+04 4.17e-02 4.80e+04 1.38e+02 0
geophys. misfits: 8279.7 (target 30.0 [False]); 4070.6 (target 30.0 [False])
2 4.29e+04 4.49e+03 1.06e-01 9.04e+03 1.31e+02 0 Skip BFGS
geophys. misfits: 557.0 (target 30.0 [False]); 253.2 (target 30.0 [False])
3 8.59e+03 2.83e+02 1.41e-01 1.49e+03 1.05e+02 0 Skip BFGS
geophys. misfits: 46.1 (target 30.0 [False]); 31.7 (target 30.0 [False])
4 1.72e+03 3.31e+01 1.51e-01 2.92e+02 7.68e+01 0 Skip BFGS
geophys. misfits: 18.3 (target 30.0 [True]); 20.3 (target 30.0 [True])
All targets have been reached
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 2.5962e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 7.6736e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 7.6736e+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,
maxIterCG=100,
tolCG=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,
maxIterCG=100,
tolCG=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,
maxIterCG=100,
tolCG=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 47.591 seconds)
Estimated memory usage: 289 MB