Repository: keurfonluu/evodcinv (language: python)

Issue:
Output models 
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



Discussion / hints:
Hi @79seismo,

All models and misfits are output by default (attributes `models` and `misfits` of `evodcinv.InversionResult`).
Please check the script `.github/sample.py` to see how the figures in the README have been generated.

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 (`density` option in method `configure`).

Not sure what you mean by "benchmark", but dispersion curves calculated by `disba` are consistent with the ones from CPS for the velocity models I have tested. `stochopy` (the optimization library used here) has been thoroughly tested during my PhD. Since `evodcinv` is mainly relying on these two libraries, the forward modeling and the optimization algorithms are assumed to be working, thus the inversion code here should work properly. If it means anything, the example in the README is actually a sample problem from Geopsy's Wiki.
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 `model.add(Layer([1, 20], [1, 5.0]))`

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):

```python
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:

```python
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.

```python
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 misfit
```

Note that you likely will never get the exact true model due to the non-uniqueness of the solution.

```
--------------------------------------------------------------------------------
Best model out of 200000 models (1 run)

Velocity model                                    Model parameters
----------------------------------------          ------------------------------
         d        vp        vs       rho                   d        vs        nu
      [km]    [km/s]    [km/s]   [g/cm3]                [km]    [km/s]       [-]
----------------------------------------          ------------------------------
    0.0050    0.3372    0.2065    2.0000              0.0050    0.2065    0.2000
    0.0050    0.3506    0.2147    2.0000              0.0050    0.2147    0.2000
    0.0050    0.5330    0.2263    2.0000              0.0050    0.2263    0.3901
    0.0050    0.5984    0.2443    2.0000              0.0050    0.2443    0.4000
    0.0050    0.2664    0.1632    2.0000              0.0050    0.1632    0.2000
    0.0050    1.3206    0.5391    2.0000              0.0050    0.5391    0.4000
    0.0050    7.3476    3.0000    2.0000              0.0050    3.0000    0.4000
    0.0050    5.8546    2.3907    2.0000              0.0050    2.3907    0.3999
    0.0050    0.3920    0.1600    2.0000              0.0050    0.1600    0.4000
    1.0000    2.0985    0.8567    2.0000                   -    0.8567    0.4000
----------------------------------------          ------------------------------

Number of layers: 10
Number of parameters: 29
Best model misfit: 0.0006
--------------------------------------------------------------------------------
```

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?
<img width="1101" alt="Screen Shot 2023-08-03 at 1 20 25 pm" src="https://github.com/keurfonluu/evodcinv/assets/31010305/6eae4b5f-0ba1-4fd4-81bb-7eeb8a24a639">

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?
`res = res.threshold(0.1)`

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?
`res.plot_model(
    "vs",
    zmax=zmax,
    show="all",
    ax=ax[0],
    plot_args={"cmap": cmap},`

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 ``threshold`` does exactly what you want: it extracts models that have misfits lower than the input value. If you want the threshold to be 20% of the minimum misfit:

```python
value = 1.2 * res.misfit
```

Note that:

```python
res = res.threshold(value)
```

overwrites the result variable ``res`` with the result filtered by the method ``threshold``. You may want to save the filtered results in another variable, e.g.:

```
res2 = res.threshold(value)
```

Then, invoking the method ``plot_model`` from the filtered result will only show the models that satisfied your threshold criterion. Alternatively, you can pipe the methods as follows:

```python
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!

`velocity_model = np.array([
    [ 0.001 , 0.346 , 0.2 , 1.4099 ],
[ 0.001 , 0.4786 , 0.2767 , 1.4116 ],
[ 0.001 , 0.6113 , 0.3533 , 1.4132 ],
[ 0.001 , 0.7439 , 0.43 , 1.4145 ],
[ 0.001 , 0.8765 , 0.5067 , 1.4158 ],
[ 0.001 , 1.0092 , 0.5833 , 1.4169 ],
[ 0.001 , 1.1418 , 0.66 , 1.418 ],
[ 0.001 , 1.2744 , 0.7367 , 1.419 ],
[ 0.001 , 1.4071 , 0.8133 , 1.42 ],
[ 0.001 , 1.5397 , 0.89 , 1.4209 ],
[ 0.001 , 1.6723 , 0.9667 , 1.4218 ],
[ 0.001 , 1.805 , 1.0433 , 1.4226 ],
[ 0.001 , 1.9376 , 1.12 , 1.4234 ],
[ 0.001 , 2.0702 , 1.1967 , 1.4242 ],
[ 0.001 , 2.2029 , 1.2733 , 1.425 ],
[ 0.001 , 2.3355 , 1.35 , 1.4257 ],
[ 0.001 , 2.4681 , 1.4267 , 1.4264 ],
[ 0.001 , 2.6008 , 1.5033 , 1.4271 ],
[ 0.001 , 2.7334 , 1.58 , 1.4278 ],
[ 0.001 , 2.866 , 1.6567 , 1.4285 ],
[ 0.001 , 2.9987 , 1.7333 , 1.4291 ],
[ 0.001 , 3.1313 , 1.81 , 1.4298 ],
[ 0.001 , 3.2639 , 1.8867 , 1.4304 ],
[ 0.001 , 3.3966 , 1.9633 , 1.431 ],
[ 0.001 , 3.5292 , 2.04 , 1.4316 ],
[ 0.001 , 3.6618 , 2.1167 , 1.4322 ],
[ 0.001 , 3.7945 , 2.1933 , 1.4328 ],
[ 0.001 , 3.9271 , 2.27 , 1.4334 ],
[ 0.001 , 4.0597 , 2.3467 , 1.4339 ],
[ 0.001 , 4.1924 , 2.4233 , 1.4345 ],
[ 0.001 , 4.325 , 2.5 , 1.435 ],
]) `

Fundamental mode Rayleigh Group dispersion curve from disba for this model:

Period:
`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.        ])`

Velocity:
`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])`
I inverted your group velocities with 20 layers:

```python
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:
![output](https://github.com/keurfonluu/evodcinv/assets/16566206/4b744c82-3848-4938-b7b6-172384ddfb4c)


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).

`res.threshold` simply filters out all the models that do not satisfy your threshold criterion. All the acceptable models are then accessible in `res.models`.


Brilliant, thank you!
I've updated the code and improved the option ``increasing_velocity`` (as usual, ``pip install evodcinv -U``). The models are now initialized following [this paper](https://onlinelibrary.wiley.com/doi/10.1111/1365-2478.12866). Interestingly, this options seems to work quite well with Differential Evolution for the three sample problems I have tested (again, that may not be true for all problems).

```python
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)
```

```
--------------------------------------------------------------------------------
Best model out of 75000 models (3 runs)

Velocity model                                    Model parameters
----------------------------------------          ------------------------------
         d        vp        vs       rho                   d        vs        nu
      [km]    [km/s]    [km/s]   [g/cm3]                [km]    [km/s]       [-]
----------------------------------------          ------------------------------
    0.0015    0.4773    0.2389    2.0000              0.0015    0.2389    0.3329
    0.0015    0.6175    0.3148    2.0000              0.0015    0.3148    0.3243
    0.0015    0.7540    0.4099    2.0000              0.0015    0.4099    0.2903
    0.0015    0.9715    0.5724    2.0000              0.0015    0.5724    0.2341
    0.0015    1.2553    0.7678    2.0000              0.0015    0.7678    0.2011
    0.0015    1.4520    0.7960    2.0000              0.0015    0.7960    0.2851
    0.0015    1.6085    0.8779    2.0000              0.0015    0.8779    0.2879
    0.0015    1.8379    1.0229    2.0000              0.0015    1.0229    0.2756
    0.0015    2.2715    1.0593    2.0000              0.0015    1.0593    0.3611
    0.0015    2.0950    1.1280    2.0000              0.0015    1.1280    0.2959
    0.0015    2.5391    1.1726    2.0000              0.0015    1.1726    0.3644
    0.0015    2.4278    1.3305    2.0000              0.0015    1.3305    0.2854
    0.0015    2.6118    1.4435    2.0000              0.0015    1.4435    0.2801
    0.0015    2.6806    1.5059    2.0000              0.0015    1.5059    0.2694
    0.0015    3.2776    1.5632    2.0000              0.0015    1.5632    0.3528
    0.0015    3.0416    1.6856    2.0000              0.0015    1.6856    0.2784
    0.0015    3.3495    1.8566    2.0000              0.0015    1.8566    0.2783
    0.0015    3.7398    2.1915    2.0000              0.0015    2.1915    0.2385
    0.0015    3.9894    2.3113    2.0000              0.0015    2.3113    0.2474
    1.0000    4.4713    2.4880    2.0000                   -    2.4880    0.2758
----------------------------------------          ------------------------------

Number of layers: 20
Number of parameters: 59
Best model misfit: 0.0002
--------------------------------------------------------------------------------
```
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 ``pip install evodcinv -U`` to update the package.
I guess that can be a strategy as well. I would also recommend not to invert for that many layers due to the non-unicity of the inversion problem. In general, simpler models are preferable (e.g., 10 would be enough).
> You can look at the option "extra_terms" in "configure", this option allows to add additional misfit terms (e.g., to penalize models).

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):

```python
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 ``"norm2"`` misfit function.

```
--------------------------------------------------------------------------------
Best model out of 60000 models (3 runs)

Velocity model                                    Model parameters
----------------------------------------          ------------------------------
         d        vp        vs       rho                   d        vs        nu
      [km]    [km/s]    [km/s]   [g/cm3]                [km]    [km/s]       [-]
----------------------------------------          ------------------------------
    0.0030    0.4892    0.2925    2.0000              0.0030    0.2925    0.2217
    0.0030    0.9116    0.5113    2.0000              0.0030    0.5113    0.2705
    0.0030    1.2563    0.7523    2.0000              0.0030    0.7523    0.2205
    0.0030    1.7016    0.9394    2.0000              0.0030    0.9394    0.2808
    0.0030    2.1432    1.0716    2.0000              0.0030    1.0716    0.3333
    0.0030    2.3475    1.2724    2.0000              0.0030    1.2724    0.2920
    0.0030    2.9727    1.5176    2.0000              0.0030    1.5176    0.3238
    0.0030    3.1543    1.6221    2.0000              0.0030    1.6221    0.3203
    0.0030    3.5203    1.9024    2.0000              0.0030    1.9024    0.2938
    1.0000    4.4075    2.4936    2.0000                   -    2.4936    0.2646
----------------------------------------          ------------------------------

Number of layers: 10
Number of parameters: 29
Best model misfit: 0.0004
--------------------------------------------------------------------------------
```
I am not sure if I am following what's happening within `def prior(x)`. If you can add some comments, that'd be great. Ultimately, I'm guessing you are creating a misfit term that is `norm2 + extra_terms`, correct? In this case, rather than fixing the velocity profile to a prior, I think it's useful to add a penalty term such that the velocity jumps across layers are restricted. For example see the second term in eq(2) here: https://www.sciencedirect.com/science/article/pii/S0040195117302482. The philosophy behind this is that drastic velocity changes don't occur across layers particularly if you have a model parameterisation with fine layers.  

Fix patch:
diff --git a/.github/sample.png b/.github/sample.png
index d4c546a..0c5ae87 100644
Binary files a/.github/sample.png and b/.github/sample.png differ
diff --git a/.github/sample.py b/.github/sample.py
index 14c5e04..372fc6e 100644
--- a/.github/sample.py
+++ b/.github/sample.py
@@ -41,7 +41,7 @@
 model.configure(
     optimizer="cpso",
     misfit="rmse",
-    density=lambda vp: 2.0,
+    density=density,
     optimizer_args={
         "popsize": 10,
         "maxiter": 100,
@@ -56,7 +56,7 @@
 # Run inversion
 # See stochopy's documentation for optimizer options <https://keurfonluu.github.io/stochopy/>
 res = model.invert(curves)
-res = res.threshold(0.1)
+res = res.threshold(0.02)
 
 # Plot results
 fig, ax = plt.subplots(1, 3, figsize=(15, 6))
diff --git a/README.rst b/README.rst
index f95c149..d9319c3 100644
--- a/README.rst
+++ b/README.rst
@@ -97,20 +97,20 @@ Expected output:
 .. code-block::
 
     --------------------------------------------------------------------------------
-    Best model out of 1000 models (1 run)
-
+    Best model out of 501 models (1 run)
+  
     Velocity model                                    Model parameters
     ----------------------------------------          ------------------------------
              d        vp        vs       rho                   d        vs        nu
           [km]    [km/s]    [km/s]   [g/cm3]                [km]    [km/s]       [-]
     ----------------------------------------          ------------------------------
-        0.0296    0.5033    0.2055    2.0000              0.0296    0.2055    0.4000
-        1.0000    1.8191    1.0080    2.0000                   -    1.0080    0.2785
+       0.0298    0.5033    0.2055    2.0000              0.0298    0.2055    0.4000
+       1.0000    2.0586    0.9935    2.0000                   -    0.9935    0.3482
     ----------------------------------------          ------------------------------
-
+  
     Number of layers: 2
     Number of parameters: 5
-    Best model misfit: 0.0153
+    Best model misfit: 0.0038
     --------------------------------------------------------------------------------
 
 Contributing
diff --git a/evodcinv/VERSION b/evodcinv/VERSION
index 8f9174b..e3a4f19 100644
--- a/evodcinv/VERSION
+++ b/evodcinv/VERSION
@@ -1,1 +1,1 @@
-2.1.2
\ No newline at end of file
+2.2.0
\ No newline at end of file
diff --git a/evodcinv/__init__.py b/evodcinv/__init__.py
index a880586..a185a1e 100644
--- a/evodcinv/__init__.py
+++ b/evodcinv/__init__.py
@@ -1,3 +1,4 @@
+from . import factory
 from .__about__ import __version__
 from ._curve import Curve
 from ._io import read
@@ -10,6 +11,7 @@
     "EarthModel",
     "Layer",
     "InversionResult",
+    "factory",
     "read",
     "__version__",
 ]
diff --git a/evodcinv/_curve.py b/evodcinv/_curve.py
index ca28675..e1fe74c 100644
--- a/evodcinv/_curve.py
+++ b/evodcinv/_curve.py
@@ -31,7 +31,7 @@ def __init__(
         weight : scalar, optional, default 1.0
             Overall weight applied to the misfit error associated to this data set.
         uncertainties : scalar, array_like or None, optional, default None
-            Uncertainties associated to data points. If None, error will be normalized by `data`.
+            Uncertainties associated to data points.
 
         """
         if len(period) != len(data):
@@ -44,6 +44,12 @@ def __init__(
             raise ValueError()
         if not is_sorted(period):
             raise ValueError()
+        if (
+            uncertainties is not None
+            and np.ndim(uncertainties) == 1
+            and len(uncertainties) != len(data)
+        ):
+            raise ValueError()
 
         self._period = np.asarray(period)
         self._data = np.asarray(data)
diff --git a/evodcinv/_model.py b/evodcinv/_model.py
index dbe566d..985f803 100644
--- a/evodcinv/_model.py
+++ b/evodcinv/_model.py
@@ -3,6 +3,7 @@
 from disba._common import ifunc
 from stochopy.optimize import minimize
 
+from . import factory
 from ._common import itype
 from ._curve import Curve
 from ._helpers import nafe_drake
@@ -284,12 +285,7 @@ def invert(self, curves, maxrun=1, split_results=False):
                     f"Option `increasing_velocity` is not compatible yet with optimizer `{method}`."
                 )
 
-            def constraint(x):
-                vs = x[self.n_layers - 1 : 2 * self.n_layers - 1]
-
-                return 0.0 if (vs[1:] >= vs[:-1]).all() else np.Inf
-
-            self._configuration["extra_terms"].append(constraint)
+            self._configuration["extra_terms"].append(factory.increasing_velocity)
 
         # Run maxrun inversion
         results = []
@@ -337,7 +333,9 @@ def constraint(x):
             with ProgressBar(prefix, max=maxiter) as bar:
 
                 def callback(X, res):
-                    bar.misfit = res.fun
+                    bar.misfit = (
+                        f"{res.fun:.4f}" if res.fun >= 1.0e-4 else f"{res.fun:.4e}"
+                    )
                     bar.next()
 
                 x = minimize(
@@ -452,10 +450,9 @@ def _misfit_function(self, x, curves):
                 n = len(dcalc)
                 if n > 0:
                     sigma = (
-                        curve.data[:n]
-                        if curve.uncertainties is None
-                        else curve.uncertainties[:n]
+                        curve.uncertainties if curve.uncertainties is not None else 1.0
                     )
+                    sigma = sigma[:n] if np.ndim(sigma) == 1 else sigma
                     error += curve.weight * misfit((dcalc - curve.data[:n]) / sigma)
                     weights_sum += curve.weight
 
diff --git a/evodcinv/_progress.py b/evodcinv/_progress.py
index 0ef0948..47808b2 100644
--- a/evodcinv/_progress.py
+++ b/evodcinv/_progress.py
@@ -9,7 +9,7 @@ def __init__(self, *args, **kwargs):
 
         self.width = 20
         self.suffix = (
-            "%(percent)3d%% [%(elapsed_td)s / %(eta_td)s] - Misfit: %(misfit).4f"
+            "%(percent)3d%% [%(elapsed_td)s / %(eta_td)s] - Misfit: %(misfit)s"
         )
         self.check_tty = False
         self.misfit = np.Inf
diff --git a/evodcinv/_result.py b/evodcinv/_result.py
index ce7749c..5ca6745 100644
--- a/evodcinv/_result.py
+++ b/evodcinv/_result.py
@@ -66,9 +66,10 @@ def __repr__(self):
         out += [f"{40 * '-'}{'':10}{30 * '-'}\n"]
 
         # Misc
+        misfit = f"{self.misfit:.4f}" if self.misfit >= 1.0e-4 else f"{self.misfit:.4e}"
         out += [f"Number of layers: {n_layers}"]
         out += [f"Number of parameters: {n_layers * 3 - 1}"]
-        out += [f"Best model misfit: {self.misfit:>.4f}"]
+        out += [f"Best model misfit: {misfit}"]
 
         out += [f"{80 * '-'}"]
 
diff --git a/evodcinv/factory/__init__.py b/evodcinv/factory/__init__.py
new file mode 100644
index 0000000..c4a5cbe
--- /dev/null
+++ b/evodcinv/factory/__init__.py
@@ -0,0 +1,7 @@
+from ._extra_terms import increasing_velocity, prior, smooth
+
+__all__ = [
+    "increasing_velocity",
+    "prior",
+    "smooth",
+]
diff --git a/evodcinv/factory/_extra_terms.py b/evodcinv/factory/_extra_terms.py
new file mode 100644
index 0000000..a9e0198
--- /dev/null
+++ b/evodcinv/factory/_extra_terms.py
@@ -0,0 +1,82 @@
+import numpy as np
+
+
+def increasing_velocity(x, penalty=np.inf):
+    """
+    Strictly increasing velocity models.
+
+    Parameters
+    ----------
+    x : array_like
+        Model parameter vector.
+    penalty : scalar, optional, default Inf
+        Penalty applied to misfit of models with velocities that are not strictly increasing.
+
+    Returns
+    -------
+    scalar
+        Misfit value.
+
+    """
+    n_layers = (len(x) + 1) // 3
+    vs = x[n_layers - 1 : 2 * n_layers - 1]
+
+    return 0.0 if (vs[1:] >= vs[:-1]).all() else penalty
+
+
+def prior(x, depth, velocity_s, uncertainties=None, alpha=1.0e-3):
+    """
+    Prior S-wave velocity model.
+
+    Parameters
+    ----------
+    x : array_like
+        Model parameter vector.
+    depth : array_like
+        Prior model's depths (in km).
+    velocity_s : array_like.
+        Prior model's S-wave velocities (in km/s).
+    uncertainties : scalar, array_like or None, optional, default None
+        Uncertainties associated to prior model S-wave velocities (in km/s).
+    alpha : scalar, optional, default 1.0e-3
+        Regularization factor.
+
+    Returns
+    -------
+    scalar
+        Misfit value.
+
+    """
+    n_layers = (len(x) + 1) // 3
+    d = x[: n_layers - 1]
+    z = np.insert(d.cumsum(), 0, 0.0)
+    vs = x[n_layers - 1 : 2 * n_layers - 1]
+    vprior = np.interp(z, depth, velocity_s)
+
+    sigma = uncertainties if uncertainties is not None else 1.0
+    sigma = np.interp(z, depth, uncertainties) if np.ndim(sigma) == 1 else sigma
+
+    return alpha * np.square((vs - vprior) / sigma).sum()
+
+
+def smooth(x, alpha=1.0e-3):
+    """
+    Smooth velocity models.
+
+    Parameters
+    ----------
+    x : array_like
+        Model parameter vector.
+    alpha : scalar, optional, default 1.0e-3
+        Regularization factor.
+
+    Returns
+    -------
+    scalar
+        Misfit value.
+
+    """
+    n_layers = (len(x) + 1) // 3
+    vs = x[n_layers - 1 : 2 * n_layers - 1]
+
+    return alpha * n_layers * np.square(vs[:-2] - 2.0 * vs[1:-1] + vs[2:]).sum()


Test patch:
diff --git a/tests/test_model.py b/tests/test_model.py
index be58694..78e1f85 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -6,7 +6,7 @@
 
 
 @pytest.mark.parametrize(
-    "xref, increasing_velocity", [(8.90561595, False), (8.88578982, True)]
+    "xref, increasing_velocity", [(8.90561595, False), (9.13114157, True)]
 )
 def test_invert(xref, increasing_velocity):
     model = copy.deepcopy(helpers.model)