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.

Problem 1, Problem 2, Petrophysical Distribution, Problem 1, Problem 2, Petrophysical Distribution, Problem 1, Problem 2, Petro Distribution
/home/vsts/work/1/s/simpeg/optimization.py:1535: FutureWarning:

InexactCG.tolCG has been deprecated, please use InexactCG.cg_atol. It will be removed in version 0.26.0 of SimPEG.

/home/vsts/work/1/s/simpeg/optimization.py:1061: FutureWarning:

InexactCG.maxIterCG has been deprecated, please use InexactCG.cg_maxiter. It will be removed in version 0.26.0 of SimPEG.


Running inversion with SimPEG v0.24.1.dev37+gf5e5be547
Alpha scales: [np.float64(3.482392103611884), np.float64(0.0), np.float64(3.480814983191258e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers:  [0.09546346 0.90453654]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales:  [0.09546346 0.90453654]
================================================= 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.91e+01  3.00e+05  0.00e+00  3.00e+05                         0           inf          inf
   1  1.91e+01  1.76e+02  1.79e+02  3.59e+03    1.40e+02      0     100       1.57e-06     1.44e+01
geophys. misfits: 1059.9 (target 30.0 [False]); 82.9 (target 30.0 [False]) | smallness misfit: 2916.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [1059.9   82.9]; minimum progress targets: [240000. 240000.]
   2  1.91e+01  7.50e+01  3.99e+01  8.37e+02    7.64e+01      0     100       2.45e-02     1.84e+01   Skip BFGS
geophys. misfits: 493.9 (target 30.0 [False]); 30.8 (target 30.0 [False]) | smallness misfit: 1267.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [493.9  30.8]; minimum progress targets: [847.9  66.3]
   3  1.91e+01  7.44e+01  3.99e+01  8.37e+02    2.42e+01      0     100       1.89e-01     4.58e+00
geophys. misfits: 487.5 (target 30.0 [False]); 30.8 (target 30.0 [False]) | smallness misfit: 1254.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [487.5  30.8]; minimum progress targets: [395.1  30. ]
Decreasing beta to counter data misfit decrase plateau.
   4  9.54e+00  4.51e+01  4.21e+01  4.47e+02    8.34e+01      0     100       7.49e+00     2.29e+03
geophys. misfits: 200.7 (target 30.0 [False]); 28.6 (target 30.0 [True]) | smallness misfit: 1249.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [200.7  28.6]; minimum progress targets: [390.  30.]
Updating scaling for data misfits by  1.0481755163463566
New scales: [0.09960438 0.90039562]
   5  9.54e+00  4.45e+01  4.22e+01  4.47e+02    8.45e+01      0     100       2.27e-02     5.41e+01   Skip BFGS
geophys. misfits: 189.2 (target 30.0 [False]); 28.5 (target 30.0 [True]) | smallness misfit: 1246.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [189.2  28.5]; minimum progress targets: [160.6  30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by  1.051583822579721
New scales: [0.10420694 0.89579306]
   6  4.77e+00  3.15e+01  4.41e+01  2.42e+02    7.91e+01      0     100       1.17e+01     2.41e+03
geophys. misfits: 72.1 (target 30.0 [False]); 26.8 (target 30.0 [True]) | smallness misfit: 1363.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [72.1 26.8]; minimum progress targets: [151.4  30. ]
Updating scaling for data misfits by  1.1204465149359941
New scales: [0.11531099 0.88468901]
   7  4.77e+00  3.13e+01  4.42e+01  2.42e+02    7.68e+01      0     100       9.95e-02     2.37e+02   Skip BFGS
geophys. misfits: 67.5 (target 30.0 [False]); 26.6 (target 30.0 [True]) | smallness misfit: 1366.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [67.5 26.6]; minimum progress targets: [57.7 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by  1.1259511011925538
New scales: [0.12797587 0.87202413]
   8  2.39e+00  2.53e+01  4.61e+01  1.35e+02    7.60e+01      0     100       4.33e-01     1.28e+02
geophys. misfits: 31.5 (target 30.0 [False]); 24.4 (target 30.0 [True]) | smallness misfit: 1736.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [31.5 24.4]; minimum progress targets: [54. 30.]
Updating scaling for data misfits by  1.231097889055021
New scales: [0.15302512 0.84697488]
   9  2.39e+00  2.52e+01  4.60e+01  1.35e+02    7.03e+01      0     100       2.34e+00     3.20e+02
geophys. misfits: 24.8 (target 30.0 [True]); 25.2 (target 30.0 [True]) | smallness misfit: 1552.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [24.8 25.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  1.1982239889674502
  10  2.39e+00  2.50e+01  4.67e+01  1.36e+02    6.37e+01      0     100       8.09e-01     2.59e+02
geophys. misfits: 22.6 (target 30.0 [True]); 25.4 (target 30.0 [True]) | smallness misfit: 1484.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.6 25.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  1.5002132993239634
  11  2.39e+00  2.49e+01  4.77e+01  1.39e+02    6.72e+01      0     100       6.12e-01     1.59e+02
geophys. misfits: 20.2 (target 30.0 [True]); 25.8 (target 30.0 [True]) | smallness misfit: 1407.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [20.2 25.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  1.9857818862394554
  12  2.39e+00  2.51e+01  4.90e+01  1.42e+02    6.39e+01      0     100       2.83e+00     4.56e+02
geophys. misfits: 18.7 (target 30.0 [True]); 26.2 (target 30.0 [True]) | smallness misfit: 1303.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [18.7 26.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  2.727073988357054
  13  2.39e+00  2.56e+01  5.09e+01  1.47e+02    8.25e+01      0     100       7.07e-01     3.23e+02
geophys. misfits: 18.6 (target 30.0 [True]); 26.8 (target 30.0 [True]) | smallness misfit: 1187.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [18.6 26.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  3.7259976449987384
  14  2.39e+00  2.60e+01  5.31e+01  1.53e+02    8.41e+01      0     100       2.28e+00     7.48e+02
geophys. misfits: 17.6 (target 30.0 [True]); 27.6 (target 30.0 [True]) | smallness misfit: 1078.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [17.6 27.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  5.19735754885325
  15  2.39e+00  2.70e+01  5.60e+01  1.61e+02    9.21e+01      0     100       1.01e+00     7.57e+02
geophys. misfits: 19.4 (target 30.0 [True]); 28.4 (target 30.0 [True]) | smallness misfit: 963.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.4 28.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  6.75602389950755
  16  2.39e+00  2.61e+01  5.84e+01  1.66e+02    9.37e+01      0     100       1.61e+00     1.23e+03
geophys. misfits: 18.0 (target 30.0 [True]); 27.5 (target 30.0 [True]) | smallness misfit: 821.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [18.  27.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  9.323429886829778
  17  2.39e+00  2.62e+01  6.06e+01  1.71e+02    1.10e+02      0     100       1.24e+00     1.52e+03
geophys. misfits: 18.1 (target 30.0 [True]); 27.6 (target 30.0 [True]) | smallness misfit: 637.0 (target: 200.0 [False])
Beta cooling evaluation: progress: [18.1 27.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  12.768148249815132
  18  2.39e+00  2.73e+01  6.28e+01  1.77e+02    9.68e+01      0     100       1.35e+00     2.07e+03
geophys. misfits: 19.9 (target 30.0 [True]); 28.7 (target 30.0 [True]) | smallness misfit: 475.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [19.9 28.7]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  16.296788715164183
  19  2.39e+00  2.83e+01  6.58e+01  1.85e+02    9.04e+01      0     100       8.80e-01     1.82e+03
geophys. misfits: 22.7 (target 30.0 [True]); 29.3 (target 30.0 [True]) | smallness misfit: 407.9 (target: 200.0 [False])
Beta cooling evaluation: progress: [22.7 29.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  19.099279000327073
  20  2.39e+00  2.85e+01  6.80e+01  1.91e+02    8.95e+01      0     100       6.62e-01     1.20e+03   Skip BFGS
geophys. misfits: 23.7 (target 30.0 [True]); 29.4 (target 30.0 [True]) | smallness misfit: 369.4 (target: 200.0 [False])
Beta cooling evaluation: progress: [23.7 29.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  21.851326419109387
  21  2.39e+00  2.93e+01  6.99e+01  1.96e+02    8.83e+01      0     100       9.16e-01     1.10e+03
geophys. misfits: 26.4 (target 30.0 [True]); 29.8 (target 30.0 [True]) | smallness misfit: 347.3 (target: 200.0 [False])
Beta cooling evaluation: progress: [26.4 29.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  23.39967059663779
  22  2.39e+00  2.96e+01  7.08e+01  1.99e+02    8.35e+01      0     100       7.25e-01     8.01e+02   Skip BFGS
geophys. misfits: 27.5 (target 30.0 [True]); 30.0 (target 30.0 [True]) | smallness misfit: 330.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [27.5 30. ]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  24.47286011122914
  23  2.39e+00  3.01e+01  7.14e+01  2.00e+02    7.79e+01      0     100       1.65e+00     1.32e+03
geophys. misfits: 29.7 (target 30.0 [True]); 30.1 (target 30.0 [False]) | smallness misfit: 320.1 (target: 200.0 [False])
Beta cooling evaluation: progress: [29.7 30.1]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by  1.0109803306408114
New scales: [0.15161509 0.84838491]
  24  1.19e+00  2.54e+01  7.41e+01  1.14e+02    9.10e+01      0     100       1.80e-01     2.42e+02
geophys. misfits: 14.8 (target 30.0 [True]); 27.3 (target 30.0 [True]) | smallness misfit: 352.8 (target: 200.0 [False])
Beta cooling evaluation: progress: [14.8 27.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  38.255222613993766
  25  1.19e+00  2.53e+01  8.43e+01  1.26e+02    7.74e+01      0     100       3.41e+00     8.77e+02
geophys. misfits: 14.8 (target 30.0 [True]); 27.2 (target 30.0 [True]) | smallness misfit: 302.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [14.8 27.2]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  60.00658983603364
  26  1.19e+00  2.50e+01  9.85e+01  1.42e+02    1.01e+02      0     100       7.56e-01     6.76e+02
geophys. misfits: 15.1 (target 30.0 [True]); 26.8 (target 30.0 [True]) | smallness misfit: 285.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [15.1 26.8]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  93.40296755294365
  27  1.19e+00  2.47e+01  1.18e+02  1.65e+02    1.08e+02      0     100       9.95e-01     7.86e+02
geophys. misfits: 15.8 (target 30.0 [True]); 26.3 (target 30.0 [True]) | smallness misfit: 235.7 (target: 200.0 [False])
Beta cooling evaluation: progress: [15.8 26.3]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  142.04724826902822
  28  1.19e+00  2.45e+01  1.44e+02  1.96e+02    1.16e+02      0     100       2.73e+00     2.68e+03
geophys. misfits: 13.8 (target 30.0 [True]); 26.5 (target 30.0 [True]) | smallness misfit: 231.5 (target: 200.0 [False])
Beta cooling evaluation: progress: [13.8 26.5]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  234.78973740137803
  29  1.19e+00  2.49e+01  1.85e+02  2.46e+02    1.28e+02      1     100       1.11e+00     3.18e+03
geophys. misfits: 15.7 (target 30.0 [True]); 26.6 (target 30.0 [True]) | smallness misfit: 203.6 (target: 200.0 [False])
Beta cooling evaluation: progress: [15.7 26.6]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  356.1463839412645
  30  1.19e+00  2.82e+01  2.24e+02  2.96e+02    1.32e+02      1     100       6.23e+00     1.65e+04
geophys. misfits: 36.3 (target 30.0 [False]); 26.8 (target 30.0 [True]) | smallness misfit: 179.6 (target: 200.0 [True])
Beta cooling evaluation: progress: [36.3 26.8]; minimum progress targets: [30. 30.]
Decreasing beta to counter data misfit increase.
Updating scaling for data misfits by  1.1199561910314613
New scales: [0.16676919 0.83323081]
  31  5.97e-01  2.73e+01  2.17e+02  1.57e+02    1.19e+02      1     100       2.44e-01     1.98e+03
geophys. misfits: 27.0 (target 30.0 [True]); 27.4 (target 30.0 [True]) | smallness misfit: 169.7 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [27.  27.4]; minimum progress targets: [30. 30.]
Warming alpha_pgi to favor clustering:  392.46716129644136
------------------------- STOP! -------------------------
1 : |fc-fOld| = 5.2717e+00 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 3.2632e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x|    = 1.1856e+02 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.1856e+02 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      50    <= iter          =     31
------------------------- DONE! -------------------------
/home/vsts/work/1/s/simpeg/optimization.py:1535: FutureWarning:

InexactCG.tolCG has been deprecated, please use InexactCG.cg_atol. It will be removed in version 0.26.0 of SimPEG.

/home/vsts/work/1/s/simpeg/optimization.py:1061: FutureWarning:

InexactCG.maxIterCG has been deprecated, please use InexactCG.cg_maxiter. It will be removed in version 0.26.0 of SimPEG.


Running inversion with SimPEG v0.24.1.dev37+gf5e5be547
Alpha scales: [np.float64(0.0004701811806052387), np.float64(0.0), np.float64(3.559609317027757e-06), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers:  [0.09546346 0.90453654]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales:  [0.09546346 0.90453654]
================================================= 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.98e+03  3.00e+05  0.00e+00  3.00e+05                         0           inf          inf
   1  1.98e+03  6.63e+04  2.22e+01  1.10e+05    1.41e+02      0      48       6.18e-12     5.67e-05
geophys. misfits: 90074.6 (target 30.0 [False]); 63781.2 (target 30.0 [False]) | smallness misfit: 280.2 (target: 200.0 [False])
Beta cooling evaluation: progress: [90074.6 63781.2]; minimum progress targets: [240000. 240000.]
   2  1.98e+03  9.65e+01  5.57e-01  1.20e+03    8.87e+01      0     100       1.87e-02     5.11e+02   Skip BFGS
geophys. misfits: 693.8 (target 30.0 [False]); 33.4 (target 30.0 [False]) | smallness misfit: 91.0 (target: 200.0 [True])
Beta cooling evaluation: progress: [693.8  33.4]; minimum progress targets: [72059.7 51025. ]
   3  1.98e+03  4.34e+01  1.35e-01  3.10e+02    8.08e+01      0     100       2.60e-01     5.28e+02   Skip BFGS
geophys. misfits: 139.0 (target 30.0 [False]); 33.3 (target 30.0 [False]) | smallness misfit: 68.2 (target: 200.0 [True])
Beta cooling evaluation: progress: [139.   33.3]; minimum progress targets: [555.1  30. ]
   4  1.98e+03  4.14e+01  1.30e-01  2.98e+02    8.20e+01      0     100       7.83e-03     1.27e+01
geophys. misfits: 139.9 (target 30.0 [False]); 31.1 (target 30.0 [False]) | smallness misfit: 54.2 (target: 200.0 [True])
Beta cooling evaluation: progress: [139.9  31.1]; minimum progress targets: [111.2  30. ]
Decreasing beta to counter data misfit decrase plateau.
   5  9.89e+02  3.10e+01  1.36e-01  1.66e+02    9.45e+01      0     100       1.11e-01     4.18e+01
geophys. misfits: 56.5 (target 30.0 [False]); 28.3 (target 30.0 [True]) | smallness misfit: 48.6 (target: 200.0 [True])
Beta cooling evaluation: progress: [56.5 28.3]; minimum progress targets: [111.9  30. ]
Updating scaling for data misfits by  1.0606099942810125
New scales: [0.10066703 0.89933297]
   6  9.89e+02  3.08e+01  1.36e-01  1.66e+02    3.77e+01      0     100       9.96e-01     4.40e+01   Skip BFGS
geophys. misfits: 51.9 (target 30.0 [False]); 28.4 (target 30.0 [True]) | smallness misfit: 48.8 (target: 200.0 [True])
Beta cooling evaluation: progress: [51.9 28.4]; minimum progress targets: [45.2 30. ]
Decreasing beta to counter data misfit decrase plateau.
Updating scaling for data misfits by  1.0563245728238762
New scales: [0.10573752 0.89426248]
   7  4.95e+02  2.66e+01  1.42e-01  9.69e+01    8.06e+01      0     100       2.27e-01     3.15e+01
geophys. misfits: 24.0 (target 30.0 [True]); 26.9 (target 30.0 [True]) | smallness misfit: 51.2 (target: 200.0 [True])
All targets have been reached
Beta cooling evaluation: progress: [24.  26.9]; minimum progress targets: [41.5 30. ]
Warming alpha_pgi to favor clustering:  1.1832422294450518
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.4170e+00 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 1.9760e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x|    = 8.0615e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 8.0615e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      50    <= iter          =      7
------------------------- DONE! -------------------------
/home/vsts/work/1/s/simpeg/optimization.py:1535: FutureWarning:

InexactCG.tolCG has been deprecated, please use InexactCG.cg_atol. It will be removed in version 0.26.0 of SimPEG.

/home/vsts/work/1/s/simpeg/optimization.py:1061: FutureWarning:

InexactCG.maxIterCG has been deprecated, please use InexactCG.cg_maxiter. It will be removed in version 0.26.0 of SimPEG.


Running inversion with SimPEG v0.24.1.dev37+gf5e5be547
Alpha scales: [np.float64(3.4332922983895255e-05), np.float64(0.0), np.float64(3.434957697677808e-05), np.float64(0.0)]
Calculating the scaling parameter.
Scale Multipliers:  [0.09546346 0.90453654]
/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.09546346 0.90453654]
================================================= 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.04e+06  3.00e+05  0.00e+00  3.00e+05                         0           inf          inf
   1  1.04e+06  3.78e+04  4.26e-02  8.22e+04    1.41e+02      0     100       2.93e-04     2.68e+03
geophys. misfits: 56906.9 (target 30.0 [False]); 35835.5 (target 30.0 [False])
   2  2.08e+05  4.33e+03  1.07e-01  2.65e+04    1.37e+02      0     100       6.66e-03     1.10e+02   Skip BFGS
geophys. misfits: 8525.3 (target 30.0 [False]); 3884.1 (target 30.0 [False])
   3  4.16e+04  2.85e+02  1.41e-01  6.16e+03    1.31e+02      0     100       2.15e-02     1.10e+02   Skip BFGS
geophys. misfits: 571.0 (target 30.0 [False]); 254.7 (target 30.0 [False])
   4  8.32e+03  4.04e+01  1.51e-01  1.30e+03    1.02e+02      0     100       6.83e+00     8.02e+03   Skip BFGS
geophys. misfits: 35.9 (target 30.0 [False]); 40.9 (target 30.0 [False])
   5  1.66e+03  2.42e+01  1.55e-01  2.82e+02    9.53e+01      0     100       2.57e-03     2.06e+01   Skip BFGS
geophys. misfits: 8.3 (target 30.0 [True]); 25.9 (target 30.0 [True])
All targets have been reached
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.0265e+01 <= tolF*(1+|f0|) = 3.0000e+04
0 : |xc-x_last| = 5.0672e-01 <= tolX*(1+|x0|) = 1.0000e-06
0 : |proj(x-g)-x|    = 9.5323e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 9.5323e+01 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      50    <= iter          =      5
------------------------- DONE! -------------------------
/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:302: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "b.-" (-> marker='.'). The keyword argument will take precedence.

/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:309: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "r.-" (-> marker='.'). The keyword argument will take precedence.

/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:353: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "b.-" (-> marker='.'). The keyword argument will take precedence.

/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:360: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "r.-" (-> marker='.'). The keyword argument will take precedence.

/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:367: UserWarning:

The following kwargs were not used by contour: 'label'

/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:412: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "b.-" (-> marker='.'). The keyword argument will take precedence.

/home/vsts/work/1/s/examples/10-pgi/plot_inv_1_PGI_Linear_1D_joint_WithRelationships.py:419: UserWarning:

marker is redundantly defined by the 'marker' keyword argument and the fmt string "r.-" (-> marker='.'). The keyword argument will take precedence.

import discretize as Mesh
import matplotlib.pyplot as plt
import matplotlib.lines as mlines
import numpy as np
from simpeg import (
    data_misfit,
    directives,
    inverse_problem,
    inversion,
    maps,
    optimization,
    regularization,
    simulation,
    utils,
)

# Random seed for reproductibility
np.random.seed(1)
# Mesh
N = 100
mesh = Mesh.TensorMesh([N])

# Survey design parameters
nk = 30
jk = np.linspace(1.0, 59.0, nk)
p = -0.25
q = 0.25


# Physics
def g(k):
    return np.exp(p * jk[k] * mesh.cell_centers_x) * np.cos(
        np.pi * q * jk[k] * mesh.cell_centers_x
    )


G = np.empty((nk, mesh.nC))

for i in range(nk):
    G[i, :] = g(i)

m0 = np.zeros(mesh.nC)
m0[20:41] = np.linspace(0.0, 1.0, 21)
m0[41:57] = np.linspace(-1, 0.0, 16)

poly0 = maps.PolynomialPetroClusterMap(coeffyx=np.r_[0.0, -4.0, 4.0])
poly1 = maps.PolynomialPetroClusterMap(coeffyx=np.r_[-0.0, 3.0, 6.0, 6.0])
poly0_inverse = maps.PolynomialPetroClusterMap(coeffyx=-np.r_[0.0, -4.0, 4.0])
poly1_inverse = maps.PolynomialPetroClusterMap(coeffyx=-np.r_[0.0, 3.0, 6.0, 6.0])
cluster_mapping = [maps.IdentityMap(), poly0_inverse, poly1_inverse]

m1 = np.zeros(100)
m1[20:41] = 1.0 + (poly0 * np.vstack([m0[20:41], m1[20:41]]).T)[:, 1]
m1[41:57] = -1.0 + (poly1 * np.vstack([m0[41:57], m1[41:57]]).T)[:, 1]

model2d = np.vstack([m0, m1]).T
m = utils.mkvc(model2d)

clfmapping = utils.GaussianMixtureWithNonlinearRelationships(
    mesh=mesh,
    n_components=3,
    covariance_type="full",
    tol=1e-8,
    reg_covar=1e-3,
    max_iter=1000,
    n_init=100,
    init_params="kmeans",
    random_state=None,
    warm_start=False,
    means_init=np.array(
        [
            [0, 0],
            [m0[20:41].mean(), m1[20:41].mean()],
            [m0[41:57].mean(), m1[41:57].mean()],
        ]
    ),
    verbose=0,
    verbose_interval=10,
    cluster_mapping=cluster_mapping,
)
clfmapping = clfmapping.fit(model2d)

clfnomapping = utils.WeightedGaussianMixture(
    mesh=mesh,
    n_components=3,
    covariance_type="full",
    tol=1e-8,
    reg_covar=1e-3,
    max_iter=1000,
    n_init=100,
    init_params="kmeans",
    random_state=None,
    warm_start=False,
    verbose=0,
    verbose_interval=10,
)
clfnomapping = clfnomapping.fit(model2d)

wires = maps.Wires(("m1", mesh.nC), ("m2", mesh.nC))

relatrive_error = 0.01
noise_floor = 0.0

prob1 = simulation.LinearSimulation(mesh, G=G, model_map=wires.m1)
survey1 = prob1.make_synthetic_data(
    m, relative_error=relatrive_error, noise_floor=noise_floor, add_noise=True
)

prob2 = simulation.LinearSimulation(mesh, G=G, model_map=wires.m2)
survey2 = prob2.make_synthetic_data(
    m, relative_error=relatrive_error, noise_floor=noise_floor, add_noise=True
)


dmis1 = data_misfit.L2DataMisfit(simulation=prob1, data=survey1)
dmis2 = data_misfit.L2DataMisfit(simulation=prob2, data=survey2)
dmis = dmis1 + dmis2
minit = np.zeros_like(m)

# Distance weighting
wr1 = np.sum(prob1.G**2.0, axis=0) ** 0.5 / mesh.cell_volumes
wr1 = wr1 / np.max(wr1)
wr2 = np.sum(prob2.G**2.0, axis=0) ** 0.5 / mesh.cell_volumes
wr2 = wr2 / np.max(wr2)

reg_simple = regularization.PGI(
    mesh=mesh,
    gmmref=clfmapping,
    gmm=clfmapping,
    approx_gradient=True,
    wiresmap=wires,
    non_linear_relationships=True,
    weights_list=[wr1, wr2],
)

opt = optimization.ProjectedGNCG(
    maxIter=50,
    tolX=1e-6,
    maxIterCG=100,
    tolCG=1e-3,
    lower=-10,
    upper=10,
)

invProb = inverse_problem.BaseInvProblem(dmis, reg_simple, opt)

# directives
scales = directives.ScalingMultipleDataMisfits_ByEig(
    chi0_ratio=np.r_[1.0, 1.0], verbose=True, n_pw_iter=10
)
scaling_schedule = directives.JointScalingSchedule(verbose=True)
alpha0_ratio = np.r_[1e6, 1e4, 1, 1]
alphas = directives.AlphasSmoothEstimate_ByEig(
    alpha0_ratio=alpha0_ratio, n_pw_iter=10, verbose=True
)
beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-5, n_pw_iter=10)
betaIt = directives.PGI_BetaAlphaSchedule(
    verbose=True,
    coolingFactor=2.0,
    progress=0.2,
)
targets = directives.MultiTargetMisfits(verbose=True)
petrodir = directives.PGI_UpdateParameters(update_gmm=False)

# Setup Inversion
inv = inversion.BaseInversion(
    invProb,
    directiveList=[alphas, scales, beta, petrodir, targets, betaIt, scaling_schedule],
)

mcluster_map = inv.run(minit)

# Inversion with no nonlinear mapping
reg_simple_no_map = regularization.PGI(
    mesh=mesh,
    gmmref=clfnomapping,
    gmm=clfnomapping,
    approx_gradient=True,
    wiresmap=wires,
    non_linear_relationships=False,
    weights_list=[wr1, wr2],
)

opt = optimization.ProjectedGNCG(
    maxIter=50,
    tolX=1e-6,
    maxIterCG=100,
    tolCG=1e-3,
    lower=-10,
    upper=10,
)

invProb = inverse_problem.BaseInvProblem(dmis, reg_simple_no_map, opt)

# directives
scales = directives.ScalingMultipleDataMisfits_ByEig(
    chi0_ratio=np.r_[1.0, 1.0], verbose=True, n_pw_iter=10
)
scaling_schedule = directives.JointScalingSchedule(verbose=True)
alpha0_ratio = np.r_[100.0 * np.ones(2), 1, 1]
alphas = directives.AlphasSmoothEstimate_ByEig(
    alpha0_ratio=alpha0_ratio, n_pw_iter=10, verbose=True
)
beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-5, n_pw_iter=10)
betaIt = directives.PGI_BetaAlphaSchedule(
    verbose=True,
    coolingFactor=2.0,
    progress=0.2,
)
targets = directives.MultiTargetMisfits(
    chiSmall=1.0, TriggerSmall=True, TriggerTheta=False, verbose=True
)
petrodir = directives.PGI_UpdateParameters(update_gmm=False)

# Setup Inversion
inv = inversion.BaseInversion(
    invProb,
    directiveList=[alphas, scales, beta, petrodir, targets, betaIt, scaling_schedule],
)

mcluster_no_map = inv.run(minit)

# WeightedLeastSquares Inversion

reg1 = regularization.WeightedLeastSquares(
    mesh, alpha_s=1.0, alpha_x=1.0, mapping=wires.m1, weights={"cell_weights": wr1}
)

reg2 = regularization.WeightedLeastSquares(
    mesh, alpha_s=1.0, alpha_x=1.0, mapping=wires.m2, weights={"cell_weights": wr2}
)

reg = reg1 + reg2

opt = optimization.ProjectedGNCG(
    maxIter=50,
    tolX=1e-6,
    maxIterCG=100,
    tolCG=1e-3,
    lower=-10,
    upper=10,
)

invProb = inverse_problem.BaseInvProblem(dmis, reg, opt)

# directives
alpha0_ratio = np.r_[1, 1, 1, 1]
alphas = directives.AlphasSmoothEstimate_ByEig(
    alpha0_ratio=alpha0_ratio, n_pw_iter=10, verbose=True
)
scales = directives.ScalingMultipleDataMisfits_ByEig(
    chi0_ratio=np.r_[1.0, 1.0], verbose=True, n_pw_iter=10
)
scaling_schedule = directives.JointScalingSchedule(verbose=True)
beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-5, n_pw_iter=10)
beta_schedule = directives.BetaSchedule(coolingFactor=5.0, coolingRate=1)
targets = directives.MultiTargetMisfits(
    TriggerSmall=False,
    verbose=True,
)

# Setup Inversion
inv = inversion.BaseInversion(
    invProb,
    directiveList=[alphas, scales, beta, targets, beta_schedule, scaling_schedule],
)

mtik = inv.run(minit)


# Final Plot
fig, axes = plt.subplots(3, 4, figsize=(25, 15))
axes = axes.reshape(12)
left, width = 0.25, 0.5
bottom, height = 0.25, 0.5
right = left + width
top = bottom + height

axes[0].set_axis_off()
axes[0].text(
    0.5 * (left + right),
    0.5 * (bottom + top),
    ("Using true nonlinear\npetrophysical relationships"),
    horizontalalignment="center",
    verticalalignment="center",
    fontsize=20,
    color="black",
    transform=axes[0].transAxes,
)

axes[1].plot(mesh.cell_centers_x, wires.m1 * mcluster_map, "b.-", ms=5, marker="v")
axes[1].plot(mesh.cell_centers_x, wires.m1 * m, "k--")
axes[1].set_title("Problem 1")
axes[1].legend(["Recovered Model", "True Model"], loc=1)
axes[1].set_xlabel("X")
axes[1].set_ylabel("Property 1")

axes[2].plot(mesh.cell_centers_x, wires.m2 * mcluster_map, "r.-", ms=5, marker="v")
axes[2].plot(mesh.cell_centers_x, wires.m2 * m, "k--")
axes[2].set_title("Problem 2")
axes[2].legend(["Recovered Model", "True Model"], loc=1)
axes[2].set_xlabel("X")
axes[2].set_ylabel("Property 2")

x, y = np.mgrid[-1:1:0.01, -4:2:0.01]
pos = np.empty(x.shape + (2,))
pos[:, :, 0] = x
pos[:, :, 1] = y
CS = axes[3].contour(
    x,
    y,
    np.exp(clfmapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
    100,
    alpha=0.25,
    cmap="viridis",
)
cs_proxy = mlines.Line2D([], [], label="True Petrophysical Distribution")

ps = axes[3].scatter(
    wires.m1 * mcluster_map,
    wires.m2 * mcluster_map,
    marker="v",
    label="Recovered model crossplot",
)
axes[3].set_title("Petrophysical Distribution")
axes[3].legend(handles=[cs_proxy, ps])
axes[3].set_xlabel("Property 1")
axes[3].set_ylabel("Property 2")

axes[4].set_axis_off()
axes[4].text(
    0.5 * (left + right),
    0.5 * (bottom + top),
    ("Using a pure\nGaussian distribution"),
    horizontalalignment="center",
    verticalalignment="center",
    fontsize=20,
    color="black",
    transform=axes[4].transAxes,
)

axes[5].plot(mesh.cell_centers_x, wires.m1 * mcluster_no_map, "b.-", ms=5, marker="v")
axes[5].plot(mesh.cell_centers_x, wires.m1 * m, "k--")
axes[5].set_title("Problem 1")
axes[5].legend(["Recovered Model", "True Model"], loc=1)
axes[5].set_xlabel("X")
axes[5].set_ylabel("Property 1")

axes[6].plot(mesh.cell_centers_x, wires.m2 * mcluster_no_map, "r.-", ms=5, marker="v")
axes[6].plot(mesh.cell_centers_x, wires.m2 * m, "k--")
axes[6].set_title("Problem 2")
axes[6].legend(["Recovered Model", "True Model"], loc=1)
axes[6].set_xlabel("X")
axes[6].set_ylabel("Property 2")

CSF = axes[7].contour(
    x,
    y,
    np.exp(clfmapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
    100,
    alpha=0.5,
    label="True Petro. Distribution",
)
CS = axes[7].contour(
    x,
    y,
    np.exp(clfnomapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
    500,
    cmap="viridis",
    linestyles="--",
)
axes[7].scatter(
    wires.m1 * mcluster_no_map,
    wires.m2 * mcluster_no_map,
    marker="v",
    label="Recovered model crossplot",
)
cs_modeled_proxy = mlines.Line2D(
    [], [], linestyle="--", label="Modeled Petro. Distribution"
)

axes[7].set_title("Petrophysical Distribution")
axes[7].legend(handles=[cs_proxy, cs_modeled_proxy, ps])
axes[7].set_xlabel("Property 1")
axes[7].set_ylabel("Property 2")

# Tikonov

axes[8].set_axis_off()
axes[8].text(
    0.5 * (left + right),
    0.5 * (bottom + top),
    ("Least-Squares\n~Using a single cluster"),
    horizontalalignment="center",
    verticalalignment="center",
    fontsize=20,
    color="black",
    transform=axes[8].transAxes,
)

axes[9].plot(mesh.cell_centers_x, wires.m1 * mtik, "b.-", ms=5, marker="v")
axes[9].plot(mesh.cell_centers_x, wires.m1 * m, "k--")
axes[9].set_title("Problem 1")
axes[9].legend(["Recovered Model", "True Model"], loc=1)
axes[9].set_xlabel("X")
axes[9].set_ylabel("Property 1")

axes[10].plot(mesh.cell_centers_x, wires.m2 * mtik, "r.-", ms=5, marker="v")
axes[10].plot(mesh.cell_centers_x, wires.m2 * m, "k--")
axes[10].set_title("Problem 2")
axes[10].legend(["Recovered Model", "True Model"], loc=1)
axes[10].set_xlabel("X")
axes[10].set_ylabel("Property 2")

CS = axes[11].contour(
    x,
    y,
    np.exp(clfmapping.score_samples(pos.reshape(-1, 2)).reshape(x.shape)),
    100,
    alpha=0.25,
    cmap="viridis",
)
axes[11].scatter(wires.m1 * mtik, wires.m2 * mtik, marker="v")
axes[11].set_title("Petro Distribution")
axes[11].legend(handles=[cs_proxy, ps])
axes[11].set_xlabel("Property 1")
axes[11].set_ylabel("Property 2")
plt.subplots_adjust(wspace=0.3, hspace=0.3, top=0.85)
plt.show()

Total running time of the script: (0 minutes 38.440 seconds)

Estimated memory usage: 298 MB

Gallery generated by Sphinx-Gallery