Merge pull request #17 from keurfonluu/devel

Add factory module for extra misfit terms

.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))

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

evodcinv/VERSION

@@ -1 +1 @@
-2.1.2
+2.2.0

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__",
 ]

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)

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
 

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

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 * '-'}"]
 

evodcinv/factory/__init__.py

@@ -0,0 +1,7 @@
+from ._extra_terms import increasing_velocity, prior, smooth
+
+__all__ = [
+    "increasing_velocity",
+    "prior",
+    "smooth",
+]

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()

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)