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.1.dev27+ge1af3e81a
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.4732325724342), np.float64(0.0), np.float64(3.4873412732152422e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09674557 0.90325443]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09674557 0.90325443]
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 2.51e+01 3.00e+05 0.00e+00 3.00e+05 1.41e+02 0
geophys. misfits: 1604.8 (target 30.0 [False]); 129.6 (target 30.0 [False]) | smallness misfit: 3102.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [1604.8 129.6]; minimum progress targets: [240000. 240000.]
1 2.51e+01 2.72e+02 3.96e+01 1.27e+03 9.64e+01 0
geophys. misfits: 681.2 (target 30.0 [False]); 41.0 (target 30.0 [False]) | smallness misfit: 1471.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [681.2 41. ]; minimum progress targets: [1283.8 103.7]
2 2.51e+01 1.03e+02 3.90e+01 1.08e+03 5.81e+01 0 Skip BFGS
geophys. misfits: 648.4 (target 30.0 [False]); 40.3 (target 30.0 [False]) | smallness misfit: 1437.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [648.4 40.3]; minimum progress targets: [545. 32.8]
Decreasing beta to counter data misfit decrase plateau.
3 1.25e+01 9.91e+01 3.91e+01 5.90e+02 9.02e+01 0
geophys. misfits: 274.5 (target 30.0 [False]); 33.1 (target 30.0 [False]) | smallness misfit: 1501.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [274.5 33.1]; minimum progress targets: [518.7 32.2]
4 1.25e+01 5.64e+01 4.15e+01 5.77e+02 5.52e+01 0
geophys. misfits: 272.5 (target 30.0 [False]); 32.9 (target 30.0 [False]) | smallness misfit: 1497.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [272.5 32.9]; minimum progress targets: [219.6 30. ]
Decreasing beta to counter data misfit decrase plateau.
5 6.27e+00 5.61e+01 4.15e+01 3.16e+02 8.36e+01 0 Skip BFGS
geophys. misfits: 112.6 (target 30.0 [False]); 27.0 (target 30.0 [True]) | smallness misfit: 1774.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [112.6 27. ]; minimum progress targets: [218. 30.]
Updating scaling for data misfits by 1.1091383263554575
New scales: [0.10618307 0.89381693]
6 6.27e+00 3.61e+01 4.38e+01 3.11e+02 7.50e+01 0
geophys. misfits: 103.7 (target 30.0 [False]); 27.2 (target 30.0 [True]) | smallness misfit: 1745.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [103.7 27.2]; minimum progress targets: [90.1 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.1039453915914745
New scales: [0.11594065 0.88405935]
7 3.14e+00 3.61e+01 4.39e+01 1.74e+02 7.69e+01 0 Skip BFGS
geophys. misfits: 58.3 (target 30.0 [False]); 22.2 (target 30.0 [True]) | smallness misfit: 2147.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [58.3 22.2]; minimum progress targets: [83. 30.]
Updating scaling for data misfits by 1.3528811144913078
New scales: [0.15068874 0.84931126]
8 3.14e+00 2.76e+01 4.60e+01 1.72e+02 8.74e+01 0
geophys. misfits: 41.0 (target 30.0 [False]); 23.6 (target 30.0 [True]) | smallness misfit: 2071.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [41. 23.6]; minimum progress targets: [46.7 30. ]
Updating scaling for data misfits by 1.2700607693011543
New scales: [0.18390003 0.81609997]
9 3.14e+00 2.68e+01 4.60e+01 1.71e+02 7.12e+01 0
geophys. misfits: 33.0 (target 30.0 [False]); 23.2 (target 30.0 [True]) | smallness misfit: 2126.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [33. 23.2]; minimum progress targets: [32.8 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.2912921516689961
New scales: [0.22539458 0.77460542]
10 1.57e+00 2.54e+01 4.65e+01 9.84e+01 7.31e+01 0 Skip BFGS
geophys. misfits: 24.4 (target 30.0 [True]); 19.8 (target 30.0 [True]) | smallness misfit: 2583.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [24.4 19.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.3715168255447403
11 1.57e+00 2.09e+01 5.01e+01 9.95e+01 5.93e+01 0
geophys. misfits: 24.4 (target 30.0 [True]); 19.9 (target 30.0 [True]) | smallness misfit: 1770.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [24.4 19.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.8790423790606863
12 1.57e+00 2.09e+01 5.01e+01 9.95e+01 6.69e+01 0
geophys. misfits: 23.4 (target 30.0 [True]); 20.6 (target 30.0 [True]) | smallness misfit: 1579.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [23.4 20.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.574551292390448
13 1.57e+00 2.12e+01 5.20e+01 1.03e+02 4.69e+01 0
geophys. misfits: 23.7 (target 30.0 [True]); 21.7 (target 30.0 [True]) | smallness misfit: 1334.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [23.7 21.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3.4150979235328194
14 1.57e+00 2.21e+01 5.36e+01 1.06e+02 6.15e+01 0
geophys. misfits: 23.9 (target 30.0 [True]); 22.8 (target 30.0 [True]) | smallness misfit: 1141.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [23.9 22.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 4.3878966767236545
15 1.57e+00 2.31e+01 5.51e+01 1.09e+02 6.40e+01 0
geophys. misfits: 24.0 (target 30.0 [True]); 24.0 (target 30.0 [True]) | smallness misfit: 998.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [24. 24.]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 5.4819497322587125
16 1.57e+00 2.40e+01 5.66e+01 1.13e+02 7.22e+01 0
geophys. misfits: 24.2 (target 30.0 [True]); 24.8 (target 30.0 [True]) | smallness misfit: 835.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [24.2 24.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 6.706018143679931
17 1.57e+00 2.47e+01 5.76e+01 1.15e+02 8.87e+01 0
geophys. misfits: 24.5 (target 30.0 [True]); 25.9 (target 30.0 [True]) | smallness misfit: 708.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [24.5 25.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 7.993041105984709
18 1.57e+00 2.56e+01 5.82e+01 1.17e+02 9.02e+01 0
geophys. misfits: 25.0 (target 30.0 [True]); 25.7 (target 30.0 [True]) | smallness misfit: 644.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [25. 25.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 9.465397841366356
19 1.57e+00 2.55e+01 5.99e+01 1.19e+02 8.32e+01 0 Skip BFGS
geophys. misfits: 25.3 (target 30.0 [True]); 26.3 (target 30.0 [True]) | smallness misfit: 573.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [25.3 26.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 11.025262115573797
20 1.57e+00 2.60e+01 6.12e+01 1.22e+02 9.13e+01 0
geophys. misfits: 24.8 (target 30.0 [True]); 25.5 (target 30.0 [True]) | smallness misfit: 492.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [24.8 25.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 13.148297439265052
21 1.57e+00 2.54e+01 6.34e+01 1.25e+02 8.96e+01 0
geophys. misfits: 25.5 (target 30.0 [True]); 26.1 (target 30.0 [True]) | smallness misfit: 431.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [25.5 26.1]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 15.284762373267132
22 1.57e+00 2.60e+01 6.46e+01 1.27e+02 8.76e+01 0
geophys. misfits: 26.8 (target 30.0 [True]); 26.3 (target 30.0 [True]) | smallness misfit: 389.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.8 26.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 17.287182280043
23 1.57e+00 2.64e+01 6.59e+01 1.30e+02 8.07e+01 0 Skip BFGS
geophys. misfits: 27.0 (target 30.0 [True]); 26.8 (target 30.0 [True]) | smallness misfit: 359.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [27. 26.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 19.253922652703928
24 1.57e+00 2.69e+01 6.70e+01 1.32e+02 9.49e+01 0
geophys. misfits: 27.3 (target 30.0 [True]); 27.2 (target 30.0 [True]) | smallness misfit: 337.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [27.3 27.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 21.198791736796903
25 1.57e+00 2.72e+01 6.82e+01 1.34e+02 7.70e+01 0 Skip BFGS
geophys. misfits: 28.2 (target 30.0 [True]); 27.3 (target 30.0 [True]) | smallness misfit: 317.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [28.2 27.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 22.9244372012882
26 1.57e+00 2.75e+01 6.91e+01 1.36e+02 8.39e+01 0
geophys. misfits: 28.5 (target 30.0 [True]); 27.3 (target 30.0 [True]) | smallness misfit: 299.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [28.5 27.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 24.672049665599026
27 1.57e+00 2.75e+01 7.02e+01 1.38e+02 9.20e+01 0
geophys. misfits: 29.3 (target 30.0 [True]); 27.5 (target 30.0 [True]) | smallness misfit: 277.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [29.3 27.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 26.104656606462232
28 1.57e+00 2.79e+01 7.07e+01 1.39e+02 9.07e+01 0
geophys. misfits: 30.0 (target 30.0 [True]); 27.3 (target 30.0 [True]) | smallness misfit: 267.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [30. 27.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 27.408764877385455
29 1.57e+00 2.79e+01 7.14e+01 1.40e+02 8.55e+01 0 Skip BFGS
geophys. misfits: 30.5 (target 30.0 [False]); 27.4 (target 30.0 [True]) | smallness misfit: 256.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [30.5 27.4]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.0941076921213695
New scales: [0.24148374 0.75851626]
30 7.84e-01 2.82e+01 7.13e+01 8.41e+01 9.51e+01 0
geophys. misfits: 26.0 (target 30.0 [True]); 25.8 (target 30.0 [True]) | smallness misfit: 295.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [26. 25.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 31.803487663320233
31 7.84e-01 2.58e+01 7.62e+01 8.56e+01 9.13e+01 0
geophys. misfits: 26.2 (target 30.0 [True]); 25.4 (target 30.0 [True]) | smallness misfit: 270.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.2 25.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 36.949415944533285
32 7.84e-01 2.56e+01 7.92e+01 8.77e+01 8.83e+01 0
geophys. misfits: 26.2 (target 30.0 [True]); 25.1 (target 30.0 [True]) | smallness misfit: 255.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.2 25.1]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 43.27495806708066
33 7.84e-01 2.53e+01 8.27e+01 9.02e+01 8.99e+01 0
geophys. misfits: 26.1 (target 30.0 [True]); 24.6 (target 30.0 [True]) | smallness misfit: 239.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.1 24.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 51.24387351187036
34 7.84e-01 2.50e+01 8.70e+01 9.32e+01 9.45e+01 0
geophys. misfits: 25.8 (target 30.0 [True]); 23.9 (target 30.0 [True]) | smallness misfit: 218.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [25.8 23.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 61.84113317462196
35 7.84e-01 2.44e+01 9.21e+01 9.67e+01 1.02e+02 0
geophys. misfits: 25.6 (target 30.0 [True]); 23.6 (target 30.0 [True]) | smallness misfit: 206.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [25.6 23.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 75.62682308365137
36 7.84e-01 2.41e+01 9.84e+01 1.01e+02 1.05e+02 0
geophys. misfits: 25.5 (target 30.0 [True]); 23.3 (target 30.0 [True]) | smallness misfit: 193.0 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [25.5 23.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 93.03694468843317
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 1.6714e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 1.0512e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 1.0512e+02 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 37
------------------------- DONE! -------------------------
Running inversion with SimPEG v0.24.1.dev27+ge1af3e81a
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.00035005933966265846), np.float64(0.0), np.float64(3.5002561904909406e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09674557 0.90325443]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09674557 0.90325443]
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.93e+03 3.00e+05 0.00e+00 3.00e+05 1.41e+02 0
geophys. misfits: 87550.5 (target 30.0 [False]); 62675.8 (target 30.0 [False]) | smallness misfit: 273.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [87550.5 62675.8]; minimum progress targets: [240000. 240000.]
1 1.93e+03 6.51e+04 6.62e-01 6.64e+04 9.27e+01 0
geophys. misfits: 624.2 (target 30.0 [False]); 27.9 (target 30.0 [True]) | smallness misfit: 101.2 (target: 200.0 [True])
Beta cooling evaluation: progress: [624.2 27.9]; minimum progress targets: [70040.4 50140.6]
Updating scaling for data misfits by 1.076979442178532
New scales: [0.10342276 0.89657724]
2 1.93e+03 8.95e+01 2.81e-01 6.34e+02 8.11e+01 0 Skip BFGS
geophys. misfits: 135.4 (target 30.0 [False]); 24.5 (target 30.0 [True]) | smallness misfit: 60.4 (target: 200.0 [True])
Beta cooling evaluation: progress: [135.4 24.5]; minimum progress targets: [499.4 30. ]
Updating scaling for data misfits by 1.222571131456211
New scales: [0.12359662 0.87640338]
3 1.93e+03 3.82e+01 1.52e-01 3.32e+02 8.91e+01 0 Skip BFGS
geophys. misfits: 100.0 (target 30.0 [False]); 24.3 (target 30.0 [True]) | smallness misfit: 46.6 (target: 200.0 [True])
Beta cooling evaluation: progress: [100. 24.3]; minimum progress targets: [108.3 30. ]
Updating scaling for data misfits by 1.2364718737042402
New scales: [0.14848398 0.85151602]
4 1.93e+03 3.55e+01 1.26e-01 2.80e+02 7.57e+01 0
geophys. misfits: 79.4 (target 30.0 [False]); 24.3 (target 30.0 [True]) | smallness misfit: 47.3 (target: 200.0 [True])
Beta cooling evaluation: progress: [79.4 24.3]; minimum progress targets: [80. 30.]
Updating scaling for data misfits by 1.2345165869470793
New scales: [0.17713766 0.82286234]
5 1.93e+03 3.41e+01 1.28e-01 2.81e+02 5.57e+01 0 Skip BFGS
geophys. misfits: 65.2 (target 30.0 [False]); 24.6 (target 30.0 [True]) | smallness misfit: 47.8 (target: 200.0 [True])
Beta cooling evaluation: progress: [65.2 24.6]; minimum progress targets: [63.5 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.2206111666570432
New scales: [0.20808456 0.79191544]
6 9.67e+02 3.30e+01 1.29e-01 1.58e+02 9.52e+01 0 Skip BFGS
geophys. misfits: 33.8 (target 30.0 [False]); 22.6 (target 30.0 [True]) | smallness misfit: 50.4 (target: 200.0 [True])
Beta cooling evaluation: progress: [33.8 22.6]; minimum progress targets: [52.2 30. ]
Updating scaling for data misfits by 1.324677844118581
New scales: [0.25820084 0.74179916]
7 9.67e+02 2.55e+01 1.34e-01 1.55e+02 4.89e+01 0
geophys. misfits: 30.5 (target 30.0 [False]); 22.8 (target 30.0 [True]) | smallness misfit: 50.8 (target: 200.0 [True])
Beta cooling evaluation: progress: [30.5 22.8]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.3147163052041795
New scales: [0.31394933 0.68605067]
8 4.84e+02 2.52e+01 1.35e-01 9.04e+01 7.78e+01 0 Skip BFGS
geophys. misfits: 24.9 (target 30.0 [True]); 21.8 (target 30.0 [True]) | smallness misfit: 52.4 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [24.9 21.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.289682010661877
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 9.1210e-02 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 7.7812e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 7.7812e+01 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 9
------------------------- DONE! -------------------------
Running inversion with SimPEG v0.24.1.dev27+ge1af3e81a
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.50701058124783e-05), np.float64(0.0), np.float64(3.5283698871926646e-05), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09674557 0.90325443]
/home/vsts/work/1/s/simpeg/directives/_directives.py:333: UserWarning:
There is no PGI regularization. Smallness target is turned off (TriggerSmall flag)
Initial data misfit scales: [0.09674557 0.90325443]
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.04e+06 3.00e+05 0.00e+00 3.00e+05 1.41e+02 0
geophys. misfits: 55963.5 (target 30.0 [False]); 35936.3 (target 30.0 [False])
1 2.09e+05 3.79e+04 4.26e-02 4.68e+04 1.38e+02 0
geophys. misfits: 8280.2 (target 30.0 [False]); 3903.2 (target 30.0 [False])
2 4.18e+04 4.33e+03 1.07e-01 8.78e+03 1.31e+02 0 Skip BFGS
geophys. misfits: 556.2 (target 30.0 [False]); 272.2 (target 30.0 [False])
3 8.35e+03 3.00e+02 1.41e-01 1.47e+03 1.03e+02 0 Skip BFGS
geophys. misfits: 49.9 (target 30.0 [False]); 44.3 (target 30.0 [False])
4 1.67e+03 4.49e+01 1.51e-01 2.98e+02 7.58e+01 0 Skip BFGS
geophys. misfits: 25.9 (target 30.0 [True]); 22.6 (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| = 6.5113e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 7.5752e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 7.5752e+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 40.539 seconds)
Estimated memory usage: 290 MB