Add linearity tests for forward models

test/test_linearity.py

@@ -0,0 +1,372 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Test linearity of the forward models.
+
+The forward models are linear on the sources: the field generated by a set of
+sources should be equal to the sum of the fields generated by each source
+separately.
+"""
+
+import bordado as bd
+import boule
+import numpy as np
+import numpy.testing as npt
+import pytest
+
+from harmonica import (
+    dipole_magnetic,
+    point_gravity,
+    prism_gravity,
+    prism_magnetic,
+    tesseroid_gravity,
+)
+from harmonica._forward.dipole import VALID_FIELDS as DIPOLE_FIELDS
+from harmonica._forward.prisms.magnetic import VALID_FIELDS as PRISM_MAGNETIC_FIELDS
+
+GRAVITY_FIELDS = (
+    "potential",
+    "g_e",
+    "g_n",
+    "g_z",
+    "g_ee",
+    "g_nn",
+    "g_zz",
+    "g_en",
+    "g_ez",
+    "g_nz",
+)
+
+
+class TestPointGravityLinearity:
+    """
+    Test linearity of ``point_gravity``.
+    """
+
+    @pytest.fixture()
+    def sample_points(self):
+        """
+        Return four sample point masses in Cartesian coordinates.
+        """
+        points = (
+            np.array([-50.0, 0.0, 55.0, 10.0]),
+            np.array([20.0, 0.0, -35.0, 15.0]),
+            np.array([-100.0, -150.0, -50.0, -20.0]),
+        )
+        return points
+
+    @pytest.fixture()
+    def sample_masses(self):
+        """
+        Return masses for the sample points, including negative ones.
+        """
+        return np.array([2.0e3, -3.0e3, 5.0e3, -1.5e3])
+
+    @pytest.fixture()
+    def sample_coordinates(self):
+        """
+        Return a grid of observation points above the sources.
+        """
+        return bd.grid_coordinates(
+            region=(-100, 100, -100, 100), spacing=20, non_dimensional_coords=10
+        )
+
+    @pytest.mark.parametrize("field", GRAVITY_FIELDS)
+    def test_linearity(self, sample_coordinates, sample_points, sample_masses, field):
+        """
+        Compare the field of all points vs. the sum of their individual fields.
+        """
+        result = point_gravity(sample_coordinates, sample_points, sample_masses, field)
+        expected = sum(
+            point_gravity(
+                sample_coordinates,
+                tuple(p[i : i + 1] for p in sample_points),
+                sample_masses[i : i + 1],
+                field,
+            )
+            for i in range(sample_masses.size)
+        )
+        npt.assert_allclose(result, expected)
+
+    @pytest.mark.parametrize("field", ["potential", "g_z"])
+    def test_linearity_spherical(self, sample_masses, field):
+        """
+        Compare joint vs. summed individual fields in spherical coordinates.
+        """
+        radius = boule.WGS84.mean_radius
+        points = (
+            np.array([-50.0, 0.0, 55.0, 10.0]),
+            np.array([20.0, 0.0, -35.0, 15.0]),
+            radius - np.array([1e3, 2e3, 5e3, 10e3]),
+        )
+        coordinates = bd.grid_coordinates(
+            region=(-100, 100, -70, 70), spacing=20, non_dimensional_coords=radius
+        )
+        result = point_gravity(
+            coordinates, points, sample_masses, field, coordinate_system="spherical"
+        )
+        expected = sum(
+            point_gravity(
+                coordinates,
+                tuple(p[i : i + 1] for p in points),
+                sample_masses[i : i + 1],
+                field,
+                coordinate_system="spherical",
+            )
+            for i in range(sample_masses.size)
+        )
+        npt.assert_allclose(result, expected)
+
+
+class TestDipoleMagneticLinearity:
+    """
+    Test linearity of ``dipole_magnetic``.
+    """
+
+    @pytest.fixture()
+    def sample_dipoles(self):
+        """
+        Return four sample dipoles.
+        """
+        dipoles = (
+            np.array([-10.0, 0.0, 10.0, 5.0]),
+            np.array([2.0, 0.0, -1.0, 1.5]),
+            np.array([-10.0, -15.0, -5.0, -2.0]),
+        )
+        return dipoles
+
+    @pytest.fixture()
+    def sample_magnetic_moments(self):
+        """
+        Return sample magnetic moment vectors for the dipoles.
+        """
+        magnetic_moments = (
+            np.array([1.0, -1.0, 3.0, -1.5]),
+            np.array([1.0, -0.5, -1.0, 2.0]),
+            np.array([1.0, 2.0, -3.0, -2.0]),
+        )
+        return magnetic_moments
+
+    @pytest.fixture()
+    def sample_coordinates(self):
+        """
+        Return a grid of observation points above the dipoles.
+        """
+        return bd.grid_coordinates(
+            region=(-20, 20, -20, 20), spacing=4, non_dimensional_coords=10
+        )
+
+    @pytest.mark.parametrize("field", DIPOLE_FIELDS)
+    def test_linearity(
+        self, sample_coordinates, sample_dipoles, sample_magnetic_moments, field
+    ):
+        """
+        Compare the field of all dipoles vs. the sum of their individual fields.
+        """
+        result = dipole_magnetic(
+            sample_coordinates, sample_dipoles, sample_magnetic_moments, field=field
+        )
+        n_dipoles = sample_dipoles[0].size
+        expected = sum(
+            np.asarray(
+                dipole_magnetic(
+                    sample_coordinates,
+                    tuple(d[i : i + 1] for d in sample_dipoles),
+                    tuple(m[i : i + 1] for m in sample_magnetic_moments),
+                    field=field,
+                )
+            )
+            for i in range(n_dipoles)
+        )
+        npt.assert_allclose(result, expected)
+
+
+class TestPrismGravityLinearity:
+    """
+    Test linearity of ``prism_gravity``.
+    """
+
+    @pytest.fixture()
+    def sample_prisms(self):
+        """
+        Return four non-overlapping sample prisms.
+        """
+        prisms = np.array(
+            [
+                [-100.0, -50.0, -100.0, -50.0, -30.0, -10.0],
+                [-50.0, 50.0, 50.0, 100.0, -50.0, -20.0],
+                [50.0, 100.0, -100.0, 50.0, -20.0, -5.0],
+                [-100.0, 50.0, -50.0, 50.0, -60.0, -40.0],
+            ]
+        )
+        return prisms
+
+    @pytest.fixture()
+    def sample_densities(self):
+        """
+        Return densities for the sample prisms, including negative ones.
+        """
+        return np.array([2670.0, -1500.0, 3300.0, -2900.0])
+
+    @pytest.fixture()
+    def sample_coordinates(self):
+        """
+        Return a grid of observation points above the prisms.
+        """
+        return bd.grid_coordinates(
+            region=(-100, 100, -100, 100), spacing=20, non_dimensional_coords=10
+        )
+
+    @pytest.mark.parametrize("field", GRAVITY_FIELDS)
+    def test_linearity(
+        self, sample_coordinates, sample_prisms, sample_densities, field
+    ):
+        """
+        Compare the field of all prisms vs. the sum of their individual fields.
+        """
+        result = prism_gravity(
+            sample_coordinates, sample_prisms, sample_densities, field
+        )
+        expected = sum(
+            prism_gravity(
+                sample_coordinates,
+                sample_prisms[i : i + 1],
+                sample_densities[i : i + 1],
+                field,
+            )
+            for i in range(sample_densities.size)
+        )
+        npt.assert_allclose(result, expected)
+
+
+class TestPrismMagneticLinearity:
+    """
+    Test linearity of ``prism_magnetic``.
+    """
+
+    @pytest.fixture()
+    def sample_prisms(self):
+        """
+        Return four non-overlapping sample prisms.
+        """
+        prisms = np.array(
+            [
+                [-100.0, -50.0, -100.0, -50.0, -30.0, -10.0],
+                [-50.0, 50.0, 50.0, 100.0, -50.0, -20.0],
+                [50.0, 100.0, -100.0, 50.0, -20.0, -5.0],
+                [-100.0, 50.0, -50.0, 50.0, -60.0, -40.0],
+            ]
+        )
+        return prisms
+
+    @pytest.fixture()
+    def sample_magnetizations(self):
+        """
+        Return sample magnetization vectors for the prisms.
+        """
+        magnetizations = (
+            np.array([1.0, -1.0, 3.0, -1.5]),
+            np.array([1.0, -0.5, -1.0, 2.0]),
+            np.array([1.0, 2.0, -3.0, -2.0]),
+        )
+        return magnetizations
+
+    @pytest.fixture()
+    def sample_coordinates(self):
+        """
+        Return a grid of observation points above the prisms.
+        """
+        return bd.grid_coordinates(
+            region=(-100, 100, -100, 100), spacing=20, non_dimensional_coords=10
+        )
+
+    @pytest.mark.parametrize("field", PRISM_MAGNETIC_FIELDS)
+    def test_linearity(
+        self, sample_coordinates, sample_prisms, sample_magnetizations, field
+    ):
+        """
+        Compare the field of all prisms vs. the sum of their individual fields.
+        """
+        result = prism_magnetic(
+            sample_coordinates, sample_prisms, sample_magnetizations, field=field
+        )
+        expected = sum(
+            np.asarray(
+                prism_magnetic(
+                    sample_coordinates,
+                    sample_prisms[i : i + 1],
+                    tuple(m[i : i + 1] for m in sample_magnetizations),
+                    field=field,
+                )
+            )
+            for i in range(sample_prisms.shape[0])
+        )
+        npt.assert_allclose(result, expected)
+
+
+class TestTesseroidGravityLinearity:
+    """
+    Test linearity of ``tesseroid_gravity``.
+    """
+
+    @pytest.fixture()
+    def sample_tesseroids(self):
+        """
+        Return three non-overlapping sample tesseroids.
+        """
+        ellipsoid = boule.WGS84
+        top = ellipsoid.mean_radius
+        bottom = top - 1e3
+        tesseroids = np.array(
+            [
+                [-10.0, 10.0, -10.0, 10.0, bottom, top],
+                [20.0, 30.0, 20.0, 30.0, bottom - 1e3, top - 500],
+                [-50.0, -40.0, -30.0, -20.0, bottom - 2e3, top],
+            ]
+        )
+        return tesseroids
+
+    @pytest.fixture()
+    def sample_densities(self):
+        """
+        Return densities for the sample tesseroids, including a negative one.
+        """
+        return np.array([2670.0, -1500.0, 3300.0])
+
+    @pytest.fixture()
+    def sample_coordinates(self):
+        """
+        Return a grid of observation points above the tesseroids.
+        """
+        return bd.grid_coordinates(
+            region=(-60, 40, -40, 40),
+            spacing=20,
+            non_dimensional_coords=boule.WGS84.mean_radius + 100,
+        )
+
+    @pytest.mark.use_numba
+    @pytest.mark.parametrize("field", ["potential", "g_z"])
+    def test_linearity(
+        self, sample_coordinates, sample_tesseroids, sample_densities, field
+    ):
+        """
+        Compare the field of all tesseroids vs. the sum of their individual
+        fields.
+        """
+        result = tesseroid_gravity(
+            sample_coordinates, sample_tesseroids, sample_densities, field=field
+        )
+        expected = sum(
+            tesseroid_gravity(
+                sample_coordinates,
+                sample_tesseroids[i : i + 1],
+                sample_densities[i : i + 1],
+                field=field,
+            )
+            for i in range(sample_densities.size)
+        )
+        npt.assert_allclose(result, expected)