From 0e00d0d443c784ba2b8caa9be62a5b360ff2a808 Mon Sep 17 00:00:00 2001 From: Keurfon Luu Date: Sat, 11 Nov 2023 10:44:28 +0100 Subject: [PATCH 1/4] handle water layer --- evodcinv/_layer.py | 5 +++++ evodcinv/_model.py | 20 ++++++++++++++++---- 2 files changed, 21 insertions(+), 4 deletions(-) diff --git a/evodcinv/_layer.py b/evodcinv/_layer.py index df4ed88..6031a9e 100644 --- a/evodcinv/_layer.py +++ b/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. diff --git a/evodcinv/_model.py b/evodcinv/_model.py index 2e0d5f3..e8d6f72 100644 --- a/evodcinv/_model.py +++ b/evodcinv/_model.py @@ -287,7 +287,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 +389,24 @@ 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 + cond = np.flatnonzero(velocity_s <= 0.0) + + if cond.size: + if cond[0] != 0 or cond.size != 1: + raise ValueError("only the first layer can be a water layer.") + + else: + 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) ) From 4132a71f728f7db22c3c97eb102d28159554dea7 Mon Sep 17 00:00:00 2001 From: Keurfon Luu Date: Sat, 11 Nov 2023 10:44:49 +0100 Subject: [PATCH 2/4] fix number of parameters --- evodcinv/_result.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/evodcinv/_result.py b/evodcinv/_result.py index 5ca6745..736dd66 100644 --- a/evodcinv/_result.py +++ b/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 * '-'}"] From c74d3dccc7fc62e13b3070a58c05df284eb69eeb Mon Sep 17 00:00:00 2001 From: Keurfon Luu Date: Sat, 11 Nov 2023 10:48:31 +0100 Subject: [PATCH 3/4] version bump --- evodcinv/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/evodcinv/VERSION b/evodcinv/VERSION index fae692e..7e541ae 100644 --- a/evodcinv/VERSION +++ b/evodcinv/VERSION @@ -1 +1 @@ -2.2.1 \ No newline at end of file +2.2.2 \ No newline at end of file From decc2180ae8d5e4584b51bb99c1021e44a07e319 Mon Sep 17 00:00:00 2001 From: Keurfon Luu Date: Sat, 11 Nov 2023 11:00:02 +0100 Subject: [PATCH 4/4] handle error when adding layer --- evodcinv/_model.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/evodcinv/_model.py b/evodcinv/_model.py index e8d6f72..418dfca 100644 --- a/evodcinv/_model.py +++ b/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): @@ -396,16 +399,10 @@ def transform(self, x): density = self._get_density(velocity_p) # Handle water layer - cond = np.flatnonzero(velocity_s <= 0.0) - - if cond.size: - if cond[0] != 0 or cond.size != 1: - raise ValueError("only the first layer can be a water layer.") - - else: - 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 + 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)