From 244632c866921825c70ac8dc891bf96674d92733 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 14 Apr 2026 13:02:55 -0700 Subject: [PATCH 1/5] Extend type hints for dia_arrays in mesh-based regs --- src/inversion_ideas/regularization/_mesh_based.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/inversion_ideas/regularization/_mesh_based.py b/src/inversion_ideas/regularization/_mesh_based.py index bda8c78..dce081d 100644 --- a/src/inversion_ideas/regularization/_mesh_based.py +++ b/src/inversion_ideas/regularization/_mesh_based.py @@ -229,7 +229,7 @@ def hessian(self, model: Model): # noqa: ARG002 ) @property - def weights_matrix(self) -> dia_array: + def weights_matrix(self) -> dia_array[np.float64]: """ Diagonal matrix with the square root of regularization weights on cells. """ @@ -243,7 +243,7 @@ def weights_matrix(self) -> dia_array: return diags_array(np.sqrt(cell_weights)) @property - def _volumes_sqrt_matrix(self) -> dia_array: + def _volumes_sqrt_matrix(self) -> dia_array[np.float64]: """ Diagonal matrix with the square root of cell volumes. """ @@ -427,7 +427,7 @@ def hessian(self, model: Model): # noqa: ARG002 ) @property - def weights_matrix(self) -> dia_array: + def weights_matrix(self) -> dia_array[np.float64]: """ Diagonal matrix with the square root of cell weights averaged on faces. """ @@ -441,7 +441,7 @@ def weights_matrix(self) -> dia_array: return diags_array(self._average_cells_to_faces @ np.sqrt(cell_weights)) @property - def _volumes_sqrt_matrix(self) -> dia_array: + def _volumes_sqrt_matrix(self) -> dia_array[np.float64]: """ Diagonal matrix with the square root of cell volumes averaged on faces. """ From 6852f7fff5ef82475f59a525a06e313f8ddd0c87 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 14 Apr 2026 13:03:15 -0700 Subject: [PATCH 2/5] Add a `n_active` property to all mesh-based regs --- src/inversion_ideas/regularization/_mesh_based.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/inversion_ideas/regularization/_mesh_based.py b/src/inversion_ideas/regularization/_mesh_based.py index dce081d..bed7423 100644 --- a/src/inversion_ideas/regularization/_mesh_based.py +++ b/src/inversion_ideas/regularization/_mesh_based.py @@ -31,6 +31,10 @@ class _MeshBasedRegularization(Objective): @property def n_params(self) -> int: + return self.n_active + + @property + def n_active(self) -> int: return int(np.sum(self.active_cells)) @property From e417d6110d9c22ac40ca4a3d595c7a82908794fe Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 14 Apr 2026 13:47:20 -0700 Subject: [PATCH 3/5] Replace `n_params` for `n_active` in `Smallness` and `Flatness` Make the distinction between number of model parameters and number of active cells in the mesh: the number of parameters is the size of the model that will get passed to the regularization. In case of joint inversions, the size of the model will not necessarily match the number of active cells. --- .../regularization/_mesh_based.py | 27 ++++++++++--------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/src/inversion_ideas/regularization/_mesh_based.py b/src/inversion_ideas/regularization/_mesh_based.py index bed7423..d81c6dd 100644 --- a/src/inversion_ideas/regularization/_mesh_based.py +++ b/src/inversion_ideas/regularization/_mesh_based.py @@ -59,20 +59,20 @@ def cell_weights( "It must be an array or a dictionary." ) raise TypeError(msg) - if isinstance(value, np.ndarray) and value.size != self.n_params: + if isinstance(value, np.ndarray) and value.size != self.n_active: msg = ( f"Invalid cell_weights array with '{value.size}' elements. " - f"It must have '{self.n_params}' elements, " + f"It must have '{self.n_active}' elements, " "equal to the number of active cells." ) raise ValueError(msg) if isinstance(value, dict): for key, array in value.items(): - if array.size != self.n_params: + if array.size != self.n_active: msg = ( f"Invalid cell_weights array '{key}' with " f"'{array.size}' elements. " - f"It must have '{self.n_params}' elements, " + f"It must have '{self.n_active}' elements, " "equal to the number of active cells." ) raise ValueError(msg) @@ -92,13 +92,14 @@ class Smallness(_MeshBasedRegularization): active_cells : (n_cells) array or None, optional Array full of bools that indicate the active cells in the mesh. It must have the same amount of elements as cells in the mesh. - cell_weights : (n_params) array or dict of (n_params) arrays or None, optional + cell_weights : (n_active) array or dict of (n_active) arrays or None, optional Array with cell weights. For multiple cell weights, pass a dictionary where keys are strings and values are the different weights arrays. If None, no cell weights are going to be used. - reference_model : (n_params) array or None, optional - Array with values for the reference model. + reference_model : (n_active) array or None, optional + Array with values for the reference model. It must have the same number of + elements as active cells in the mesh. Notes ----- @@ -161,13 +162,13 @@ def __init__( self.cell_weights = ( cell_weights if cell_weights is not None - else np.ones(self.n_params, dtype=np.float64) + else np.ones(self.n_active, dtype=np.float64) ) self.reference_model = ( reference_model if reference_model is not None - else np.zeros(self.n_params, dtype=np.float64) + else np.zeros(self.n_active, dtype=np.float64) ) self.set_name("s") @@ -270,12 +271,12 @@ class Flatness(_MeshBasedRegularization): active_cells : (n_cells) array or None, optional Array full of bools that indicate the active cells in the mesh. It must have the same amount of elements as cells in the mesh. - cell_weights : (n_params) array or dict of (n_params) arrays or None, optional + cell_weights : (n_active) array or dict of (n_active) arrays or None, optional Array with cell weights. For multiple cell weights, pass a dictionary where keys are strings and values are the different weights arrays. If None, no cell weights are going to be used. - reference_model : (n_params) array or None, optional + reference_model : (n_active) array or None, optional Array with values for the reference model. Notes @@ -350,13 +351,13 @@ def __init__( self.cell_weights = ( cell_weights if cell_weights is not None - else np.ones(self.n_params, dtype=np.float64) + else np.ones(self.n_active, dtype=np.float64) ) self.reference_model = ( reference_model if reference_model is not None - else np.zeros(self.n_params, dtype=np.float64) + else np.zeros(self.n_active, dtype=np.float64) ) self.set_name(direction) From 0af266be14f1fd98dd9ceee47baf483fe0bb0f8e Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 14 Apr 2026 13:48:42 -0700 Subject: [PATCH 4/5] Start adding tests for mesh-based regularizations --- tests/test_regularizations.py | 67 +++++++++++++++++++++++++++++++++-- 1 file changed, 65 insertions(+), 2 deletions(-) diff --git a/tests/test_regularizations.py b/tests/test_regularizations.py index 0bc2642..5fbce88 100644 --- a/tests/test_regularizations.py +++ b/tests/test_regularizations.py @@ -2,11 +2,12 @@ Test regularization classes. """ +import numpy as np import pytest from discretize.tensor_mesh import TensorMesh -from scipy.sparse import sparray +from scipy.sparse import dia_array, sparray -from inversion_ideas import Flatness +from inversion_ideas import Flatness, Smallness class TestBugfixFlatness: @@ -24,3 +25,65 @@ def mesh(self): def test_cell_gradient_type(self, mesh, direction): flatness = Flatness(mesh, direction=direction) assert isinstance(flatness._cell_gradient, sparray) + + +class TestSmallness: + """ + Test the :class:`inversion_ideas.Smallness` regularization class. + """ + + @pytest.fixture + def mesh(self): + hx = [(1.0, 10)] + h = [hx, hx, hx] + return TensorMesh(h, origin="CCN") + + @pytest.fixture + def active_cells(self, mesh: TensorMesh): + active_cells = np.ones(mesh.n_cells, dtype=bool) + _, _, z = mesh.cell_centers.T + active_cells[z > -1.0] = False + assert not active_cells.all() + return active_cells + + def test_smallness(self, mesh, active_cells): + n_active = active_cells.sum() + cell_weights = np.full(n_active, fill_value=0.1) + reference_model = np.full(n_active, 1e-8) + smallness = Smallness( + mesh, + active_cells=active_cells, + cell_weights=cell_weights, + reference_model=reference_model, + ) + + model = np.random.default_rng(seed=12312).uniform(size=n_active) + + # Test call + result = smallness(model) + assert np.isscalar(result) + expected = np.sum( + mesh.cell_volumes[active_cells] + * cell_weights + * (model - reference_model) ** 2 + ) + np.testing.assert_allclose(result, expected) + + # Test gradient + gradient = smallness.gradient(model) + expected = ( + 2 + * mesh.cell_volumes[active_cells] + * cell_weights + * (model - reference_model) + ) + assert gradient.size == model.size + np.testing.assert_allclose(gradient, expected) + + # Test hessian + hessian = smallness.hessian(model) + assert hessian.shape == (model.size, model.size) + assert isinstance(hessian, dia_array) + assert hessian.offsets == 0 # should be a diagonal matrix (only main diag) + expected_diagonal = 2 * mesh.cell_volumes[active_cells] * cell_weights + np.testing.assert_allclose(hessian.diagonal(), expected_diagonal) From 6cf0b92985dd2b71ec6dae972213473fc133b158 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 14 Apr 2026 13:52:01 -0700 Subject: [PATCH 5/5] Update `n_params` for `n_active` in `SparseSmallness` --- .../regularization/_mesh_based.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/inversion_ideas/regularization/_mesh_based.py b/src/inversion_ideas/regularization/_mesh_based.py index d81c6dd..5decaf9 100644 --- a/src/inversion_ideas/regularization/_mesh_based.py +++ b/src/inversion_ideas/regularization/_mesh_based.py @@ -502,21 +502,22 @@ class SparseSmallness(_MeshBasedRegularization): active_cells : (n_cells) array or None, optional Array full of bools that indicate the active cells in the mesh. It must have the same amount of elements as cells in the mesh. - cell_weights : (n_params) array or dict of (n_params) arrays or None, optional + cell_weights : (n_active) array or dict of (n_active) arrays or None, optional Array with cell weights. For multiple cell weights, pass a dictionary where keys are strings and values are the different weights arrays. If None, no cell weights are going to be used. - reference_model : (n_params) array, optional - Reference model used in the regularization. + reference_model : (n_active) array or None, optional + Array with values for the reference model. It must have the same number of + elements as active cells in the mesh. threshold : float, optional IRLS threshold. Symbolized with :math:`\epsilon` in Fournier and Oldenburg (2019). cooling_factor : float, optional Factor used to cool down the ``threshold`` when updating the IRLS. - model_previous : (n_params) array + model_previous : (n_params) array or None, optional Array with previous model in the iterations. This model is used to build the - ``R`` matrix. + ``R`` matrix. If None, an array full of zeros will be assigned. irls : bool, optional Flag to activate or deactivate IRLS. If False, the class would work as an L2 smallness term. If True, the R matrix will be built using the @@ -552,13 +553,13 @@ def __init__( self.cell_weights = ( cell_weights if cell_weights is not None - else np.ones(self.n_params, dtype=np.float64) + else np.ones(self.n_active, dtype=np.float64) ) self.reference_model = ( reference_model if reference_model is not None - else np.zeros(self.n_params, dtype=np.float64) + else np.zeros(self.n_active, dtype=np.float64) ) self.threshold = threshold self.cooling_factor = cooling_factor @@ -607,7 +608,7 @@ def R(self) -> dia_array: R matrix to approximate lp norm using Lawson's algorithm. """ if not self.irls: - return eye_array(self.n_params) + return eye_array(self.n_active) power = self.norm / 4 - 0.5 diagonal = (self.model_previous**2 + self.threshold**2) ** power return diags_array(diagonal)