Note
Go to the end to download the full example code.
Petrophysically guided inversion: Joint linear example with nonlinear relationships#
We do a comparison between the classic least-squares inversion and our formulation of a petrophysically guided inversion. We explore it through coupling two linear problems whose respective physical properties are linked by polynomial relationships that change between rock units.
Running inversion with SimPEG v0.25.0
Alpha scales: [np.float64(3.4024880074053807), np.float64(0.0), np.float64(3.87793496498065e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09420521 0.90579479]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09420521 0.90579479]
================================================= 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.95e+01 3.00e+05 0.00e+00 3.00e+05 0 inf inf
1 1.95e+01 1.26e+03 1.70e+02 4.59e+03 1.41e+02 0 21 8.51e-04 7.83e+03
geophys. misfits: 7515.7 (target 30.0 [False]); 609.2 (target 30.0 [False]) | smallness misfit: 3977.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [7515.7 609.2]; minimum progress targets: [240000. 240000.]
2 1.95e+01 6.37e+01 3.98e+01 8.41e+02 1.39e+02 0 100 1.88e-02 1.47e+02 Skip BFGS
geophys. misfits: 519.4 (target 30.0 [False]); 16.3 (target 30.0 [True]) | smallness misfit: 1426.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [519.4 16.3]; minimum progress targets: [6012.6 487.4]
Updating scaling for data misfits by 1.8371604159321522
New scales: [0.16041869 0.83958131]
3 1.95e+01 5.50e+01 4.05e+01 8.46e+02 9.71e+01 0 100 1.92e-01 8.30e+01 Skip BFGS
geophys. misfits: 258.0 (target 30.0 [False]); 16.2 (target 30.0 [True]) | smallness misfit: 1170.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [258. 16.2]; minimum progress targets: [415.5 30. ]
Updating scaling for data misfits by 1.8516231258157523
New scales: [0.26133266 0.73866734]
4 1.95e+01 4.88e+01 4.17e+01 8.63e+02 7.88e+01 0 100 4.78e-02 2.11e+01
geophys. misfits: 140.1 (target 30.0 [False]); 16.4 (target 30.0 [True]) | smallness misfit: 1106.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [140.1 16.4]; minimum progress targets: [206.4 30. ]
Updating scaling for data misfits by 1.8247688867322356
New scales: [0.39231295 0.60768705]
5 1.95e+01 4.57e+01 4.25e+01 8.76e+02 7.86e+01 0 100 1.01e-01 3.56e+01 Skip BFGS
geophys. misfits: 89.6 (target 30.0 [False]); 17.4 (target 30.0 [True]) | smallness misfit: 1052.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [89.6 17.4]; minimum progress targets: [112.1 30. ]
Updating scaling for data misfits by 1.7209301546770703
New scales: [0.52629192 0.47370808]
6 1.95e+01 4.49e+01 4.29e+01 8.84e+02 7.48e+01 0 100 1.18e-01 2.99e+01 Skip BFGS
geophys. misfits: 67.9 (target 30.0 [False]); 19.3 (target 30.0 [True]) | smallness misfit: 1004.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [67.9 19.3]; minimum progress targets: [71.7 30. ]
Updating scaling for data misfits by 1.5534358396625516
New scales: [0.63314526 0.36685474]
7 1.95e+01 4.50e+01 4.32e+01 8.88e+02 7.24e+01 0 100 5.60e-02 8.87e+00 Skip BFGS
geophys. misfits: 58.3 (target 30.0 [False]); 22.0 (target 30.0 [True]) | smallness misfit: 965.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [58.3 22. ]; minimum progress targets: [54.3 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.3627451540780815
New scales: [0.70166416 0.29833584]
8 9.77e+00 2.94e+01 4.44e+01 4.63e+02 8.21e+01 0 100 1.15e-02 5.40e+00
geophys. misfits: 34.9 (target 30.0 [False]); 16.6 (target 30.0 [True]) | smallness misfit: 995.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [34.9 16.6]; minimum progress targets: [46.7 30. ]
Updating scaling for data misfits by 1.8069566093300384
New scales: [0.80951765 0.19048235]
9 9.77e+00 3.04e+01 4.45e+01 4.65e+02 5.12e+01 0 100 1.33e+00 8.55e+01 Skip BFGS
geophys. misfits: 32.7 (target 30.0 [False]); 20.4 (target 30.0 [True]) | smallness misfit: 938.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [32.7 20.4]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.4673921937762933
New scales: [0.86180531 0.13819469]
10 4.88e+00 2.37e+01 4.55e+01 2.46e+02 7.23e+01 0 100 3.03e-02 7.78e+00
geophys. misfits: 24.8 (target 30.0 [True]); 16.6 (target 30.0 [True]) | smallness misfit: 981.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [24.8 16.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.5069569346014748
11 4.88e+00 2.42e+01 4.65e+01 2.51e+02 2.54e+01 0 100 5.74e-01 1.53e+01
geophys. misfits: 24.9 (target 30.0 [True]); 19.8 (target 30.0 [True]) | smallness misfit: 904.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [24.9 19.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.0489585519210505
12 4.88e+00 2.48e+01 4.74e+01 2.57e+02 3.08e+01 0 100 6.58e-01 2.23e+01 Skip BFGS
geophys. misfits: 25.0 (target 30.0 [True]); 23.5 (target 30.0 [True]) | smallness misfit: 841.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [25. 23.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.5367240818607706
13 4.88e+00 2.54e+01 4.83e+01 2.61e+02 3.11e+01 0 100 1.77e+00 6.26e+01 Skip BFGS
geophys. misfits: 25.2 (target 30.0 [True]); 26.9 (target 30.0 [True]) | smallness misfit: 795.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [25.2 26.9]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 2.9272987408902322
14 4.88e+00 2.59e+01 4.89e+01 2.65e+02 4.77e+01 0 100 3.47e-01 2.36e+01 Skip BFGS
geophys. misfits: 25.3 (target 30.0 [True]); 29.8 (target 30.0 [True]) | smallness misfit: 763.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [25.3 29.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 3.209889180578086
15 4.88e+00 2.62e+01 4.93e+01 2.67e+02 2.68e+01 0 100 8.46e-01 2.34e+01 Skip BFGS
geophys. misfits: 25.3 (target 30.0 [True]); 31.9 (target 30.0 [False]) | smallness misfit: 742.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [25.3 31.9]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.1853682969539452
New scales: [0.84027993 0.15972007]
16 2.44e+00 2.13e+01 5.07e+01 1.45e+02 7.46e+01 0 100 1.58e-01 2.04e+01
geophys. misfits: 21.9 (target 30.0 [True]); 17.7 (target 30.0 [True]) | smallness misfit: 842.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [21.9 17.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 4.914676418226267
17 2.44e+00 2.22e+01 5.35e+01 1.53e+02 4.10e+01 0 100 9.60e-01 4.72e+01
geophys. misfits: 22.3 (target 30.0 [True]); 21.7 (target 30.0 [True]) | smallness misfit: 735.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.3 21.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 6.700152019092635
18 2.44e+00 2.34e+01 5.60e+01 1.60e+02 5.50e+01 0 100 2.18e+00 1.42e+02 Skip BFGS
geophys. misfits: 22.7 (target 30.0 [True]); 27.0 (target 30.0 [True]) | smallness misfit: 655.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.7 27. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 8.154819196460231
19 2.44e+00 2.42e+01 5.79e+01 1.65e+02 5.66e+01 0 100 4.70e-01 6.84e+01 Skip BFGS
geophys. misfits: 22.9 (target 30.0 [True]); 30.9 (target 30.0 [False]) | smallness misfit: 602.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.9 30.9]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.3102842333940516
New scales: [0.80060314 0.19939686]
20 1.22e+00 2.01e+01 6.01e+01 9.35e+01 7.54e+01 0 100 2.79e+00 3.37e+02
geophys. misfits: 20.9 (target 30.0 [True]); 17.0 (target 30.0 [True]) | smallness misfit: 711.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.9 17. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 13.050785660361534
21 1.22e+00 2.15e+01 6.63e+01 1.03e+02 6.06e+01 0 100 3.40e-01 1.16e+02
geophys. misfits: 21.7 (target 30.0 [True]); 21.1 (target 30.0 [True]) | smallness misfit: 575.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [21.7 21.1]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 18.313976601586425
22 1.22e+00 2.27e+01 7.19e+01 1.11e+02 6.36e+01 0 100 3.24e+00 4.11e+02 Skip BFGS
geophys. misfits: 22.0 (target 30.0 [True]); 25.5 (target 30.0 [True]) | smallness misfit: 471.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [22. 25.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 23.246956912835785
23 1.22e+00 2.39e+01 7.63e+01 1.17e+02 7.24e+01 0 100 1.80e-01 7.43e+01 Skip BFGS
geophys. misfits: 22.4 (target 30.0 [True]); 30.1 (target 30.0 [False]) | smallness misfit: 410.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.4 30.1]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.3394400303679477
New scales: [0.74985078 0.25014922]
24 6.10e-01 1.94e+01 8.09e+01 6.88e+01 8.00e+01 0 100 9.06e-01 1.25e+02
geophys. misfits: 20.8 (target 30.0 [True]); 15.1 (target 30.0 [True]) | smallness misfit: 498.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.8 15.1]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 39.87443360748594
25 6.10e-01 2.11e+01 9.57e+01 7.95e+01 6.60e+01 0 100 1.45e+00 2.07e+02
geophys. misfits: 21.8 (target 30.0 [True]); 19.2 (target 30.0 [True]) | smallness misfit: 379.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [21.8 19.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 58.56502837212406
26 6.10e-01 2.21e+01 1.10e+02 8.90e+01 8.19e+01 0 100 3.50e-01 7.21e+01
geophys. misfits: 22.1 (target 30.0 [True]); 22.2 (target 30.0 [True]) | smallness misfit: 307.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.1 22.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 79.26969261648453
27 6.10e-01 2.29e+01 1.21e+02 9.69e+01 6.87e+01 0 100 2.85e+02 3.26e+04
geophys. misfits: 22.9 (target 30.0 [True]); 22.8 (target 30.0 [True]) | smallness misfit: 253.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.9 22.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 104.19696849477768
28 6.10e-01 2.66e+01 1.30e+02 1.06e+02 8.84e+01 0 100 8.08e-03 2.63e+02
geophys. misfits: 23.8 (target 30.0 [True]); 34.9 (target 30.0 [False]) | smallness misfit: 235.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [23.8 34.9]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.2609509219863182
New scales: [0.70390238 0.29609762]
29 3.05e-01 2.12e+01 1.41e+02 6.41e+01 9.35e+01 0 100 2.14e+01 6.71e+03
geophys. misfits: 22.4 (target 30.0 [True]); 18.1 (target 30.0 [True]) | smallness misfit: 271.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.4 18.1]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 155.89616147158043
30 3.05e-01 2.16e+01 1.71e+02 7.38e+01 8.76e+01 0 100 8.01e-01 5.38e+03
geophys. misfits: 22.7 (target 30.0 [True]); 19.2 (target 30.0 [True]) | smallness misfit: 221.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.7 19.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 224.95586256788266
31 3.05e-01 2.31e+01 2.01e+02 8.45e+01 9.08e+01 0 100 9.74e-01 5.24e+03
geophys. misfits: 23.3 (target 30.0 [True]); 22.5 (target 30.0 [True]) | smallness misfit: 203.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [23.3 22.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 294.67029670650703
32 3.05e-01 2.35e+01 2.35e+02 9.52e+01 9.31e+01 0 100 2.33e-01 1.22e+03
geophys. misfits: 24.1 (target 30.0 [True]); 22.1 (target 30.0 [True]) | smallness misfit: 188.1 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [24.1 22.1]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 383.75969813761947
------------------------- STOP! -------------------------
1 : |fc-fOld| = 6.6976e-01 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 3.6445e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 9.3091e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 9.3091e+01 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 32
------------------------- DONE! -------------------------
Running inversion with SimPEG v0.25.0
Alpha scales: [np.float64(0.0003446943237708252), np.float64(0.0), np.float64(3.3795805153555175e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09420521 0.90579479]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09420521 0.90579479]
================================================= Projected GNCG =================================================
# beta phi_d phi_m f |proj(x-g)-x| LS iter_CG CG |Ax-b|/|b| CG |Ax-b| Comment
-----------------------------------------------------------------------------------------------------------------
0 2.17e+03 3.00e+05 0.00e+00 3.00e+05 0 inf inf
1 2.17e+03 7.26e+04 1.94e+01 1.15e+05 1.41e+02 0 12 9.40e-04 8.65e+03
geophys. misfits: 101527.3 (target 30.0 [False]); 69617.2 (target 30.0 [False]) | smallness misfit: 272.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [101527.3 69617.2]; minimum progress targets: [240000. 240000.]
2 2.17e+03 1.04e+02 5.66e-01 1.33e+03 1.37e+02 0 100 1.67e+00 4.91e+04 Skip BFGS
geophys. misfits: 876.6 (target 30.0 [False]); 23.5 (target 30.0 [True]) | smallness misfit: 105.4 (target: 200.0 [True])
Beta cooling evaluation: progress: [876.6 23.5]; minimum progress targets: [81221.9 55693.7]
Updating scaling for data misfits by 1.2745756474097902
New scales: [0.11704415 0.88295585]
3 2.17e+03 3.62e+01 1.41e-01 3.43e+02 1.15e+02 0 80 9.71e-04 4.65e+01 Skip BFGS
geophys. misfits: 162.4 (target 30.0 [False]); 19.4 (target 30.0 [True]) | smallness misfit: 59.2 (target: 200.0 [True])
Beta cooling evaluation: progress: [162.4 19.4]; minimum progress targets: [701.3 30. ]
Updating scaling for data misfits by 1.5425425184144526
New scales: [0.16976524 0.83023476]
4 2.17e+03 3.02e+01 1.42e-01 3.38e+02 8.51e+01 0 100 4.79e-01 6.25e+02
geophys. misfits: 100.6 (target 30.0 [False]); 15.8 (target 30.0 [True]) | smallness misfit: 54.7 (target: 200.0 [True])
Beta cooling evaluation: progress: [100.6 15.8]; minimum progress targets: [129.9 30. ]
Updating scaling for data misfits by 1.892914470768749
New scales: [0.27905092 0.72094908]
5 2.17e+03 2.96e+01 1.44e-01 3.41e+02 8.94e+01 0 100 4.19e-01 5.07e+02
geophys. misfits: 60.7 (target 30.0 [False]); 17.5 (target 30.0 [True]) | smallness misfit: 59.1 (target: 200.0 [True])
Beta cooling evaluation: progress: [60.7 17.5]; minimum progress targets: [80.5 30. ]
Updating scaling for data misfits by 1.710322516663376
New scales: [0.39831467 0.60168533]
6 2.17e+03 3.11e+01 1.43e-01 3.42e+02 9.53e+01 0 100 5.39e-01 3.88e+02
geophys. misfits: 47.1 (target 30.0 [False]); 20.6 (target 30.0 [True]) | smallness misfit: 54.1 (target: 200.0 [True])
Beta cooling evaluation: progress: [47.1 20.6]; minimum progress targets: [48.6 30. ]
Updating scaling for data misfits by 1.4592161516408186
New scales: [0.49135264 0.50864736]
7 2.17e+03 3.27e+01 1.43e-01 3.44e+02 8.47e+01 0 100 8.39e+00 2.93e+03 Skip BFGS
geophys. misfits: 42.0 (target 30.0 [False]); 23.7 (target 30.0 [True]) | smallness misfit: 54.2 (target: 200.0 [True])
Beta cooling evaluation: progress: [42. 23.7]; minimum progress targets: [37.6 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.2641946658421352
New scales: [0.54979499 0.45020501]
8 1.08e+03 2.53e+01 1.49e-01 1.87e+02 1.04e+02 0 100 3.46e-01 1.12e+03
geophys. misfits: 33.0 (target 30.0 [False]); 15.9 (target 30.0 [True]) | smallness misfit: 65.6 (target: 200.0 [True])
Beta cooling evaluation: progress: [33. 15.9]; minimum progress targets: [33.6 30. ]
Updating scaling for data misfits by 1.883327312650989
New scales: [0.69696407 0.30303593]
9 1.08e+03 2.80e+01 1.47e-01 1.87e+02 8.72e+01 0 100 8.38e-01 7.04e+02
geophys. misfits: 31.3 (target 30.0 [False]); 20.5 (target 30.0 [True]) | smallness misfit: 65.8 (target: 200.0 [True])
Beta cooling evaluation: progress: [31.3 20.5]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.4663055067721396
New scales: [0.77129332 0.22870668]
10 5.42e+02 2.75e+01 1.47e-01 1.07e+02 9.33e+01 0 100 2.15e+00 1.28e+03
geophys. misfits: 29.0 (target 30.0 [True]); 22.5 (target 30.0 [True]) | smallness misfit: 55.2 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [29. 22.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering: 1.1854875450938418
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.0552e+01 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 1.6569e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 9.3261e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 9.3261e+01 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 10
------------------------- DONE! -------------------------
Running inversion with SimPEG v0.25.0
Alpha scales: [np.float64(3.38037264854586e-05), np.float64(0.0), np.float64(3.409166991765572e-05), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers: [0.09420521 0.90579479]
/home/vsts/work/1/s/simpeg/directives/_directives.py:334: UserWarning:
There is no PGI regularization. Smallness target is turned off (TriggerSmall flag)
Initial data misfit scales: [0.09420521 0.90579479]
================================================= 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.05e+06 3.00e+05 0.00e+00 3.00e+05 0 inf inf
1 1.05e+06 4.13e+04 4.20e-02 8.54e+04 1.41e+02 0 23 9.24e-04 8.50e+03
geophys. misfits: 60749.7 (target 30.0 [False]); 39266.3 (target 30.0 [False])
2 2.10e+05 4.36e+03 1.06e-01 2.67e+04 1.37e+02 0 99 7.16e-04 1.32e+01 Skip BFGS
geophys. misfits: 8662.8 (target 30.0 [False]); 3909.8 (target 30.0 [False])
3 4.20e+04 2.76e+02 1.41e-01 6.19e+03 1.30e+02 0 100 1.77e-02 9.12e+01 Skip BFGS
geophys. misfits: 593.9 (target 30.0 [False]); 242.8 (target 30.0 [False])
4 8.40e+03 2.66e+01 1.51e-01 1.29e+03 1.04e+02 0 100 7.27e-03 8.58e+00 Skip BFGS
geophys. misfits: 53.5 (target 30.0 [False]); 23.8 (target 30.0 [True])
Updating scaling for data misfits by 1.2600077666096456
New scales: [0.11586138 0.88413862]
5 1.68e+03 1.15e+01 1.54e-01 2.71e+02 7.80e+01 0 100 5.68e+00 1.42e+03 Skip BFGS
geophys. misfits: 22.0 (target 30.0 [True]); 10.1 (target 30.0 [True])
All targets have been reached
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.0255e+01 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 4.0791e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 7.7979e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 7.7979e+01 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 5
------------------------- DONE! -------------------------
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:302: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "b.-" (-> marker='.'). The keyword argument will take precedence.
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:309: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "r.-" (-> marker='.'). The keyword argument will take precedence.
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:353: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "b.-" (-> marker='.'). The keyword argument will take precedence.
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:360: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "r.-" (-> marker='.'). The keyword argument will take precedence.
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:367: UserWarning:
The following kwargs were not used by contour: 'label'
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:412: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "b.-" (-> marker='.'). The keyword argument will take precedence.
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:419: UserWarning:
marker is redundantly defined by the 'marker' keyword argument and the fmt string "r.-" (-> marker='.'). The keyword argument will take precedence.
import discretize as Mesh
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
import numpy as np
from simpeg import (
data_misfit,
directives,
inverse_problem,
inversion,
maps,
optimization,
regularization,
simulation,
utils,
)
# Random seed for reproductibility
np.random.seed(1)
# Mesh
N = 100
mesh = Mesh.TensorMesh([N])
# Survey design parameters
nk = 30
jk = np.linspace(1.0, 59.0, nk)
p = -0.25
q = 0.25
# Physics
def g(k):
return np.exp(p * jk[k] * mesh.cell_centers_x) * np.cos(
np.pi * q * jk[k] * mesh.cell_centers_x
)
G = np.empty((nk, mesh.nC))
for i in range(nk):
G[i, :] = g(i)
m0 = np.zeros(mesh.nC)
m0[20:41] = np.linspace(0.0, 1.0, 21)
m0[41:57] = np.linspace(-1, 0.0, 16)
poly0 = maps.PolynomialPetroClusterMap(coeffyx=np.r_[0.0, -4.0, 4.0])
poly1 = maps.PolynomialPetroClusterMap(coeffyx=np.r_[-0.0, 3.0, 6.0, 6.0])
poly0_inverse = maps.PolynomialPetroClusterMap(coeffyx=-np.r_[0.0, -4.0, 4.0])
poly1_inverse = maps.PolynomialPetroClusterMap(coeffyx=-np.r_[0.0, 3.0, 6.0, 6.0])
cluster_mapping = [maps.IdentityMap(), poly0_inverse, poly1_inverse]
m1 = np.zeros(100)
m1[20:41] = 1.0 + (poly0 * np.vstack([m0[20:41], m1[20:41]]).T)[:, 1]
m1[41:57] = -1.0 + (poly1 * np.vstack([m0[41:57], m1[41:57]]).T)[:, 1]
model2d = np.vstack([m0, m1]).T
m = utils.mkvc(model2d)
clfmapping = utils.GaussianMixtureWithNonlinearRelationships(
mesh=mesh,
n_components=3,
covariance_type="full",
tol=1e-8,
reg_covar=1e-3,
max_iter=1000,
n_init=100,
init_params="kmeans",
random_state=None,
warm_start=False,
means_init=np.array(
[
[0, 0],
[m0[20:41].mean(), m1[20:41].mean()],
[m0[41:57].mean(), m1[41:57].mean()],
]
),
verbose=0,
verbose_interval=10,
cluster_mapping=cluster_mapping,
)
clfmapping = clfmapping.fit(model2d)
clfnomapping = utils.WeightedGaussianMixture(
mesh=mesh,
n_components=3,
covariance_type="full",
tol=1e-8,
reg_covar=1e-3,
max_iter=1000,
n_init=100,
init_params="kmeans",
random_state=None,
warm_start=False,
verbose=0,
verbose_interval=10,
)
clfnomapping = clfnomapping.fit(model2d)
wires = maps.Wires(("m1", mesh.nC), ("m2", mesh.nC))
relatrive_error = 0.01
noise_floor = 0.0
prob1 = simulation.LinearSimulation(mesh, G=G, model_map=wires.m1)
survey1 = prob1.make_synthetic_data(
m, relative_error=relatrive_error, noise_floor=noise_floor, add_noise=True
)
prob2 = simulation.LinearSimulation(mesh, G=G, model_map=wires.m2)
survey2 = prob2.make_synthetic_data(
m, relative_error=relatrive_error, noise_floor=noise_floor, add_noise=True
)
dmis1 = data_misfit.L2DataMisfit(simulation=prob1, data=survey1)
dmis2 = data_misfit.L2DataMisfit(simulation=prob2, data=survey2)
dmis = dmis1 + dmis2
minit = np.zeros_like(m)
# Distance weighting
wr1 = np.sum(prob1.G**2.0, axis=0) ** 0.5 / mesh.cell_volumes
wr1 = wr1 / np.max(wr1)
wr2 = np.sum(prob2.G**2.0, axis=0) ** 0.5 / mesh.cell_volumes
wr2 = wr2 / np.max(wr2)
reg_simple = regularization.PGI(
mesh=mesh,
gmmref=clfmapping,
gmm=clfmapping,
approx_gradient=True,
wiresmap=wires,
non_linear_relationships=True,
weights_list=[wr1, wr2],
)
opt = optimization.ProjectedGNCG(
maxIter=50,
tolX=1e-6,
cg_maxiter=100,
cg_rtol=1e-3,
lower=-10,
upper=10,
)
invProb = inverse_problem.BaseInvProblem(dmis, reg_simple, opt)
# directives
scales = directives.ScalingMultipleDataMisfits_ByEig(
chi0_ratio=np.r_[1.0, 1.0], verbose=True, n_pw_iter=10
)
scaling_schedule = directives.JointScalingSchedule(verbose=True)
alpha0_ratio = np.r_[1e6, 1e4, 1, 1]
alphas = directives.AlphasSmoothEstimate_ByEig(
alpha0_ratio=alpha0_ratio, n_pw_iter=10, verbose=True
)
beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-5, n_pw_iter=10)
betaIt = directives.PGI_BetaAlphaSchedule(
verbose=True,
coolingFactor=2.0,
progress=0.2,
)
targets = directives.MultiTargetMisfits(verbose=True)
petrodir = directives.PGI_UpdateParameters(update_gmm=False)
# Setup Inversion
inv = inversion.BaseInversion(
invProb,
directiveList=[alphas, scales, beta, petrodir, targets, betaIt, scaling_schedule],
)
mcluster_map = inv.run(minit)
# Inversion with no nonlinear mapping
reg_simple_no_map = regularization.PGI(
mesh=mesh,
gmmref=clfnomapping,
gmm=clfnomapping,
approx_gradient=True,
wiresmap=wires,
non_linear_relationships=False,
weights_list=[wr1, wr2],
)
opt = optimization.ProjectedGNCG(
maxIter=50,
tolX=1e-6,
cg_maxiter=100,
cg_rtol=1e-3,
lower=-10,
upper=10,
)
invProb = inverse_problem.BaseInvProblem(dmis, reg_simple_no_map, opt)
# directives
scales = directives.ScalingMultipleDataMisfits_ByEig(
chi0_ratio=np.r_[1.0, 1.0], verbose=True, n_pw_iter=10
)
scaling_schedule = directives.JointScalingSchedule(verbose=True)
alpha0_ratio = np.r_[100.0 * np.ones(2), 1, 1]
alphas = directives.AlphasSmoothEstimate_ByEig(
alpha0_ratio=alpha0_ratio, n_pw_iter=10, verbose=True
)
beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-5, n_pw_iter=10)
betaIt = directives.PGI_BetaAlphaSchedule(
verbose=True,
coolingFactor=2.0,
progress=0.2,
)
targets = directives.MultiTargetMisfits(
chiSmall=1.0, TriggerSmall=True, TriggerTheta=False, verbose=True
)
petrodir = directives.PGI_UpdateParameters(update_gmm=False)
# Setup Inversion
inv = inversion.BaseInversion(
invProb,
directiveList=[alphas, scales, beta, petrodir, targets, betaIt, scaling_schedule],
)
mcluster_no_map = inv.run(minit)
# WeightedLeastSquares Inversion
reg1 = regularization.WeightedLeastSquares(
mesh, alpha_s=1.0, alpha_x=1.0, mapping=wires.m1, weights={"cell_weights": wr1}
)
reg2 = regularization.WeightedLeastSquares(
mesh, alpha_s=1.0, alpha_x=1.0, mapping=wires.m2, weights={"cell_weights": wr2}
)
reg = reg1 + reg2
opt = optimization.ProjectedGNCG(
maxIter=50,
tolX=1e-6,
cg_maxiter=100,
cg_rtol=1e-3,
lower=-10,
upper=10,
)
invProb = inverse_problem.BaseInvProblem(dmis, reg, opt)
# directives
alpha0_ratio = np.r_[1, 1, 1, 1]
alphas = directives.AlphasSmoothEstimate_ByEig(
alpha0_ratio=alpha0_ratio, n_pw_iter=10, verbose=True
)
scales = directives.ScalingMultipleDataMisfits_ByEig(
chi0_ratio=np.r_[1.0, 1.0], verbose=True, n_pw_iter=10
)
scaling_schedule = directives.JointScalingSchedule(verbose=True)
beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-5, n_pw_iter=10)
beta_schedule = directives.BetaSchedule(coolingFactor=5.0, coolingRate=1)
targets = directives.MultiTargetMisfits(
TriggerSmall=False,
verbose=True,
)
# Setup Inversion
inv = inversion.BaseInversion(
invProb,
directiveList=[alphas, scales, beta, targets, beta_schedule, scaling_schedule],
)
mtik = inv.run(minit)
# Final Plot
fig, axes = plt.subplots(3, 4, figsize=(25, 15))
axes = axes.reshape(12)
left, width = 0.25, 0.5
bottom, height = 0.25, 0.5
right = left + width
top = bottom + height
axes[0].set_axis_off()
axes[0].text(
0.5 * (left + right),
0.5 * (bottom + top),
("Using true nonlinear\npetrophysical relationships"),
horizontalalignment="center",
verticalalignment="center",
fontsize=20,
color="black",
transform=axes[0].transAxes,
)
axes[1].plot(mesh.cell_centers_x, wires.m1 * mcluster_map, "b.-", ms=5, marker="v")
axes[1].plot(mesh.cell_centers_x, wires.m1 * m, "k--")
axes[1].set_title("Problem 1")
axes[1].legend(["Recovered Model", "True Model"], loc=1)
axes[1].set_xlabel("X")
axes[1].set_ylabel("Property 1")
axes[2].plot(mesh.cell_centers_x, wires.m2 * mcluster_map, "r.-", ms=5, marker="v")
axes[2].plot(mesh.cell_centers_x, wires.m2 * m, "k--")
axes[2].set_title("Problem 2")
axes[2].legend(["Recovered Model", "True Model"], loc=1)
axes[2].set_xlabel("X")
axes[2].set_ylabel("Property 2")
x, y = np.mgrid[-1:1:0.01, -4:2:0.01]
pos = np.empty(x.shape + (2,))
pos[:, :, 0] = x
pos[:, :, 1] = y
CS = axes[3].contour(
x,
y,
np.exp(clfmapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
100,
alpha=0.25,
cmap="viridis",
)
cs_proxy = mlines.Line2D([], [], label="True Petrophysical Distribution")
ps = axes[3].scatter(
wires.m1 * mcluster_map,
wires.m2 * mcluster_map,
marker="v",
label="Recovered model crossplot",
)
axes[3].set_title("Petrophysical Distribution")
axes[3].legend(handles=[cs_proxy, ps])
axes[3].set_xlabel("Property 1")
axes[3].set_ylabel("Property 2")
axes[4].set_axis_off()
axes[4].text(
0.5 * (left + right),
0.5 * (bottom + top),
("Using a pure\nGaussian distribution"),
horizontalalignment="center",
verticalalignment="center",
fontsize=20,
color="black",
transform=axes[4].transAxes,
)
axes[5].plot(mesh.cell_centers_x, wires.m1 * mcluster_no_map, "b.-", ms=5, marker="v")
axes[5].plot(mesh.cell_centers_x, wires.m1 * m, "k--")
axes[5].set_title("Problem 1")
axes[5].legend(["Recovered Model", "True Model"], loc=1)
axes[5].set_xlabel("X")
axes[5].set_ylabel("Property 1")
axes[6].plot(mesh.cell_centers_x, wires.m2 * mcluster_no_map, "r.-", ms=5, marker="v")
axes[6].plot(mesh.cell_centers_x, wires.m2 * m, "k--")
axes[6].set_title("Problem 2")
axes[6].legend(["Recovered Model", "True Model"], loc=1)
axes[6].set_xlabel("X")
axes[6].set_ylabel("Property 2")
CSF = axes[7].contour(
x,
y,
np.exp(clfmapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
100,
alpha=0.5,
label="True Petro. Distribution",
)
CS = axes[7].contour(
x,
y,
np.exp(clfnomapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
500,
cmap="viridis",
linestyles="--",
)
axes[7].scatter(
wires.m1 * mcluster_no_map,
wires.m2 * mcluster_no_map,
marker="v",
label="Recovered model crossplot",
)
cs_modeled_proxy = mlines.Line2D(
[], [], linestyle="--", label="Modeled Petro. Distribution"
)
axes[7].set_title("Petrophysical Distribution")
axes[7].legend(handles=[cs_proxy, cs_modeled_proxy, ps])
axes[7].set_xlabel("Property 1")
axes[7].set_ylabel("Property 2")
# Tikonov
axes[8].set_axis_off()
axes[8].text(
0.5 * (left + right),
0.5 * (bottom + top),
("Least-Squares\n~Using a single cluster"),
horizontalalignment="center",
verticalalignment="center",
fontsize=20,
color="black",
transform=axes[8].transAxes,
)
axes[9].plot(mesh.cell_centers_x, wires.m1 * mtik, "b.-", ms=5, marker="v")
axes[9].plot(mesh.cell_centers_x, wires.m1 * m, "k--")
axes[9].set_title("Problem 1")
axes[9].legend(["Recovered Model", "True Model"], loc=1)
axes[9].set_xlabel("X")
axes[9].set_ylabel("Property 1")
axes[10].plot(mesh.cell_centers_x, wires.m2 * mtik, "r.-", ms=5, marker="v")
axes[10].plot(mesh.cell_centers_x, wires.m2 * m, "k--")
axes[10].set_title("Problem 2")
axes[10].legend(["Recovered Model", "True Model"], loc=1)
axes[10].set_xlabel("X")
axes[10].set_ylabel("Property 2")
CS = axes[11].contour(
x,
y,
np.exp(clfmapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
100,
alpha=0.25,
cmap="viridis",
)
axes[11].scatter(wires.m1 * mtik, wires.m2 * mtik, marker="v")
axes[11].set_title("Petro Distribution")
axes[11].legend(handles=[cs_proxy, ps])
axes[11].set_xlabel("Property 1")
axes[11].set_ylabel("Property 2")
plt.subplots_adjust(wspace=0.3, hspace=0.3, top=0.85)
plt.show()
Total running time of the script: (0 minutes 36.980 seconds)
Estimated memory usage: 319 MB