Repository: simpeg/simpeg
Pull request #1571 (2024-11-02 18:14:11)

Remove empymod dependency
#### Summary
`empymod` is a pretty heavy dependency for the functionality we are using from it. This removes the dependency and instead uses `libdlf` to get the filters, and a streamlined `get_splined_dlf_points` for our internal use.

#### PR Checklist
* [x] If this is a work in progress PR, set as a Draft PR
* [x] Linted my code according to the [style guides](https://docs.simpeg.xyz/latest/content/getting_started/contributing/code-style.html).
* [x] Added [tests](https://docs.simpeg.xyz/latest/content/getting_started/contributing/testing.html) to verify changes to the code.
* [x] Added necessary documentation to any new functions/classes following the
      expect [style](https://docs.simpeg.xyz/latest/content/getting_started/contributing/documentation.html).
* [x] Marked as ready for review (if this is was a draft PR), and converted 
      to a Pull Request
* [x] Tagged ``@simpeg/simpeg-developers`` when ready for review.

#### Additional information
I say `empymod` is a heavy dependency because it in turn relies on `numba`.


codecov[bot]: ## [Codecov](https://app.codecov.io/gh/simpeg/simpeg/pull/1571?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) Report
Attention: Patch coverage is `89.47368%` with `4 lines` in your changes missing coverage. Please review.
> Project coverage is 86.42%. Comparing base [(`30b3b42`)](https://app.codecov.io/gh/simpeg/simpeg/commit/30b3b420ea0eabace1420760dfb9d9e3e2ec4900?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) to head [(`530ca9f`)](https://app.codecov.io/gh/simpeg/simpeg/commit/530ca9f05a4987fa7b55ed4664109f5dfb3585ab?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg).
> Report is 1 commits behind head on main.

| [Files with missing lines](https://app.codecov.io/gh/simpeg/simpeg/pull/1571?dropdown=coverage&src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) | Patch % | Lines |
|---|---|---|
| [...mpeg/electromagnetics/time\_domain/simulation\_1d.py](https://app.codecov.io/gh/simpeg/simpeg/pull/1571?src=pr&el=tree&filepath=simpeg%2Felectromagnetics%2Ftime_domain%2Fsimulation_1d.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg#diff-c2ltcGVnL2VsZWN0cm9tYWduZXRpY3MvdGltZV9kb21haW4vc2ltdWxhdGlvbl8xZC5weQ==) | 81.81% | [1 Missing and 1 partial :warning: ](https://app.codecov.io/gh/simpeg/simpeg/pull/1571?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) |
| [simpeg/electromagnetics/utils/em1d\_utils.py](https://app.codecov.io/gh/simpeg/simpeg/pull/1571?src=pr&el=tree&filepath=simpeg%2Felectromagnetics%2Futils%2Fem1d_utils.py&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg#diff-c2ltcGVnL2VsZWN0cm9tYWduZXRpY3MvdXRpbHMvZW0xZF91dGlscy5weQ==) | 81.81% | [1 Missing and 1 partial :warning: ](https://app.codecov.io/gh/simpeg/simpeg/pull/1571?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg) |

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #1571      +/-   ##
==========================================
+ Coverage   86.25%   86.42%   +0.16%     
==========================================
  Files         405      406       +1     
  Lines       52162    52429     +267     
  Branches     4980     4990      +10     
==========================================
+ Hits        44994    45312     +318     
+ Misses       5761     5689      -72     
- Partials     1407     1428      +21     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/simpeg/simpeg/pull/1571?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=simpeg).


Patch:
diff --git a/environment.yml b/environment.yml
index 44a312829..5c017b7e1 100644
--- a/environment.yml
+++ b/environment.yml
@@ -10,7 +10,7 @@ dependencies:
   - matplotlib-base
   - discretize>=0.11
   - geoana>=0.7
-  - empymod>=2.0.0
+  - libdlf
 
 # solver
 # uncomment the next line if you are on an intel platform
@@ -33,6 +33,7 @@ dependencies:
   - sphinx-gallery>=0.1.13
   - sphinxcontrib-apidoc
   - pydata-sphinx-theme
+  - empymod>=2.0.0
   - nbsphinx
   - numpydoc
   - pillow
diff --git a/pyproject.toml b/pyproject.toml
index 024a31499..9076ef11f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -21,7 +21,7 @@ dependencies = [
     "matplotlib",
     "discretize>=0.11",
     "geoana>=0.7",
-    "empymod>=2.0.0",
+    "libdlf",
 ]
 classifiers = [
     "Development Status :: 4 - Beta",
@@ -74,6 +74,7 @@ docs = [
     "sphinxcontrib-apidoc",
     "pydata-sphinx-theme",
     "nbsphinx",
+    "empymod>=2.0.0",
     "numpydoc",
     "pillow",
     "sympy",
diff --git a/simpeg/electromagnetics/base_1d.py b/simpeg/electromagnetics/base_1d.py
index f1c85f44e..cf73c0ce8 100644
--- a/simpeg/electromagnetics/base_1d.py
+++ b/simpeg/electromagnetics/base_1d.py
@@ -1,8 +1,9 @@
+from collections import namedtuple
+
 from scipy.constants import mu_0
 import numpy as np
 from scipy import sparse as sp
 from scipy.special import roots_legendre
-from empymod.transform import get_dlf_points
 
 from ..simulation import BaseSimulation
 
@@ -17,10 +18,16 @@ from ..utils import (
     validate_integer,
 )
 from .. import props
-from empymod.utils import check_hankel
+import libdlf
 
 __all__ = ["BaseEM1DSimulation"]
 
+HANKEL_FILTERS = {}
+for filter_name in libdlf.hankel.__all__:
+    hankel_filter = getattr(libdlf.hankel, filter_name)
+    if "j0" in hankel_filter.values and "j1" in hankel_filter.values:
+        HANKEL_FILTERS[filter_name] = hankel_filter
+
 ###############################################################################
 #                                                                             #
 #                             Base EM1D Simulation                            #
@@ -137,30 +144,33 @@ class BaseEM1DSimulation(BaseSimulation):
 
         self.hankel_filter = hankel_filter
         self.fix_Jmatrix = fix_Jmatrix
-        # Check input arguments. If self.hankel_filter is not a valid filter,
-        # it will set it to the default (key_201_2009).
-        ht, htarg = check_hankel(
-            "dlf", {"dlf": self.hankel_filter, "pts_per_dec": 0}, 1
-        )
-
-        self._fhtfilt = htarg["dlf"]  # Store filter
-        # self.hankel_pts_per_dec = htarg["pts_per_dec"]  # Store pts_per_dec
         if self.verbose:
             print(">> Use " + self.hankel_filter + " filter for Hankel Transform")
 
     @property
     def hankel_filter(self):
-        """The hankely filter to use.
+        """The hankel filter used.
 
         Returns
         -------
         str
+
+        See Also
+        --------
+        libdlf.hankel
+            The package housing the filter values.
         """
         return self._hankel_filter
 
     @hankel_filter.setter
     def hankel_filter(self, value):
-        self._hankel_filter = validate_string("hankel_filter", value)
+        self._hankel_filter = validate_string(
+            "hankel_filter", value, list(HANKEL_FILTERS.keys())
+        )
+        base, j0, j1 = HANKEL_FILTERS[self._hankel_filter]()
+        hank = namedtuple("HankelFilter", "base j0 j1")
+        self._fhtfilt = hank(base, j0, j1)
+        self._coefficients_set = False
 
     _hankel_pts_per_dec = 0  # Default: Standard DLF
 
@@ -430,9 +440,8 @@ class BaseEM1DSimulation(BaseSimulation):
                     offsets = src.radius * np.ones(rx.locations.shape[0])
 
                 # computations for hankel transform...
-                lambd, _ = get_dlf_points(
-                    self._fhtfilt, offsets, self._hankel_pts_per_dec
-                )
+                lambd = self._fhtfilt.base / offsets[:, None]
+
                 # calculate the source-rx coefficients for the hankel transform
                 C0 = 0.0
                 C1 = 0.0
diff --git a/simpeg/electromagnetics/static/resistivity/simulation_1d.py b/simpeg/electromagnetics/static/resistivity/simulation_1d.py
index 7f913ea3e..b57ea3591 100644
--- a/simpeg/electromagnetics/static/resistivity/simulation_1d.py
+++ b/simpeg/electromagnetics/static/resistivity/simulation_1d.py
@@ -1,28 +1,22 @@
+from collections import namedtuple
+
+import libdlf
 import numpy as np
 
+from ...utils.em1d_utils import get_splined_dlf_points
 from ....simulation import BaseSimulation
 from .... import props
 
 from .survey import Survey
 
-from empymod.transform import get_dlf_points
-from empymod import filters
 from ....utils import validate_type, validate_string
 from scipy.interpolate import InterpolatedUnivariateSpline as iuSpline
 
-
-HANKEL_FILTERS = [
-    "kong_61_2007",
-    "kong_241_2007",
-    "key_101_2009",
-    "key_201_2009",
-    "key_401_2009",
-    "anderson_801_1982",
-    "key_51_2012",
-    "key_101_2012",
-    "key_201_2012",
-    "wer_201_2018",
-]
+HANKEL_FILTERS = {}
+for filter_name in libdlf.hankel.__all__:
+    hankel_filter = getattr(libdlf.hankel, filter_name)
+    if "j0" in hankel_filter.values:
+        HANKEL_FILTERS[filter_name] = hankel_filter
 
 
 def _phi_tilde(rho, thicknesses, lambdas):
@@ -200,11 +194,12 @@ class Simulation1DLayers(BaseSimulation):
     @hankel_filter.setter
     def hankel_filter(self, value):
         self._hankel_filter = validate_string(
-            "hankel_filter",
-            value,
-            HANKEL_FILTERS,
+            "hankel_filter", value, list(HANKEL_FILTERS.keys())
         )
-        self._fhtfilt = getattr(filters, self._hankel_filter)()
+        filt = HANKEL_FILTERS[self._hankel_filter]()
+        hank = namedtuple("HankelFilter", "base j0")
+        self._fhtfilt = hank(filt[0], filt[1])
+        self._coefficients_set = False
 
     @property
     def fix_Jmatrix(self):
@@ -241,9 +236,8 @@ class Simulation1DLayers(BaseSimulation):
                     r_max = max(off.max(), r_max)
         self.survey.set_geometric_factor()
 
-        lambdas, r_spline_points = get_dlf_points(
-            self._fhtfilt, np.r_[r_min, r_max], -1
-        )
+        lambdas, r_spline_points = get_splined_dlf_points(self._fhtfilt, r_min, r_max)
+
         lambdas = lambdas.reshape(-1)
         n_lambda = len(lambdas)
         n_r = len(r_spline_points)
diff --git a/simpeg/electromagnetics/time_domain/simulation_1d.py b/simpeg/electromagnetics/time_domain/simulation_1d.py
index 8611c6713..e2f5b8a6b 100644
--- a/simpeg/electromagnetics/time_domain/simulation_1d.py
+++ b/simpeg/electromagnetics/time_domain/simulation_1d.py
@@ -1,3 +1,5 @@
+from collections import namedtuple
+
 from ..base_1d import BaseEM1DSimulation
 from .sources import StepOffWaveform
 from .receivers import (
@@ -12,13 +14,19 @@ from scipy.constants import mu_0
 from scipy.interpolate import InterpolatedUnivariateSpline as iuSpline
 from scipy.special import roots_legendre
 
-from empymod import filters
-from empymod.transform import get_dlf_points
+import libdlf
 
 from geoana.kernels.tranverse_electric_reflections import rTE_forward, rTE_gradient
 
+from ..utils.em1d_utils import get_splined_dlf_points
 from ...utils import validate_type, validate_string
 
+COS_FILTERS = {}
+for filter_name in libdlf.fourier.__all__:
+    fourier_filter = getattr(libdlf.fourier, filter_name)
+    if "cos" in fourier_filter.values:
+        COS_FILTERS[filter_name] = fourier_filter
+
 
 class Simulation1DLayered(BaseEM1DSimulation):
     """
@@ -26,7 +34,7 @@ class Simulation1DLayered(BaseEM1DSimulation):
     for a single sounding.
     """
 
-    def __init__(self, survey=None, time_filter="key_81_CosSin_2009", **kwargs):
+    def __init__(self, survey=None, time_filter="key_81_2009", **kwargs):
         super().__init__(survey=survey, **kwargs)
         self._coefficients_set = False
         self.time_filter = time_filter
@@ -54,18 +62,20 @@ class Simulation1DLayered(BaseEM1DSimulation):
 
     @time_filter.setter
     def time_filter(self, value):
+        # translate old accepted keys to the names in libdlf for compatibility.
+        if value in [
+            "key_81_CosSin_2009",
+            "key_201_CosSin_2012",
+            "key_601_CosSin_2009",
+        ]:
+            value = value.replace("CosSin_", "")
         self._time_filter = validate_string(
-            "time_filter",
-            value,
-            ["key_81_CosSin_2009", "key_201_CosSin_2012", "key_601_CosSin_2009"],
+            "time_filter", value, list(COS_FILTERS.keys())
         )
-
-        if self._time_filter == "key_81_CosSin_2009":
-            self._fftfilt = filters.key_81_CosSin_2009()
-        elif self._time_filter == "key_201_CosSin_2012":
-            self._fftfilt = filters.key_201_CosSin_2012()
-        elif self._time_filter == "key_601_CosSin_2009":
-            self._fftfilt = filters.key_601_CosSin_2009()
+        filt = COS_FILTERS[self._time_filter]()
+        cos_filt = namedtuple("CosineFilter", "base cos")
+        self._fftfilt = cos_filt(filt[0], filt[-1])
+        self._coefficients_set = False
 
     def get_coefficients(self):
         if self._coefficients_set is False:
@@ -123,8 +133,8 @@ class Simulation1DLayered(BaseEM1DSimulation):
                             f"Unsupported source waveform object of {src.waveform}"
                         )
 
-        omegas, t_spline_points = get_dlf_points(self._fftfilt, np.r_[t_min, t_max], -1)
-        omegas = omegas.reshape(-1)
+        omegas, t_spline_points = get_splined_dlf_points(self._fftfilt, t_min, t_max)
+
         n_omega = len(omegas)
         n_t = len(t_spline_points)
 
diff --git a/simpeg/electromagnetics/utils/em1d_utils.py b/simpeg/electromagnetics/utils/em1d_utils.py
index 91b0ee3f6..c551cb311 100644
--- a/simpeg/electromagnetics/utils/em1d_utils.py
+++ b/simpeg/electromagnetics/utils/em1d_utils.py
@@ -219,3 +219,39 @@ def LogUniform(f, chi_inf=0.05, del_chi=0.05, tau1=1e-5, tau2=1e-2):
     return chi_inf + del_chi * (
         1 - np.log((1 + 1j * w * tau2) / (1 + 1j * w * tau1)) / np.log(tau2 / tau1)
     )
+
+
+def get_splined_dlf_points(filt, v_min, v_max):
+    """
+
+    Parameters
+    ----------
+    filt : namedtuple of numpy.ndarray
+        The filter parameters loaded from libdlf, in a named tuple containing a `.base` property.
+    v_min, v_max : float
+        The minimum and maximum points needed
+
+    Returns
+    -------
+    lamb, spline_points : numpy.ndarray
+        The points to evaluate at in the filter space, and the
+        corresponding spline points in the transformed space.
+    """
+
+    # the below is adapted from empymod.transform.get_dlf_points
+    outmax = filt.base[-1] / v_min
+    outmin = filt.base[0] / v_max
+
+    factor = np.around([filt.base[1] / filt.base[0]], 15)
+
+    pts_per_dec = np.squeeze(1 / np.log(factor))
+
+    nout = int(np.ceil(np.log(outmax / outmin) * pts_per_dec) + 1)
+    if nout - filt.base.size < 3:
+        nout = filt.base.size + 3
+    out = np.exp(
+        np.arange(np.log(outmin), np.log(outmin) + nout / pts_per_dec, 1 / pts_per_dec)
+    )
+
+    spline_points = v_max * np.exp(-np.arange(nout - filt.base.size + 1) / pts_per_dec)
+    return out, spline_points
diff --git a/simpeg/utils/code_utils.py b/simpeg/utils/code_utils.py
index 116215090..5a08cc103 100644
--- a/simpeg/utils/code_utils.py
+++ b/simpeg/utils/code_utils.py
@@ -478,8 +478,8 @@ class Report(ScoobyReport):
             "numpy",
             "scipy",
             "matplotlib",
-            "empymod",
             "geoana",
+            "libdlf",
         ]
 
         # Optional packages.


Test patch:
diff --git a/.ci/environment_test.yml b/.ci/environment_test.yml
index 3663e03a6..e1ca51a0a 100644
--- a/.ci/environment_test.yml
+++ b/.ci/environment_test.yml
@@ -8,7 +8,7 @@ dependencies:
   - matplotlib-base
   - discretize>=0.11
   - geoana>=0.7
-  - empymod>=2.0.0
+  - libdlf
 
 # optional dependencies
   - dask
@@ -28,6 +28,7 @@ dependencies:
   - nbsphinx
   - numpydoc
   - pillow
+  - empymod>=2.0.0
   - sympy
   - memory_profiler
   - python-kaleido
diff --git a/tests/em/em1d/test_EM1D_TD_general_fwd.py b/tests/em/em1d/test_EM1D_TD_general_fwd.py
index 1d4daf5c6..962b4e8b4 100644
--- a/tests/em/em1d/test_EM1D_TD_general_fwd.py
+++ b/tests/em/em1d/test_EM1D_TD_general_fwd.py
@@ -276,5 +276,14 @@ class EM1D_TD_LineCurrent_FwdProblemTests(unittest.TestCase):
         np.testing.assert_allclose(self.bzdt, empymod_solution, rtol=1e-2)
 
 
+def test_backwards_compatible_filter_key():
+
+    srv = tdem.Survey([])
+    sim = tdem.Simulation1DLayered(survey=srv)
+    sim.time_filter = "key_81_CosSin_2009"
+
+    assert sim.time_filter == "key_81_2009"
+
+
 if __name__ == "__main__":
     unittest.main()
diff --git a/tests/em/em1d/test_utils.py b/tests/em/em1d/test_utils.py
new file mode 100644
index 000000000..bbeb4b217
--- /dev/null
+++ b/tests/em/em1d/test_utils.py
@@ -0,0 +1,37 @@
+from collections import namedtuple
+
+import pytest
+import libdlf
+import numpy as np
+import numpy.testing as npt
+from empymod.transform import get_dlf_points
+
+from simpeg.electromagnetics.utils.em1d_utils import get_splined_dlf_points
+
+FILTERS = [f"hankel.{filt}" for filt in libdlf.hankel.__all__] + [
+    f"fourier.{filt}" for filt in libdlf.fourier.__all__
+]
+
+
+@pytest.mark.parametrize("filt", FILTERS)
+@pytest.mark.parametrize("n_points", [1, 5, 10])
+def test_splined_dlf(filt, n_points):
+    f_type, f_name = filt.split(".")
+    f_module = getattr(libdlf, f_type)
+    filt = getattr(f_module, f_name)
+    base, *vals = filt()
+    if len(vals) == 2:
+        v0, v1 = vals
+    else:
+        v0 = v1 = vals[0]
+    factor = np.around([base[1] / base[0]], 15)
+    filt_type = namedtuple("Filter", "base v0 v1 factor")
+    filt = filt_type(base, v0, v1, factor)
+
+    r_s = np.logspace(-1, 1, n_points)
+
+    out1, out2 = get_splined_dlf_points(filt, r_s.min(), r_s.max())
+    test1, test2 = get_dlf_points(filt, r_s, -1)
+
+    npt.assert_allclose(out1, test1[0])
+    npt.assert_allclose(out2, test2)
diff --git a/tests/utils/test_report.py b/tests/utils/test_report.py
index bd6238712..2f38f361c 100644
--- a/tests/utils/test_report.py
+++ b/tests/utils/test_report.py
@@ -16,8 +16,8 @@ class TestReport(unittest.TestCase):
                 "numpy",
                 "scipy",
                 "matplotlib",
-                "empymod",
                 "geoana",
+                "libdlf",
             ],
             # Optional packages.
             optional=[