Output models #15
Comments
|
Hi @79seismo, All models and misfits are output by default (attributes Without any prior information, it would be good to set rather unconstraining velocity bounds (yet still realistic, i.e., we do not expect 5 km/s at the surface) and constant thicknesses for each layer. Also, set a large population size for the optimizer (e.g., 10 times the number of layers). Vp is calculated from Vs and Poisson's ratio. By default, density is calculated using Nafe-Drake's equation, but you can set your own model ( Not sure what you mean by "benchmark", but dispersion curves calculated by |
|
Thanks very much Keurfon. This is useful information. By benchmark I meant comparing results with an established code, which is what you may have done with CPS. I tried inverting for the model with 10 layers in your example in disba, and although the misfit gets minimised well, the model fit isn't that great for CPSO and is a bit better with NA from what I can tell. Any tips on how I can improve the fit? I set the search ranges as follows for all the layers Can I get a copy of your PhD thesis as well please? Thanks! |
|
Please note that thickness is defined in km and velocity in km/s. For the sample problem, the model layers can be configured as follows (assuming 10 layers of constant thickness h=5m down to 50 m depth): for _ in range(10):
model.add(Layer([0.005, 0.005], [0.1, 3.0]))Set the population size to 100 and increase the number of iterations to 1000-2000 to be sure to reach convergence. You may also try different optimizers as you did with NA. CMA-ES seems to work quite well here: model.configure(
optimizer="cmaes",
misfit="rmse",
density=lambda vp: 2.0,
optimizer_args={
"popsize": 100,
"maxiter": 2000,
"workers": -1,
"seed": 0,
},
)You can also increase the number of runs (e.g., to 3) to start from different initial populations. res = model.invert(curves, maxrun=3, split_results=True)
res = min(res, key=lambda x: x.misfit) # Get the results for the run with the lowest misfitNote that you likely will never get the exact true model due to the non-uniqueness of the solution. Note that CMA-ES may not always work well. CMA-ES can be compared to a gradient based optimizer as it tries to generate more samples along the direction of the steepest descent. |
|
Thanks! much appreciated. Encountered the following for a VDCMA run, possible bug? I was using population size=100 and iternations = 1000 Also, is there a way to isolate an ensemble of models within a certain threshold of the minimum misfit, say within 20% of the minimum misfit? What does this do? Also, rather than showing all models with 'show=all' as given below, is there a way to isolate and show the models that were extracted within the earlier-mentioned threshold; i.e. say within 20% of the minimum misfit? Thank you! |
|
I guess VDCMA does not like when all models have infinite misfit (happens when constraints are not satisfied). I will need to look into that. The method value = 1.2 * res.misfitNote that: res = res.threshold(value)overwrites the result variable Then, invoking the method res.threshold(value).plot_models("vs", zmax=zmax, show="all", ax=ax[0], plot_args={"cmap": cmap) |
|
Neat! Thanks very much for your help! |
|
I still haven't been able to get a good inversion result for a synthetic model I'm trying and it looks like the following - a simple model with a velocity gradient. I tried all algorithms with population size = 100 to 150 and iterations 500 to 2500. Target computational costs don't let me go above these parameters. Any suggestions to improve the inversion? Also, the plotting routines are quite computationally intensive. I wonder if it can be made more efficient. Thanks!
Fundamental mode Rayleigh Group dispersion curve from disba for this model: Period: Velocity: |
|
I inverted your group velocities with 20 layers: import numpy as np
from evodcinv import EarthModel, Layer, Curve
import matplotlib.pyplot as plt
period = np.array([ 0.1 , 0.10476158, 0.10974988, 0.1149757 , 0.12045035, 0.12618569, 0.13219411, 0.13848864, 0.14508288, 0.15199111, 0.15922828, 0.16681005, 0.17475284, 0.18307383, 0.19179103, 0.2009233 , 0.21049041, 0.22051307, 0.23101297, 0.24201283, 0.25353645, 0.26560878, 0.27825594, 0.29150531, 0.30538555, 0.31992671, 0.33516027, 0.35111917, 0.36783798, 0.38535286, 0.40370173, 0.42292429, 0.44306215, 0.46415888, 0.48626016, 0.5094138 , 0.53366992, 0.55908102, 0.58570208, 0.61359073, 0.64280731, 0.67341507, 0.70548023, 0.7390722 , 0.77426368, 0.81113083, 0.84975344, 0.89021509, 0.93260335, 0.97700996, 1.02353102, 1.07226722, 1.12332403, 1.17681195, 1.23284674, 1.29154967, 1.35304777, 1.41747416, 1.48496826, 1.55567614, 1.62975083, 1.70735265, 1.78864953, 1.87381742, 1.96304065, 2.05651231, 2.15443469, 2.25701972, 2.36448941, 2.47707636, 2.59502421, 2.71858824, 2.84803587, 2.98364724, 3.12571585, 3.27454916, 3.43046929, 3.59381366, 3.76493581, 3.94420606, 4.1320124 , 4.32876128, 4.53487851, 4.75081016, 4.97702356, 5.21400829, 5.46227722, 5.72236766, 5.9948425 , 6.28029144, 6.57933225, 6.8926121 , 7.22080902, 7.56463328, 7.92482898, 8.30217568, 8.69749003, 9.11162756, 9.54548457, 10. ])
group = np.array([1.09786535, 1.20770009, 1.30751688, 1.39670049, 1.47560087, 1.54520622, 1.60664691, 1.66103997, 1.70934501, 1.7524971 , 1.79118683, 1.82605584, 1.85757589, 1.88622391, 1.91242616, 1.93634936, 1.95832355, 1.97861568, 1.99732756, 2.01468714, 2.03079722, 2.04576357, 2.05973416, 2.07277005, 2.08493231, 2.09632235, 2.10704698, 2.11707513, 2.12650911, 2.13536061, 2.14373248, 2.15163577, 2.15908177, 2.1660324 , 2.17268075, 2.1789452 , 2.18487845, 2.19044092, 2.19573366, 2.20076534, 2.20553736, 2.21001281, 2.21428938, 2.21832549, 2.2221736 , 2.22583808, 2.22927476, 2.23253889, 2.23567875, 2.23865024, 2.24150667, 2.24420385, 2.24669491, 2.24912634, 2.25145642, 2.2536359 , 2.25576521, 2.25775087, 2.25964183, 2.26144286, 2.26315383, 2.26477708, 2.26631494, 2.26776975, 2.26923784, 2.27057825, 2.27183783, 2.27302139, 2.2742255 , 2.27530642, 2.27636073, 2.27734372, 2.2783025 , 2.27923952, 2.28015721, 2.28090897, 2.28173825, 2.28250089, 2.28324658, 2.28387826, 2.28454512, 2.285195 , 2.28573319, 2.28635639, 2.28686785, 2.28741698, 2.28785676, 2.2883816 , 2.28879949, 2.28920765, 2.2896535 , 2.2900422 , 2.29042116, 2.29074295, 2.29110732, 2.2914145 , 2.29171193, 2.2920045 , 2.2922922 , 2.29257016])
model = EarthModel()
for _ in range(20):
model.add(Layer(0.0015, [0.1, 3.0]))
model.configure(
optimizer="cpso",
misfit="rmse",
density=lambda vp: 2.0,
optimizer_args={
"popsize": 20,
"maxiter": 200,
"workers": -1,
"seed": 0,
},
)
curves = [Curve(period, group, 0, "rayleigh", "group")]
res = model.invert(curves, maxrun=3)
res.plot_curve(period, 0, "rayleigh", "group", show="best")
plt.semilogx(period, group, ls="--")The results look good to me: |
|
Yes, the dispersion curve gets recovered well but not the velocity profile? I was doing this to kind of validate the inversions. |
|
The solutions of inversion problems are, in general, non-unique (there is an infinite number of solutions that can match your data). We can't expect the code to find the exact solution, even for synthetic cases. To get "better" results, we need to provide prior information (via better boundary constraints for instance, or additional regularization terms). Without any prior information, the code, as "dumb" as it is, returns one model that fits your data. |
|
Indeed, and is there a way to modify the misfit function? Alternatively, it's easier to quantify uncertainty if I can retrieve all the models within a certain threshold of the misfit. Is it possible to retrieve models after res.threshold()? |
|
You can look at the option "extra_terms" in "configure", this option allows to add additional misfit terms (e.g., to penalize models).
|
|
Brilliant, thank you! |
|
I've updated the code and improved the option model = EarthModel()
for _ in range(20):
model.add(Layer(0.0015, [0.1, 3.0]))
model.configure(
optimizer="de",
misfit="rmse",
density=lambda vp: 2.0,
optimizer_args={
"popsize": 50,
"maxiter": 500,
"workers": -1,
"seed": 0,
"strategy": "rand1bin",
},
increasing_velocity=True,
)
curves = [Curve(period, group, 0, "rayleigh", "group")]
res = model.invert(curves, maxrun=3) |
|
Thanks. I guess I have to reinstall the codes if you've updated them? increasing_velocity=True ensures that the inversion produces a best-fitting positive velocity gradient? So to ensure that I am not omitting structure with low velocity layers, I have to run a separate inversion with increasing_velocity=False and compare the misfits as a strategy? Yes, I read your paper and that's what got me interested in your codes. |
|
Yes, please run the command |
Here an example where we assume we have a linear prior model (Vs(z=0m)=0.1 km/s, Vs(z=40m)=3 km/s). with a regularization factor of 0.001 (can be tuned): model = EarthModel()
for _ in range(10):
model.add(Layer(0.003, [0.1, 3.0]))
def prior(x):
vel = model.transform(x)
z = np.insert(vel[:-1, 0].cumsum(), 0, 0.0)
vprior = np.interp(z, [0.0, 0.04], [0.1, 3.0])
return 1.0e-3 * np.square(vel[:, 2] - vprior).sum()
model.configure(
optimizer="cpso",
misfit="norm2",
density=lambda vp: 2.0,
optimizer_args={
"popsize": 40,
"maxiter": 500,
"workers": -1,
"seed": 0,
},
extra_terms=[prior],
)
curves = [Curve(period, group, 0, "rayleigh", "group")]
res = model.invert(curves, maxrun=3)
res.plot_curve(period, 0, "rayleigh", "group", show="best")
plt.semilogx(period, group, ls="--")Note that I am here using a |
|
I am not sure if I am following what's happening within |
|
with I just updated from evodcinv import factory
model = EarthModel()
for _ in range(20):
model.add(Layer(0.0015, [0.1, 3.0]))
model.configure(
optimizer="de",
misfit="norm2",
density=lambda vp: 2.0,
optimizer_args={
"popsize": 50,
"maxiter": 500,
"workers": -1,
"seed": 0,
},
increasing_velocity=True,
extra_terms=[
lambda x: factory.smooth(x, alpha=1.0e-3),
lambda x: factory.prior(x, [0.0, 0.04], [0.1, 3.0], alpha=1.0e-3)
],
)
curves = [Curve(period, group, 0, "rayleigh", "group")]
res = model.invert(curves, maxrun=3) |
|
That's brilliant! Thank you! We should write a paper together to demonstrate the Utility of evodcinv in different applications, not just geophysics! |
|
I think there should be an option to run "maxrun=" with different seeds. Does it make sense to run multiple inversions with the same seed? Shouldn't the starting point be different for different runs? Thanks! |
|
The seed is set before the "inversion" loop, so each inversion is using a different set of random numbers. It's easy to check that since each run is returning a different optimal solution. |
|
Thanks, can you please explain the functionality of "seed":0 in model.configure(). Is it for the first run? |
|
The purpose of this option is reproducibility, so that every time you run the same script, you get the same output. It's calling |
|
Hi, there are instances when an inversion run terminates before reaching 100% (i.e. population size x iterations) if the misfit approaches a small value, e.g. 10^-9, early on. In those instances, not all the model arrays in res (i.e. res.models[i]) are populated. Instead, beyond the model that approached the minimum misfit, all model arrays in res will be populated with zeros. It would be good to trim the res.models and everything else in res by deleting unused array locations that are filled with zeros, for example with np.where method (actually the nonzero value is len(res.xs)), so that res.threshold can be used with statistical measures without having to include intermediate code by hand. Thanks! |
|
I think the default minimum misfit value in |
|
Here's an example dispersion curve I've tried producing 10^-9 misfit. As you can see, this spans lower than usual periods. I just realized that the inversion is fitting only the first couple of points of the dispersion curve and hence the very low misfit. Why would it do that? length of period array: 100 length of group velocity array: 100 |
|
Hi @79seismo, Please update |
|
Thanks, that appeared to have fixed the issue. Also, I realized that dispersion curve predictions are not stored during the inversion and are recalculated within the plotting routines. This is a bit inefficient as evodcinv is then computing the same dispersion curves twice. Why not store dispersion curves within the inversion itself and use them in the plotting routines? |
|
The first reason is API limitation: most optimization libraries (e.g., |
Thanks for that explanation. |
|
Are there reasons why a run would quit half way through other than approaching a machine precision minimum as discussed before? |
|
It depends on the algorithm I guess. For instance, for Differential Evolution, it may be caused by a lack of diversity in the population at a given iteration, i.e., all the models are too similar so recombining models would be useless. |
|
@keurfonluu is it possible to assemble the minimum misfit models across from multiple runs instead of assembling the lowest misfit models from the run that has the minimum misfit? The former should give you more precise models than the latter especially if the minimum misfit is not that different across multiple runs. |
|
Hi @79seismo, I am not sure to correctly understand your question. By default, when |
|
Thanks, exactly what I was looking for. |
|
Hi,
|
|
Hi @79seismo, Hard to tell the reason without the code used to generate this figure and the outputs... |
|
Sorry, here's where I define scalebar and the |
|
When there are inf values in |
|
Is there a easy way of plotting the y axis of velocity models in meters rather than km? |
|
I have a situation where the inversion result depends on the seed value all other things held constant, i.e., I change nothing in the model and inversion controls other than the seed value. Is this expected behavior? I was of the impression that ,regardless of the seed value, the solution should converge? Increasing the number of iterations and population size doesn't seem to matter. e.g. The model has 30 layers, Poisson's ratio is part of inversion parameters, and it is a constrained inversion using the factory extra_terms for progressive velocity increase. result with seed = 0 result with seed = 10000 |
|
It's a misconception to think that evolutionary algorithms, as global optimizers, should always converge toward the global optimum. Although global optimizers are less sensitive to initial solution compared to local optimizers, they are still sensitive to it. Evolutionary algorithms can be trapped in a local minimum, and the lack of diversity in the population can make it unable to escape from it, which is likely what happens in your first case. The good practice is to fix the seed (for reproducibility) and set the number of runs to a value greater than 1 (as much as your time frame allows you). The larger, the better as you increase the probability to start with a good initial population. |
|
Thanks. Yes, this was after multiple runs. |
Hi Keurfon,
How do you output all models + misfit by iteration? I'm looking to create a figure like on the homepage with all models color-coded by misfit and Fig 7 in your 2020 paper in Prospects. Also, what's a good starting model in the event that there is no information? I am looking at using your code for near-surface investigations; < ~100 m depth
Also, how do you estimate Vp and Density of layers in the inversion? Are you using a scaling relationship between these parameters and Vs?
Has your code, both disba and evocdinv been benchmarked?
Thanks!
Januka
The text was updated successfully, but these errors were encountered: