Merge pull request #20 from keurfonluu/devel

Handle water layer

evodcinv/VERSION

@@ -1 +1 @@
-2.2.1
+2.2.2

evodcinv/_layer.py

@@ -12,6 +12,11 @@ def __init__(self, thickness, velocity_s, poisson=None):
             Layer thickness search boundary (in km).
         velocity_s : scalar or array_like
             Layer S-wave velocity search boundary (in km/s).
+            For the first layer, it can be lower or equal to zero for a water layer which P-wave velocity is determined as follows:
+
+             - Vs = 0, Vp = 1.5 km/s
+             - Vs < 0, Vp = -Vs km/s
+
         poisson : scalar, array_like or None, optional, default None
             Layer Poisson's ratio search boundary.
 

evodcinv/_model.py

@@ -79,6 +79,9 @@ def add(self, layer):
         if not isinstance(layer, Layer):
             raise TypeError()
 
+        if self.n_layers and layer.velocity_s[0] <= 0.0:
+            raise ValueError("only the first layer can be a water layer.")
+
         self.layers.append(layer)
 
     def pop(self):
@@ -287,7 +290,7 @@ def invert(self, curves, maxrun=1, split_results=False):
 
             self._configuration["extra_terms"].append(factory.increasing_velocity)
 
-        # Run maxrun inversion
+        # Run maxrun inversions
         results = []
 
         for i in range(maxrun):
@@ -389,12 +392,18 @@ def transform(self, x):
             Velocity model with shape (n_layers, 4).
 
         """
-        thickness = x[: self.n_layers - 1]
-        velocity_s = x[self.n_layers - 1 : 2 * self.n_layers - 1]
-        poisson = x[2 * self.n_layers - 1 :]
+        thickness = x[: self.n_layers - 1].copy()
+        velocity_s = x[self.n_layers - 1 : 2 * self.n_layers - 1].copy()
+        poisson = x[2 * self.n_layers - 1 :].copy()
         velocity_p = self._get_velocity_p(velocity_s, poisson)
         density = self._get_density(velocity_p)
 
+        # Handle water layer
+        if velocity_s[0] <= 0.0:
+            velocity_p[0] = abs(velocity_s[0]) if velocity_s[0] < 0.0 else 1.5
+            velocity_s[0] = 0.0
+            density[0] = 1.0
+
         return np.column_stack(
             (np.append(thickness, 1.0), velocity_p, velocity_s, density)
         )

evodcinv/_result.py

@@ -68,7 +68,7 @@ def __repr__(self):
         # 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"Number of parameters: {(~(self.xs == self.xs[0]).all(axis=0)).sum()}"]
         out += [f"Best model misfit: {misfit}"]
 
         out += [f"{80 * '-'}"]