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.
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: [3.4879629923623376, 0.0, 3.4844686613163604e-06, 0.0]
Calculating the scaling parameter.
Scale Multipliers: [0.09605361 0.90394639]
<class 'SimPEG.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09605361 0.90394639]
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.86e+01 1.50e+05 0.00e+00 1.50e+05 1.41e+02 0
geophys. misfits: 498.1 (target 15.0 [False]); 35.7 (target 15.0 [False]) | smallness misfit: 1459.2 (target: 100.0 [False])
Beta cooling evaluation: progress: [498.1 35.7]; minimum progress targets: [120000. 120000.]
1 1.86e+01 8.01e+01 2.06e+01 4.65e+02 7.96e+01 0
geophys. misfits: 235.0 (target 15.0 [False]); 10.9 (target 15.0 [True]) | smallness misfit: 644.7 (target: 100.0 [False])
Beta cooling evaluation: progress: [235. 10.9]; minimum progress targets: [398.5 28.6]
Updating scaling for data misfits by 1.374269734362692
New scales: [0.12742273 0.87257727]
2 1.86e+01 3.95e+01 2.00e+01 4.11e+02 5.66e+01 0 Skip BFGS
geophys. misfits: 156.2 (target 15.0 [False]); 10.9 (target 15.0 [True]) | smallness misfit: 609.7 (target: 100.0 [False])
Beta cooling evaluation: progress: [156.2 10.9]; minimum progress targets: [188. 15.]
Updating scaling for data misfits by 1.381664440666124
New scales: [0.16789048 0.83210952]
3 1.86e+01 3.53e+01 2.04e+01 4.15e+02 5.64e+01 0 Skip BFGS
geophys. misfits: 105.3 (target 15.0 [False]); 10.9 (target 15.0 [True]) | smallness misfit: 586.1 (target: 100.0 [False])
Beta cooling evaluation: progress: [105.3 10.9]; minimum progress targets: [125. 15.]
Updating scaling for data misfits by 1.3811013876658975
New scales: [0.21792991 0.78207009]
4 1.86e+01 3.15e+01 2.08e+01 4.19e+02 5.71e+01 0
geophys. misfits: 73.0 (target 15.0 [False]); 11.0 (target 15.0 [True]) | smallness misfit: 568.0 (target: 100.0 [False])
Beta cooling evaluation: progress: [73. 11.]; minimum progress targets: [84.3 15. ]
Updating scaling for data misfits by 1.368561549136323
New scales: [0.27607591 0.72392409]
5 1.86e+01 2.81e+01 2.11e+01 4.22e+02 6.02e+01 0 Skip BFGS
geophys. misfits: 52.7 (target 15.0 [False]); 11.1 (target 15.0 [True]) | smallness misfit: 553.5 (target: 100.0 [False])
Beta cooling evaluation: progress: [52.7 11.1]; minimum progress targets: [58.4 15. ]
Updating scaling for data misfits by 1.34723481121236
New scales: [0.33940285 0.66059715]
6 1.86e+01 2.53e+01 2.14e+01 4.24e+02 5.20e+01 0 Skip BFGS
geophys. misfits: 40.2 (target 15.0 [False]); 11.4 (target 15.0 [True]) | smallness misfit: 541.1 (target: 100.0 [False])
Beta cooling evaluation: progress: [40.2 11.4]; minimum progress targets: [42.2 15. ]
Updating scaling for data misfits by 1.3174153484248725
New scales: [0.40364876 0.59635124]
7 1.86e+01 2.30e+01 2.16e+01 4.25e+02 4.71e+01 0 Skip BFGS
geophys. misfits: 32.4 (target 15.0 [False]); 11.7 (target 15.0 [True]) | smallness misfit: 530.2 (target: 100.0 [False])
Beta cooling evaluation: progress: [32.4 11.7]; minimum progress targets: [32.2 15. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.27930152260321
New scales: [0.46406942 0.53593058]
8 9.31e+00 2.13e+01 2.17e+01 2.24e+02 7.29e+01 0 Skip BFGS
geophys. misfits: 14.1 (target 15.0 [True]); 9.7 (target 15.0 [True]) | smallness misfit: 556.7 (target: 100.0 [False])
Beta cooling evaluation: progress: [14.1 9.7]; minimum progress targets: [25.9 15. ]
Warming alpha_pgi to favor clustering: 1.3060218468359308
9 9.31e+00 1.17e+01 2.28e+01 2.24e+02 2.11e+01 0
geophys. misfits: 13.9 (target 15.0 [True]); 10.3 (target 15.0 [True]) | smallness misfit: 515.5 (target: 100.0 [False])
Beta cooling evaluation: progress: [13.9 10.3]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 1.6509719957323006
10 9.31e+00 1.20e+01 2.31e+01 2.27e+02 1.78e+01 0
geophys. misfits: 13.8 (target 15.0 [True]); 10.9 (target 15.0 [True]) | smallness misfit: 481.9 (target: 100.0 [False])
Beta cooling evaluation: progress: [13.8 10.9]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 2.0352960333864827
11 9.31e+00 1.22e+01 2.35e+01 2.31e+02 3.67e+01 0 Skip BFGS
geophys. misfits: 13.7 (target 15.0 [True]); 11.5 (target 15.0 [True]) | smallness misfit: 451.6 (target: 100.0 [False])
Beta cooling evaluation: progress: [13.7 11.5]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 2.440521801517239
12 9.31e+00 1.25e+01 2.38e+01 2.34e+02 1.95e+01 0
geophys. misfits: 13.6 (target 15.0 [True]); 12.2 (target 15.0 [True]) | smallness misfit: 424.9 (target: 100.0 [False])
Beta cooling evaluation: progress: [13.6 12.2]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 2.848206268636033
13 9.31e+00 1.28e+01 2.42e+01 2.38e+02 1.93e+01 0 Skip BFGS
geophys. misfits: 13.6 (target 15.0 [True]); 12.8 (target 15.0 [True]) | smallness misfit: 402.0 (target: 100.0 [False])
Beta cooling evaluation: progress: [13.6 12.8]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 3.2390019898577744
14 9.31e+00 1.32e+01 2.45e+01 2.41e+02 2.84e+01 0 Skip BFGS
geophys. misfits: 13.6 (target 15.0 [True]); 13.5 (target 15.0 [True]) | smallness misfit: 382.8 (target: 100.0 [False])
Beta cooling evaluation: progress: [13.6 13.5]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 3.5973567509010813
15 9.31e+00 1.35e+01 2.48e+01 2.44e+02 4.30e+01 0 Skip BFGS
geophys. misfits: 13.6 (target 15.0 [True]); 14.0 (target 15.0 [True]) | smallness misfit: 367.0 (target: 100.0 [False])
Beta cooling evaluation: progress: [13.6 14. ]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 3.9132521705024024
16 9.31e+00 1.38e+01 2.50e+01 2.46e+02 2.11e+01 0 Skip BFGS
geophys. misfits: 13.6 (target 15.0 [True]); 14.5 (target 15.0 [True]) | smallness misfit: 354.5 (target: 100.0 [False])
Beta cooling evaluation: progress: [13.6 14.5]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 4.182578893798094
17 9.31e+00 1.41e+01 2.52e+01 2.48e+02 1.22e+01 0 Skip BFGS
geophys. misfits: 13.6 (target 15.0 [True]); 14.9 (target 15.0 [True]) | smallness misfit: 344.5 (target: 100.0 [False])
Beta cooling evaluation: progress: [13.6 14.9]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 4.404910734903503
18 9.31e+00 1.43e+01 2.53e+01 2.50e+02 1.02e+01 0 Skip BFGS
geophys. misfits: 13.6 (target 15.0 [True]); 15.3 (target 15.0 [False]) | smallness misfit: 336.8 (target: 100.0 [False])
Beta cooling evaluation: progress: [13.6 15.3]; minimum progress targets: [15. 15.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.1003745588507636
New scales: [0.44037973 0.55962027]
19 4.65e+00 1.46e+01 2.53e+01 1.32e+02 7.38e+01 0 Skip BFGS
geophys. misfits: 9.9 (target 15.0 [True]); 10.9 (target 15.0 [True]) | smallness misfit: 369.3 (target: 100.0 [False])
Beta cooling evaluation: progress: [ 9.9 10.9]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 6.387706037636578
20 4.65e+00 1.04e+01 2.75e+01 1.38e+02 3.69e+01 0
geophys. misfits: 10.0 (target 15.0 [True]); 11.7 (target 15.0 [True]) | smallness misfit: 313.3 (target: 100.0 [False])
Beta cooling evaluation: progress: [10. 11.7]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 8.867527960511257
21 4.65e+00 1.10e+01 2.90e+01 1.46e+02 6.39e+01 0
geophys. misfits: 10.3 (target 15.0 [True]); 12.8 (target 15.0 [True]) | smallness misfit: 267.8 (target: 100.0 [False])
Beta cooling evaluation: progress: [10.3 12.8]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 11.61834393430031
22 4.65e+00 1.17e+01 3.03e+01 1.53e+02 8.26e+01 0
geophys. misfits: 10.5 (target 15.0 [True]); 13.1 (target 15.0 [True]) | smallness misfit: 241.4 (target: 100.0 [False])
Beta cooling evaluation: progress: [10.5 13.1]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 14.926843025096527
23 4.65e+00 1.20e+01 3.20e+01 1.61e+02 7.42e+01 0
geophys. misfits: 11.2 (target 15.0 [True]); 15.1 (target 15.0 [False]) | smallness misfit: 199.9 (target: 100.0 [False])
Beta cooling evaluation: progress: [11.2 15.1]; minimum progress targets: [15. 15.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by 1.3442902520887736
New scales: [0.36923798 0.63076202]
24 2.33e+00 1.37e+01 3.15e+01 8.70e+01 7.61e+01 0
geophys. misfits: 9.9 (target 15.0 [True]); 11.0 (target 15.0 [True]) | smallness misfit: 224.8 (target: 100.0 [False])
Beta cooling evaluation: progress: [ 9.9 11. ]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 21.50041230040881
25 2.33e+00 1.06e+01 3.57e+01 9.36e+01 8.13e+01 0
geophys. misfits: 10.3 (target 15.0 [True]); 11.2 (target 15.0 [True]) | smallness misfit: 185.3 (target: 100.0 [False])
Beta cooling evaluation: progress: [10.3 11.2]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 30.12426037525881
26 2.33e+00 1.08e+01 3.88e+01 1.01e+02 9.53e+01 0
geophys. misfits: 11.6 (target 15.0 [True]); 10.6 (target 15.0 [True]) | smallness misfit: 170.0 (target: 100.0 [False])
Beta cooling evaluation: progress: [11.6 10.6]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 40.864544419405185
27 2.33e+00 1.09e+01 4.25e+01 1.10e+02 9.70e+01 0
geophys. misfits: 11.3 (target 15.0 [True]); 10.1 (target 15.0 [True]) | smallness misfit: 144.0 (target: 100.0 [False])
Beta cooling evaluation: progress: [11.3 10.1]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 57.7058247954592
28 2.33e+00 1.05e+01 4.73e+01 1.21e+02 1.04e+02 0
geophys. misfits: 12.6 (target 15.0 [True]); 11.4 (target 15.0 [True]) | smallness misfit: 116.0 (target: 100.0 [False])
Beta cooling evaluation: progress: [12.6 11.4]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 72.30312813814554
29 2.33e+00 1.18e+01 4.95e+01 1.27e+02 1.04e+02 0
geophys. misfits: 18.3 (target 15.0 [False]); 15.5 (target 15.0 [False]) | smallness misfit: 98.6 (target: 100.0 [True])
Beta cooling evaluation: progress: [18.3 15.5]; minimum progress targets: [15. 15.]
Decreasing beta to counter data misfit increase.
30 1.16e+00 1.66e+01 4.66e+01 7.07e+01 9.94e+01 0
geophys. misfits: 12.3 (target 15.0 [True]); 10.9 (target 15.0 [True]) | smallness misfit: 115.4 (target: 100.0 [False])
Beta cooling evaluation: progress: [12.3 10.9]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 93.78457997050192
31 1.16e+00 1.14e+01 5.60e+01 7.66e+01 9.25e+01 0
geophys. misfits: 10.9 (target 15.0 [True]); 10.2 (target 15.0 [True]) | smallness misfit: 109.1 (target: 100.0 [False])
Beta cooling evaluation: progress: [10.9 10.2]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 133.62216327456488
32 1.16e+00 1.04e+01 6.64e+01 8.77e+01 1.03e+02 0
geophys. misfits: 11.8 (target 15.0 [True]); 12.0 (target 15.0 [True]) | smallness misfit: 91.0 (target: 100.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [11.8 12. ]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 168.4889550423769
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.5000e+04
0 : |xc-x_last| = 4.3858e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 1.0274e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 1.0274e+02 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 33
------------------------- DONE! -------------------------
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: [0.0003484547914829826, 0.0, 3.564944499171368e-06, 0.0]
Calculating the scaling parameter.
Scale Multipliers: [0.09605361 0.90394639]
<class 'SimPEG.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.09605361 0.90394639]
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.97e+03 1.50e+05 0.00e+00 1.50e+05 1.41e+02 0
geophys. misfits: 44567.6 (target 15.0 [False]); 31701.9 (target 15.0 [False]) | smallness misfit: 145.4 (target: 100.0 [False])
Beta cooling evaluation: progress: [44567.6 31701.9]; minimum progress targets: [120000. 120000.]
1 1.97e+03 3.29e+04 3.41e-01 3.36e+04 8.98e+01 0
geophys. misfits: 329.7 (target 15.0 [False]); 11.9 (target 15.0 [True]) | smallness misfit: 53.1 (target: 100.0 [True])
Beta cooling evaluation: progress: [329.7 11.9]; minimum progress targets: [35654.1 25361.5]
Updating scaling for data misfits by 1.2588248712215981
New scales: [0.11798153 0.88201847]
2 1.97e+03 4.94e+01 1.53e-01 3.52e+02 8.08e+01 0 Skip BFGS
geophys. misfits: 50.3 (target 15.0 [False]); 12.8 (target 15.0 [True]) | smallness misfit: 23.3 (target: 100.0 [True])
Beta cooling evaluation: progress: [50.3 12.8]; minimum progress targets: [263.7 15. ]
Updating scaling for data misfits by 1.1684921475591683
New scales: [0.13517339 0.86482661]
3 1.97e+03 1.79e+01 6.52e-02 1.47e+02 7.91e+01 0 Skip BFGS
geophys. misfits: 42.5 (target 15.0 [False]); 10.5 (target 15.0 [True]) | smallness misfit: 23.0 (target: 100.0 [True])
Beta cooling evaluation: progress: [42.5 10.5]; minimum progress targets: [40.3 15. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by 1.4221529080495827
New scales: [0.18185962 0.81814038]
4 9.87e+02 1.64e+01 6.31e-02 7.87e+01 8.36e+01 0
geophys. misfits: 15.4 (target 15.0 [False]); 9.0 (target 15.0 [True]) | smallness misfit: 24.6 (target: 100.0 [True])
Beta cooling evaluation: progress: [15.4 9. ]; minimum progress targets: [34. 15.]
Updating scaling for data misfits by 1.663521297584011
New scales: [0.26995276 0.73004724]
5 9.87e+02 1.07e+01 6.68e-02 7.67e+01 6.49e+01 0
geophys. misfits: 12.0 (target 15.0 [True]); 9.4 (target 15.0 [True]) | smallness misfit: 24.8 (target: 100.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [12. 9.4]; minimum progress targets: [15. 15.]
Warming alpha_pgi to favor clustering: 1.4238182762757154
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.5000e+04
0 : |xc-x_last| = 5.7836e-02 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 6.4921e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 6.4921e+01 <= 1e3*eps = 1.0000e-02
0 : maxIter = 50 <= iter = 6
------------------------- DONE! -------------------------
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: [3.511758442591592e-05, 0.0, 3.5141198769325145e-05, 0.0]
Calculating the scaling parameter.
Scale Multipliers: [0.09605361 0.90394639]
/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.09605361 0.90394639]
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.02e+06 1.50e+05 0.00e+00 1.50e+05 1.41e+02 0
geophys. misfits: 27709.3 (target 15.0 [False]); 17519.4 (target 15.0 [False])
1 2.05e+05 1.85e+04 2.16e-02 2.29e+04 1.37e+02 0
geophys. misfits: 4050.6 (target 15.0 [False]); 1865.0 (target 15.0 [False])
2 4.10e+04 2.07e+03 5.34e-02 4.26e+03 1.24e+02 0 Skip BFGS
geophys. misfits: 268.2 (target 15.0 [False]); 118.9 (target 15.0 [False])
3 8.20e+03 1.33e+02 7.01e-02 7.08e+02 9.19e+01 0 Skip BFGS
geophys. misfits: 20.6 (target 15.0 [False]); 15.5 (target 15.0 [False])
4 1.64e+03 1.60e+01 7.51e-02 1.39e+02 6.28e+01 0 Skip BFGS
geophys. misfits: 9.2 (target 15.0 [True]); 8.1 (target 15.0 [True])
All targets have been reached
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.5000e+04
0 : |xc-x_last| = 4.3431e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x| = 6.2729e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 6.2729e+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:301: 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:308: 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:346: 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:353: 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:360: 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:368: 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:402: 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:409: 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 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
)
reg1.cell_weights = wr1
reg2 = regularization.WeightedLeastSquares(
mesh, alpha_s=1.0, alpha_x=1.0, mapping=wires.m2
)
reg2.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",
)
axes[3].scatter(wires.m1 * mcluster_map, wires.m2 * mcluster_map, marker="v")
axes[3].set_title("Petrophysical Distribution")
CS.collections[0].set_label("")
axes[3].legend(["True Petrophysical Distribution", "Recovered model crossplot"])
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="--",
label="Modeled Petro. Distribution",
)
axes[7].scatter(
wires.m1 * mcluster_no_map,
wires.m2 * mcluster_no_map,
marker="v",
label="Recovered model crossplot",
)
axes[7].set_title("Petrophysical Distribution")
axes[7].legend()
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")
CS.collections[0].set_label("")
axes[11].legend(["True Petro Distribution", "Recovered model crossplot"])
axes[11].set_xlabel("Property 1")
axes[11].set_ylabel("Property 2")
plt.subplots_adjust(wspace=0.3, hspace=0.3, top=0.85)
plt.show()
Total running time of the script: (0 minutes 47.091 seconds)
Estimated memory usage: 8 MB