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.2.dev6+gba60041a8
Alpha scales: [np.float64(3.4901411221158445), np.float64(0.0), np.float64(3.471954523665212e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09468326 0.90531674]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09468326 0.90531674]
================================================= 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.87e+01 3.00e+05 0.00e+00 3.00e+05 0 inf inf
1 1.87e+01 1.32e+03 1.72e+02 4.54e+03 1.41e+02 0 20 7.79e-04 7.11e+03
geophys. misfits: 7530.6 (target 30.0 [False]); 670.1 (target 30.0 [False]) | smallness misfit: 3955.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [7530.6 670.1]; minimum progress targets: [240000. 240000.]
2 1.87e+01 6.79e+01 4.10e+01 8.36e+02 1.39e+02 0 100 3.81e-03 2.70e+01 Skip BFGS
geophys. misfits: 529.9 (target 30.0 [False]); 19.6 (target 30.0 [True]) | smallness misfit: 1345.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [529.9 19.6]; minimum progress targets: [6024.5 536.1]
Updating scaling for data misfits by 1.5318536698930687
New scales: [0.13808715 0.86191285]
3 1.87e+01 6.11e+01 4.13e+01 8.36e+02 7.87e+01 0 100 1.24e-01 3.48e+01 Skip BFGS
geophys. misfits: 320.3 (target 30.0 [False]); 19.6 (target 30.0 [True]) | smallness misfit: 1111.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [320.3 19.6]; minimum progress targets: [423.9 30. ]
Updating scaling for data misfits by 1.5318075736117436
New scales: [0.19705226 0.80294774]
4 1.87e+01 5.47e+01 4.24e+01 8.49e+02 7.85e+01 0 100 1.69e-02 4.68e+00
geophys. misfits: 197.6 (target 30.0 [False]); 19.6 (target 30.0 [True]) | smallness misfit: 1066.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [197.6 19.6]; minimum progress targets: [256.3 30. ]
Updating scaling for data misfits by 1.5272309730932387
New scales: [0.27262114 0.72737886]
5 1.87e+01 4.96e+01 4.32e+01 8.60e+02 7.45e+01 0 100 5.35e-02 1.40e+01 Skip BFGS
geophys. misfits: 128.7 (target 30.0 [False]); 19.9 (target 30.0 [True]) | smallness misfit: 1030.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [128.7 19.9]; minimum progress targets: [158. 30.]
Updating scaling for data misfits by 1.5058173381085256
New scales: [0.36076888 0.63923112]
6 1.87e+01 4.58e+01 4.38e+01 8.67e+02 7.32e+01 0 100 1.02e-01 2.36e+01 Skip BFGS
geophys. misfits: 90.6 (target 30.0 [False]); 20.5 (target 30.0 [True]) | smallness misfit: 1000.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [90.6 20.5]; minimum progress targets: [102.9 30. ]
Updating scaling for data misfits by 1.4661170493971438
New scales: [0.45278829 0.54721171]
7 1.87e+01 4.32e+01 4.43e+01 8.73e+02 7.07e+01 0 100 2.78e-02 5.37e+00 Skip BFGS
geophys. misfits: 69.6 (target 30.0 [False]); 21.3 (target 30.0 [True]) | smallness misfit: 974.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [69.6 21.3]; minimum progress targets: [72.5 30. ]
Updating scaling for data misfits by 1.407256331818008
New scales: [0.53798438 0.46201562]
8 1.87e+01 4.15e+01 4.45e+01 8.76e+02 6.78e+01 0 100 5.49e-02 7.87e+00 Skip BFGS
geophys. misfits: 57.7 (target 30.0 [False]); 22.5 (target 30.0 [True]) | smallness misfit: 950.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [57.7 22.5]; minimum progress targets: [55.7 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.3329400717791915
New scales: [0.60816798 0.39183202]
9 9.37e+00 2.45e+01 4.59e+01 4.54e+02 8.54e+01 0 100 1.78e-02 8.62e+00
geophys. misfits: 28.4 (target 30.0 [True]); 18.5 (target 30.0 [True]) | smallness misfit: 980.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [28.4 18.5]; minimum progress targets: [46.2 30. ]
Warming alpha_pgi to favor clustering: 1.3396859140244435
10 9.37e+00 2.50e+01 4.66e+01 4.61e+02 3.04e+01 0 100 6.39e+00 2.20e+02
geophys. misfits: 28.2 (target 30.0 [True]); 19.9 (target 30.0 [True]) | smallness misfit: 923.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [28.2 19.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.7212075817153383
11 9.37e+00 2.54e+01 4.73e+01 4.69e+02 5.63e+01 0 100 5.94e-02 1.31e+01 Skip BFGS
geophys. misfits: 28.0 (target 30.0 [True]); 21.3 (target 30.0 [True]) | smallness misfit: 877.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [28. 21.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.1356284522593376
12 9.37e+00 2.60e+01 4.81e+01 4.77e+02 3.94e+01 0 100 3.56e-01 1.69e+01
geophys. misfits: 28.0 (target 30.0 [True]); 22.9 (target 30.0 [True]) | smallness misfit: 836.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [28. 22.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.544791199610468
13 9.37e+00 2.66e+01 4.88e+01 4.84e+02 3.79e+01 0 100 5.72e-01 2.78e+01 Skip BFGS
geophys. misfits: 27.9 (target 30.0 [True]); 24.6 (target 30.0 [True]) | smallness misfit: 801.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [27.9 24.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.9233398405338336
14 9.37e+00 2.72e+01 4.94e+01 4.90e+02 4.12e+01 0 100 5.04e+00 2.71e+02 Skip BFGS
geophys. misfits: 27.9 (target 30.0 [True]); 26.2 (target 30.0 [True]) | smallness misfit: 772.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [27.9 26.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3.250338939211129
15 9.37e+00 2.77e+01 5.00e+01 4.96e+02 6.22e+01 0 100 5.41e-01 1.48e+02 Skip BFGS
geophys. misfits: 27.8 (target 30.0 [True]); 27.6 (target 30.0 [True]) | smallness misfit: 749.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [27.8 27.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3.5206517639211206
16 9.37e+00 2.82e+01 5.04e+01 5.00e+02 5.70e+01 0 100 3.37e-01 5.16e+01 Skip BFGS
geophys. misfits: 27.8 (target 30.0 [True]); 28.8 (target 30.0 [True]) | smallness misfit: 731.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [27.8 28.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3.730622545539475
17 9.37e+00 2.86e+01 5.07e+01 5.04e+02 4.01e+01 0 100 8.92e-01 5.54e+01 Skip BFGS
geophys. misfits: 27.9 (target 30.0 [True]); 29.7 (target 30.0 [True]) | smallness misfit: 718.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [27.9 29.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3.8907736130804627
18 9.37e+00 2.88e+01 5.10e+01 5.06e+02 3.82e+01 0 100 9.54e-01 5.99e+01 Skip BFGS
geophys. misfits: 27.8 (target 30.0 [True]); 30.4 (target 30.0 [False]) | smallness misfit: 708.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [27.8 30.4]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.0775622185145939
New scales: [0.59023009 0.40976991]
19 4.69e+00 2.02e+01 5.22e+01 2.65e+02 9.07e+01 0 100 8.62e+00 2.13e+03
geophys. misfits: 20.0 (target 30.0 [True]); 20.4 (target 30.0 [True]) | smallness misfit: 770.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [20. 20.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 5.7764814270440645
20 4.69e+00 2.15e+01 5.52e+01 2.80e+02 8.75e+01 0 100 5.54e-02 1.18e+02
geophys. misfits: 20.2 (target 30.0 [True]); 23.3 (target 30.0 [True]) | smallness misfit: 679.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.2 23.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 8.012206885188785
21 4.69e+00 2.28e+01 5.84e+01 2.96e+02 8.30e+01 0 100 5.59e+00 9.81e+02
geophys. misfits: 20.4 (target 30.0 [True]); 26.4 (target 30.0 [True]) | smallness misfit: 606.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.4 26.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 10.459895491576585
22 4.69e+00 2.39e+01 6.16e+01 3.13e+02 8.93e+01 0 100 2.55e+00 2.54e+03 Skip BFGS
geophys. misfits: 20.6 (target 30.0 [True]); 28.5 (target 30.0 [True]) | smallness misfit: 554.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.6 28.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 13.102817517005537
23 4.69e+00 2.54e+01 6.47e+01 3.29e+02 9.82e+01 0 100 1.52e+00 3.89e+03
geophys. misfits: 21.1 (target 30.0 [True]); 31.6 (target 30.0 [False]) | smallness misfit: 502.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [21.1 31.6]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.4247967229051466
New scales: [0.50272184 0.49727816]
24 2.34e+00 1.93e+01 6.67e+01 1.75e+02 1.01e+02 0 100 2.51e-01 8.26e+02
geophys. misfits: 18.2 (target 30.0 [True]); 20.5 (target 30.0 [True]) | smallness misfit: 548.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [18.2 20.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 20.40541480469091
25 2.34e+00 2.06e+01 7.50e+01 1.96e+02 9.07e+01 0 100 7.74e+00 6.53e+03
geophys. misfits: 18.8 (target 30.0 [True]); 22.5 (target 30.0 [True]) | smallness misfit: 452.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [18.8 22.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 29.920641694060656
26 2.34e+00 2.05e+01 8.52e+01 2.20e+02 1.02e+02 0 100 1.29e+00 8.44e+03
geophys. misfits: 19.4 (target 30.0 [True]); 21.7 (target 30.0 [True]) | smallness misfit: 398.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.4 21.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 43.81277087715237
27 2.34e+00 2.23e+01 9.72e+01 2.50e+02 1.10e+02 0 100 8.64e-01 7.30e+03
geophys. misfits: 20.7 (target 30.0 [True]); 24.0 (target 30.0 [True]) | smallness misfit: 325.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.7 24. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 59.155673792343954
28 2.34e+00 2.49e+01 1.08e+02 2.79e+02 1.14e+02 0 100 1.94e+00 1.42e+04
geophys. misfits: 22.9 (target 30.0 [True]); 27.0 (target 30.0 [True]) | smallness misfit: 307.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.9 27. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 71.73593414034096
29 2.34e+00 2.85e+01 1.16e+02 3.00e+02 1.21e+02 0 100 1.66e+00 2.37e+04
geophys. misfits: 26.4 (target 30.0 [True]); 30.8 (target 30.0 [False]) | smallness misfit: 273.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.4 30.8]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.1383218085963995
New scales: [0.47036786 0.52963214]
30 1.17e+00 2.08e+01 1.20e+02 1.61e+02 1.15e+02 0 100 6.70e-01 1.48e+04
geophys. misfits: 20.1 (target 30.0 [True]); 21.5 (target 30.0 [True]) | smallness misfit: 285.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.1 21.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 103.64406029792389
31 1.17e+00 2.01e+01 1.43e+02 1.88e+02 1.11e+02 0 100 7.96e-01 1.18e+04
geophys. misfits: 19.6 (target 30.0 [True]); 20.6 (target 30.0 [True]) | smallness misfit: 268.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.6 20.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 154.659763788171
32 1.17e+00 1.97e+01 1.78e+02 2.28e+02 1.14e+02 0 100 8.09e-01 9.58e+03
geophys. misfits: 19.4 (target 30.0 [True]); 20.0 (target 30.0 [True]) | smallness misfit: 245.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.4 20. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 235.65069233235056
33 1.17e+00 1.96e+01 2.28e+02 2.86e+02 1.23e+02 1 100 5.98e-01 5.75e+03
geophys. misfits: 19.1 (target 30.0 [True]); 20.1 (target 30.0 [True]) | smallness misfit: 239.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.1 20.1]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 360.81893375848955
34 1.17e+00 1.99e+01 3.01e+02 3.73e+02 1.31e+02 1 100 1.76e+00 1.36e+04
geophys. misfits: 18.6 (target 30.0 [True]); 21.0 (target 30.0 [True]) | smallness misfit: 227.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [18.6 21. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 548.9049144914871
35 1.17e+00 1.98e+01 4.04e+02 4.93e+02 1.35e+02 1 100 2.45e+00 2.35e+04
geophys. misfits: 18.9 (target 30.0 [True]); 20.5 (target 30.0 [True]) | smallness misfit: 217.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [18.9 20.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 836.1437207472135
36 1.17e+00 2.05e+01 5.57e+02 6.73e+02 1.37e+02 2 100 6.67e-01 1.10e+04
geophys. misfits: 19.5 (target 30.0 [True]); 21.4 (target 30.0 [True]) | smallness misfit: 212.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.5 21.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1228.5568774890037
37 1.17e+00 2.13e+01 7.62e+02 9.14e+02 1.39e+02 2 100 1.16e+00 1.35e+04
geophys. misfits: 19.8 (target 30.0 [True]); 22.6 (target 30.0 [True]) | smallness misfit: 205.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.8 22.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1746.0403875811378
38 1.17e+00 2.21e+01 1.03e+03 1.23e+03 1.39e+02 2 100 1.62e+00 2.20e+04
geophys. misfits: 20.3 (target 30.0 [True]); 23.7 (target 30.0 [True]) | smallness misfit: 204.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.3 23.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2399.2915476340286
39 1.17e+00 2.23e+01 1.37e+03 1.62e+03 1.40e+02 3 100 4.18e+00 7.40e+04
geophys. misfits: 20.9 (target 30.0 [True]); 23.5 (target 30.0 [True]) | smallness misfit: 202.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.9 23.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3256.2493413100306
40 1.17e+00 2.39e+01 1.81e+03 2.15e+03 1.41e+02 2 100 1.43e+00 2.63e+04
geophys. misfits: 23.8 (target 30.0 [True]); 24.0 (target 30.0 [True]) | smallness misfit: 202.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [23.8 24. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 4090.8880556582417
41 1.17e+00 2.41e+01 2.23e+03 2.64e+03 1.42e+02 3 100 1.53e+00 3.48e+04
geophys. misfits: 24.0 (target 30.0 [True]); 24.2 (target 30.0 [True]) | smallness misfit: 201.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [24. 24.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 5083.109463253743
42 1.17e+00 2.44e+01 2.73e+03 3.22e+03 1.42e+02 3 100 1.59e+00 4.51e+04
geophys. misfits: 23.9 (target 30.0 [True]); 24.9 (target 30.0 [True]) | smallness misfit: 199.1 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [23.9 24.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 6256.64255300541
------------------------- STOP! -------------------------
1 : |fc-fOld| = 2.3964e+01 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 1.9340e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 1.4197e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 1.4197e+02 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 42
------------------------- DONE! -------------------------
Running inversion with SimPEG v0.25.2.dev6+gba60041a8
Alpha scales: [np.float64(0.00035330249402099496), np.float64(0.0), np.float64(3.5104581278710835e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09468326 0.90531674]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09468326 0.90531674]
================================================= 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.90e+03 3.00e+05 0.00e+00 3.00e+05 0 inf inf
1 1.90e+03 6.49e+04 2.31e+01 1.09e+05 1.41e+02 0 15 5.59e-04 5.10e+03
geophys. misfits: 92469.4 (target 30.0 [False]); 62017.5 (target 30.0 [False]) | smallness misfit: 243.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [92469.4 62017.5]; minimum progress targets: [240000. 240000.]
2 1.90e+03 7.46e+01 5.59e-01 1.14e+03 1.36e+02 0 100 1.02e-02 2.77e+02 Skip BFGS
geophys. misfits: 599.1 (target 30.0 [False]); 19.7 (target 30.0 [True]) | smallness misfit: 138.0 (target: 200.0 [True])
Beta cooling evaluation: progress: [599.1 19.7]; minimum progress targets: [73975.5 49614. ]
Updating scaling for data misfits by 1.5195152232733635
New scales: [0.13712744 0.86287256]
3 1.90e+03 2.78e+01 1.35e-01 2.86e+02 1.21e+02 0 100 3.05e-01 8.98e+02 Skip BFGS
geophys. misfits: 79.9 (target 30.0 [False]); 19.5 (target 30.0 [True]) | smallness misfit: 70.6 (target: 200.0 [True])
Beta cooling evaluation: progress: [79.9 19.5]; minimum progress targets: [479.3 30. ]
Updating scaling for data misfits by 1.5358525302386359
New scales: [0.19619136 0.80380864]
4 1.90e+03 2.74e+01 1.33e-01 2.81e+02 9.61e+01 0 100 1.37e-02 2.44e+01
geophys. misfits: 52.8 (target 30.0 [False]); 21.2 (target 30.0 [True]) | smallness misfit: 49.2 (target: 200.0 [True])
Beta cooling evaluation: progress: [52.8 21.2]; minimum progress targets: [63.9 30. ]
Updating scaling for data misfits by 1.4146764885322782
New scales: [0.256666 0.743334]
5 1.90e+03 2.62e+01 1.35e-01 2.83e+02 7.03e+01 0 100 1.13e+01 1.49e+03 Skip BFGS
geophys. misfits: 39.7 (target 30.0 [False]); 21.6 (target 30.0 [True]) | smallness misfit: 49.9 (target: 200.0 [True])
Beta cooling evaluation: progress: [39.7 21.6]; minimum progress targets: [42.3 30. ]
Updating scaling for data misfits by 1.3919515159087206
New scales: [0.3246106 0.6753894]
6 1.90e+03 2.56e+01 1.35e-01 2.84e+02 9.53e+01 0 100 2.23e-02 3.21e+01
geophys. misfits: 32.1 (target 30.0 [False]); 22.5 (target 30.0 [True]) | smallness misfit: 50.2 (target: 200.0 [True])
Beta cooling evaluation: progress: [32.1 22.5]; minimum progress targets: [31.8 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.331610600951955
New scales: [0.39024704 0.60975296]
7 9.52e+02 1.92e+01 1.40e-01 1.53e+02 1.00e+02 0 100 1.73e-01 5.20e+01
geophys. misfits: 20.5 (target 30.0 [True]); 18.4 (target 30.0 [True]) | smallness misfit: 52.3 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [20.5 18.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.5444702590575077
------------------------- STOP! -------------------------
1 : |fc-fOld| = 2.4605e+00 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 1.5606e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 1.0015e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 1.0015e+02 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 7
------------------------- DONE! -------------------------
Running inversion with SimPEG v0.25.2.dev6+gba60041a8
Alpha scales: [np.float64(3.151820750435301e-05), np.float64(0.0), np.float64(3.146620206562199e-05), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09468326 0.90531674]
/home/vsts/work/1/s/simpeg/directives/_directives.py:334: UserWarning: There is no PGI regularization. Smallness target is turned off (TriggerSmall flag)
getattr(r, ruleType)()
Initial data misfit scales: [0.09468326 0.90531674]
================================================= 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.07e+06 3.00e+05 0.00e+00 3.00e+05 0 inf inf
1 1.07e+06 4.21e+04 4.12e-02 8.64e+04 1.41e+02 0 23 9.46e-04 8.63e+03
geophys. misfits: 61688.4 (target 30.0 [False]); 40101.1 (target 30.0 [False])
2 2.15e+05 4.53e+03 1.06e-01 2.72e+04 1.37e+02 0 97 9.98e-04 1.86e+01 Skip BFGS
geophys. misfits: 8929.0 (target 30.0 [False]); 4070.9 (target 30.0 [False])
3 4.29e+04 2.89e+02 1.41e-01 6.33e+03 1.31e+02 0 100 9.59e-03 5.02e+01 Skip BFGS
geophys. misfits: 607.7 (target 30.0 [False]); 255.4 (target 30.0 [False])
4 8.59e+03 3.06e+01 1.51e-01 1.33e+03 1.02e+02 0 100 1.03e-01 1.24e+02 Skip BFGS
geophys. misfits: 46.2 (target 30.0 [False]); 28.9 (target 30.0 [True])
Updating scaling for data misfits by 1.0365739111280552
New scales: [0.09780749 0.90219251]
5 1.72e+03 1.42e+01 1.54e-01 2.79e+02 9.06e+01 0 100 7.77e-01 2.18e+02 Skip BFGS
geophys. misfits: 18.4 (target 30.0 [True]); 13.7 (target 30.0 [True])
All targets have been reached
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.0585e+01 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 4.6209e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 9.0578e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 9.0578e+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.
axes[1].plot(mesh.cell_centers_x, wires.m1 * mcluster_map, "b.-", ms=5, marker="v")
/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.
axes[2].plot(mesh.cell_centers_x, wires.m2 * mcluster_map, "r.-", ms=5, marker="v")
/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.
axes[5].plot(mesh.cell_centers_x, wires.m1 * mcluster_no_map, "b.-", ms=5, marker="v")
/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.
axes[6].plot(mesh.cell_centers_x, wires.m2 * mcluster_no_map, "r.-", ms=5, marker="v")
/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'
CSF = axes[7].contour(
/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.
axes[9].plot(mesh.cell_centers_x, wires.m1 * mtik, "b.-", ms=5, marker="v")
/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.
axes[10].plot(mesh.cell_centers_x, wires.m2 * mtik, "r.-", ms=5, marker="v")
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 41.816 seconds)
Estimated memory usage: 325 MB