Make discretize optional for tests

.github/environment_ci.yml

@@ -15,7 +15,6 @@ dependencies:
   - pydata-sphinx-theme==0.15.4
   - sphinx-gallery>=0.1.13
   - numpydoc>=1.5
-  - discretize
   - jupyter
   - graphviz
   - pillow

.github/workflows/test_with_conda.yml

@@ -18,10 +18,12 @@ jobs:
       fail-fast: false
       matrix:
         os: ["ubuntu-latest"]
-        python-version: ["3.10", "3.11", "3.12"]
+        python-version: ["3.10", "3.11", "3.12", "3.13"]
 
     steps:
     - uses: actions/checkout@v4
+      with:
+        fetch-depth: 0
 
     - name: Setup Conda
       uses: conda-incubator/setup-miniconda@v3
@@ -31,6 +33,25 @@ jobs:
         activate-environment: geoana-test
         python-version: ${{ matrix.python-version }}
 
+    - name: Install numba
+      if: matrix.python-version != '3.13'
+      # Numba doesn't work on python 3.13 just yet:
+      run: |
+        conda install --yes -c conda-forge numba
+
+      # Install discretize from it's repo until wheels/conda-forge built against numpy2.0 available:
+    - name: Pull Discretize
+      uses: actions/checkout@v4
+      with:
+        fetch-depth: 0
+        repository: 'simpeg/discretize'
+        ref: 'main'
+        path: 'discretize'
+
+    - name: Install discretize
+      run: |
+        pip install ./discretize
+
     - name: Conda information
       run: |
         conda info
@@ -39,11 +60,11 @@ jobs:
 
     - name: Install Our Package
       run: |
-        pip install . --config-settings=setup-args="-Dwith_extensions=true"
+        pip install --no-build-isolation --editable . --config-settings=setup-args="-Dwith_extensions=true"
 
     - name: Run Tests
       run: |
-        pytest --cov-config=.coveragerc --cov=geoana --cov-report=xml -s -v -W ignore::DeprecationWarning
+        pytest tests --cov-config=.coveragerc --cov=geoana --cov-report=xml -s -v -W ignore::DeprecationWarning
 
     - name: "Upload coverage to Codecov"
       if: matrix.python-version == '3.11'

geoana/gravity.py

@@ -555,7 +555,7 @@ def gravitational_gradient(self, xyz):
         ind1 = r == self.radius
         g_tens[ind0] = super().gravitational_gradient(xyz[ind0])
         g_tens[~ind0] = -G * 4 / 3 * np.pi * self.rho * np.eye(3)
-        g_tens[ind1] = np.NaN
+        g_tens[ind1] = np.nan
         return g_tens
 
 

tests/test_em_fdem_electric_dipole.py

@@ -4,7 +4,6 @@
 
 from scipy.constants import mu_0, epsilon_0
 from geoana.em import fdem
-import discretize
 
 # from SimPEG.EM import FDEM
 # from SimPEG import Maps
@@ -31,7 +30,7 @@ def E_from_EDWS(
     epsilon = epsilon_0*epsr
     sig_hat = sig + 1j*fdem.omega(f)*epsilon
 
-    XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3)
+    XYZ = np.atleast_2d(XYZ)
 
     dx = XYZ[:, 0] - srcLoc[0]
     dy = XYZ[:, 1] - srcLoc[1]
@@ -133,7 +132,7 @@ def electric_dipole_e(self, orientation):
         x = np.linspace(-20., 20., 50)
         y = np.linspace(-30., 30., 50)
         z = np.linspace(-40., 40., 50)
-        xyz = discretize.utils.ndgrid([x, y, z])
+        xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3)
 
         extest, eytest, eztest = E_from_EDWS(
             xyz, edws.location, edws.sigma, edws.frequency,
@@ -170,7 +169,7 @@ def test_electric_dipole_tilted_e(self):
         y = np.linspace(-30., 30., 50)
         z = np.linspace(-40., 40., 50)
 
-        xyz = discretize.utils.ndgrid([x, y, z])
+        xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3)
 
         extest0, eytest0, eztest0 = E_from_EDWS(
             xyz, edws.location, edws.sigma, edws.frequency,

tests/test_em_fdem_layered.py

@@ -1,4 +1,6 @@
 import unittest
+from tabnanny import check
+
 import pytest
 import numpy as np
 from numpy.testing import assert_allclose
@@ -8,11 +10,13 @@
 from geoana.kernels.tranverse_electric_reflections import (
     rTE_forward, rTE_gradient, _rTE_forward, _rTE_gradient
 )
-import discretize
 from scipy.special import iv, kv
 from geoana.em.fdem.base import sigma_hat
 
-from discretize.tests import check_derivative
+try:
+    from discretize.tests import check_derivative
+except ImportError:
+    check_derivative = None
 
 
 class TestHalfSpace:
@@ -180,7 +184,7 @@ def test_magnetic_field_errors(self):
             x = np.linspace(-20., 20., 50)
             y = np.linspace(-30., 30., 50)
             z = np.linspace(-40., 40., 50)
-            xyz = discretize.utils.ndgrid([x, y, z])
+            xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3)
             mag_layer.magnetic_field(xyz)
 
     def test_magnetic_field(self):
@@ -279,64 +283,65 @@ def test_magnetic_dipole(self):
             assert_allclose(em_layer_sec, em_half_sec, rtol=1e-05, atol=1e-08)
 
 
-class TestrTEGradient(unittest.TestCase):
-    def test_rTE_jacobian(self):
-        """Test to make sure numpy and compiled give same results"""
-        n_layer = 11
-        n_frequency = 5
-        n_lambda = 8
-        frequencies = np.logspace(1, 4, 5)
-        thicknesses = np.ones(n_layer-1)
-        lamb = np.logspace(0, 3, n_lambda)
-        np.random.seed(0)
-        sigma = 1E-1*(1 + 0.1 * np.random.rand(n_layer, n_frequency))
-        mu = mu_0 * (1 + 0.1 * np.random.rand(n_layer, n_frequency))
-        dmu = mu_0 * 0.1 * np.random.rand(n_layer, n_frequency)
-
-        def rte_sigma(x):
-            sigma = x.reshape(n_layer, n_frequency)
-            rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses)
-
-            J_sigma, _, _ = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses)
-
-            def J(y):
-                y = y.reshape(n_layer, n_frequency)
-                # do summation over layers, broadcast over frequencies only
-                # equivalent to:
-                # out = np.empty_like(rTE)
-                # for i in range(n_frequency):
-                #     out[i] = J_sigma[:, i, :]).T@y[:, i]
-                return np.einsum('i...k,i...', J_sigma, y)
-            return rTE, J
-
-        def rte_h(x):
-            thicknesses = x
-            rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses)
-            _, J_h, _ = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses)
-
-            def J(y):
-                #(J_h.T@y).T
-                out = np.einsum('i...,i', J_h, y)
-                # do sum over n_layer, broadcast over all others
-                return out
-
-            return rTE, J
-
-        def rte_mu(x):
-            mu = x.reshape(n_layer, n_frequency)
-            rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses)
-
-            _, _, J_mu = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses)
-
-            def J(y):
-                y = y.reshape(n_layer, n_frequency)
-                # do summation over layers, broadcast over frequencies only
-                return np.einsum('i...k,i...', J_mu, y)
-            return rTE, J
-
-        self.assertTrue(check_derivative(rte_sigma, sigma.reshape(-1), num=4, plotIt=False))
-        self.assertTrue(check_derivative(rte_h, thicknesses, num=4, plotIt=False))
-        self.assertTrue(check_derivative(rte_mu, mu.reshape(-1), dx=dmu.reshape(-1), num=4, plotIt=False))
+if check_derivative is not None:
+    class TestrTEGradient(unittest.TestCase):
+        def test_rTE_jacobian(self):
+            """Test to make sure numpy and compiled give same results"""
+            n_layer = 11
+            n_frequency = 5
+            n_lambda = 8
+            frequencies = np.logspace(1, 4, 5)
+            thicknesses = np.ones(n_layer-1)
+            lamb = np.logspace(0, 3, n_lambda)
+            np.random.seed(0)
+            sigma = 1E-1*(1 + 0.1 * np.random.rand(n_layer, n_frequency))
+            mu = mu_0 * (1 + 0.1 * np.random.rand(n_layer, n_frequency))
+            dmu = mu_0 * 0.1 * np.random.rand(n_layer, n_frequency)
+
+            def rte_sigma(x):
+                sigma = x.reshape(n_layer, n_frequency)
+                rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses)
+
+                J_sigma, _, _ = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses)
+
+                def J(y):
+                    y = y.reshape(n_layer, n_frequency)
+                    # do summation over layers, broadcast over frequencies only
+                    # equivalent to:
+                    # out = np.empty_like(rTE)
+                    # for i in range(n_frequency):
+                    #     out[i] = J_sigma[:, i, :]).T@y[:, i]
+                    return np.einsum('i...k,i...', J_sigma, y)
+                return rTE, J
+
+            def rte_h(x):
+                thicknesses = x
+                rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses)
+                _, J_h, _ = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses)
+
+                def J(y):
+                    #(J_h.T@y).T
+                    out = np.einsum('i...,i', J_h, y)
+                    # do sum over n_layer, broadcast over all others
+                    return out
+
+                return rTE, J
+
+            def rte_mu(x):
+                mu = x.reshape(n_layer, n_frequency)
+                rTE = rTE_forward(frequencies, lamb, sigma, mu, thicknesses)
+
+                _, _, J_mu = rTE_gradient(frequencies, lamb, sigma, mu, thicknesses)
+
+                def J(y):
+                    y = y.reshape(n_layer, n_frequency)
+                    # do summation over layers, broadcast over frequencies only
+                    return np.einsum('i...k,i...', J_mu, y)
+                return rTE, J
+
+            self.assertTrue(check_derivative(rte_sigma, sigma.reshape(-1), num=4, plotIt=False))
+            self.assertTrue(check_derivative(rte_h, thicknesses, num=4, plotIt=False))
+            self.assertTrue(check_derivative(rte_mu, mu.reshape(-1), dx=dmu.reshape(-1), num=4, plotIt=False))
 
 
 class TestCompiledVsNumpy(unittest.TestCase):

tests/test_em_fdem_magnetic_dipole.py

@@ -1,7 +1,5 @@
 import unittest
 import numpy as np
-import discretize
-
 from scipy.constants import mu_0, epsilon_0
 from geoana.em import fdem
 
@@ -26,7 +24,7 @@ def H_from_MagneticDipoleWholeSpace(
 
     omega = lambda f: 2*np.pi*f
 
-    XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3)
+    XYZ = np.atleast_2d(XYZ)
     # Check
 
     dx = XYZ[:, 0]-srcLoc[0]
@@ -88,7 +86,7 @@ def E_from_MagneticDipoleWholeSpace(
 
     omega = lambda f: 2 * np.pi * f
 
-    XYZ = discretize.utils.as_array_n_by_dim(XYZ, 3)
+    XYZ = np.atleast_2d(XYZ)
 
     dx = XYZ[:, 0]-srcLoc[0]
     dy = XYZ[:, 1]-srcLoc[1]
@@ -188,7 +186,7 @@ def magnetic_dipole_b(self, orientation):
         x = np.linspace(-20., 20., 50)
         y = np.linspace(-30., 30., 50)
         z = np.linspace(-40., 40., 50)
-        xyz = discretize.utils.ndgrid([x, y, z])
+        xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3)
 
         # srcLoc, obsLoc, component, orientation='Z', moment=1., mu=mu_0
 
@@ -221,7 +219,7 @@ def magnetic_dipole_e(self, orientation):
         x = np.linspace(-20., 20., 50)
         y = np.linspace(-30., 30., 50)
         z = np.linspace(-40., 40., 50)
-        xyz = discretize.utils.ndgrid([x, y, z])
+        xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3)
 
         extest, eytest, eztest = E_from_MagneticDipoleWholeSpace(
             xyz, mdws.location, mdws.sigma, mdws.frequency,
@@ -258,7 +256,7 @@ def test_magnetic_dipole_tilted_b(self):
         y = np.linspace(-30., 30., 50)
         z = np.linspace(-40., 40., 50)
 
-        xyz = discretize.utils.ndgrid([x, y, z])
+        xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3)
 
         bxtest0, bytest0, bztest0 = B_from_MagneticDipoleWholeSpace(
             xyz, mdws.location, mdws.sigma, mdws.frequency,
@@ -312,7 +310,7 @@ def test_magnetic_dipole_tilted_e(self):
         y = np.linspace(-30., 30., 10)
         z = np.linspace(-40., 40., 10)
 
-        xyz = discretize.utils.ndgrid([x, y, z])
+        xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3)
 
         extest0, eytest0, eztest0 = E_from_MagneticDipoleWholeSpace(
             xyz, mdws.location, mdws.sigma, mdws.frequency,

tests/test_em_fdem_planewave.py

@@ -1,6 +1,5 @@
 import pytest
 import numpy as np
-import discretize
 from geoana.em import fdem
 from scipy.constants import mu_0, epsilon_0
 
@@ -48,7 +47,7 @@ def test_electric_field():
     x = np.linspace(-20., 20., 50)
     y = np.linspace(-30., 30., 50)
     z = np.linspace(-40., 40., 50)
-    xyz = discretize.utils.ndgrid([x, y, z])
+    xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3)
     z = xyz[:, 2]
 
     kz = np.outer(k, z)
@@ -90,7 +89,7 @@ def test_current_density():
     x = np.linspace(-20., 20., 50)
     y = np.linspace(-30., 30., 50)
     z = np.linspace(-40., 40., 50)
-    xyz = discretize.utils.ndgrid([x, y, z])
+    xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3)
     z = xyz[:, 2]
 
     kz = np.outer(k, z)
@@ -130,7 +129,7 @@ def test_magnetic_field():
     x = np.linspace(-20., 20., 50)
     y = np.linspace(-30., 30., 50)
     z = np.linspace(-40., 40., 50)
-    xyz = discretize.utils.ndgrid([x, y, z])
+    xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3)
     z = xyz[:, 2]
 
     kz = z * k
@@ -171,7 +170,7 @@ def test_magnetic_flux_density():
     x = np.linspace(-20., 20., 50)
     y = np.linspace(-30., 30., 50)
     z = np.linspace(-40., 40., 50)
-    xyz = discretize.utils.ndgrid([x, y, z])
+    xyz = np.stack(np.meshgrid(x, y, z), axis=-1).reshape(-1, 3)
     z = xyz[:, 2]
 
     kz = z * k