Merge pull request #115 from TeamLEGWORK/new-python

Update LEGWORK to work with the latest python

.github/workflows/pypi_upload.yml

@@ -17,7 +17,7 @@ jobs:
     - uses: actions/setup-python@v2
     - uses: actions/setup-python@v2
       with:
-        python-version: '3.9'
+        python-version: '3.11'
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip

.github/workflows/python-package.yml

@@ -22,13 +22,16 @@ on:
       - 'setup.py'
       - 'pyproject.toml'
 
+permissions:
+  pull-requests: write
+
 jobs:
   build:
 
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7, 3.8]
+        python-version: [3.11]
 
     steps:
     - uses: actions/checkout@v2
@@ -48,21 +51,16 @@ jobs:
         flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics
         # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
-    - name: Test with pytest and generate coverage report
+    - name: Test the code with pytest
       run: |
         # test with as many cores as possible for speed
         pytest legwork -n=auto --cov=./ --cov-report=xml
+    - name: Get Coverage
+      uses: orgoro/[email protected]
+      with:
+        coverageFile: ./coverage.xml
+        token: ${{ secrets.GITHUB_TOKEN }}
+    - name: Test the notebooks with pytest
+      run: |
         # test notebooks in addition to the code
         pytest -n=auto --nbmake "docs/demos"
-    - name: Upload coverage to Codecov
-      uses: codecov/codecov-action@v1
-      with:
-        token: ${{ secrets.CODECOV_TOKEN }}
-        files: ./coverage.xml
-        directory: ./
-        flags: unittests
-        env_vars: OS,PYTHON
-        name: codecov-umbrella
-        fail_ci_if_error: true
-        path_to_write_report: ./coverage/codecov_report.txt
-        verbose: true

changelog.md

@@ -173,3 +173,9 @@ having two when we could just set `f_dom=2 f_orb`
 ## 0.4.7
 *TW 12/10/23*
 - [Issue [#111](https://github.com/TeamLEGWORK/LEGWORK/issues/111)] Allow `plot_sources_on_sc` to change the underlying sensitivity curve directly with `sc_vis_settings`
+
+## 0.5.0
+*TW 22/06/24*
+- Update to Python 3.11 and various dependencies
+    - Fix deprecations based on these updates
+- Add __repr__ and __len__ functions for the Source class and its subclasses

environment.yml

@@ -7,16 +7,16 @@ channels:
   - defaults
   - conda-forge
 dependencies:
-  - python >= 3.7
-  - pip >= 21.0.0
-  - numba >= 0.50
-  - numpy >= 1.17
-  - astropy >= 4.0
-  - scipy >= 1.5.0
-  - matplotlib >= 3.3.2
-  - seaborn >= 0.11.1
-  - schwimmbad >= 0.3.2
+  - python >= 3.11
+  - pip >= 24.0
+  - numba >= 0.58
+  - numpy >= 1.26
+  - astropy >= 6.1
+  - scipy >= 1.13
+  - matplotlib >= 3.8
+  - seaborn >= 0.13.2
+  - schwimmbad >= 0.4.2
   - jupyter
   - ipython
   - pip:
-    - legwork >= 0.4.6
+    - legwork >= 0.5.0

legwork/_version.py

@@ -1 +1 @@
-__version__ = "0.4.7"
+__version__ = "0.5.0"

legwork/psd.py

@@ -6,6 +6,7 @@
 import astropy.constants as const
 from scipy.interpolate import splev, splrep
 from importlib import resources
+import os
 
 __all__ = ['load_response_function', 'approximate_response_function', 'power_spectral_density',
            'lisa_psd', 'tianqin_psd', 'get_confusion_noise', 'get_confusion_noise_robson19',
@@ -34,8 +35,7 @@ def load_response_function(f, fstar=19.09e-3):
     """
     # try to load the values for interpolating R
     try:
-        with resources.path(package="legwork", resource="R.npy") as path:
-            f_R, R = np.load(path)
+        f_R, R = np.load(os.path.join(resources.files("legwork"),"R.npy"), allow_pickle=True)
     except FileExistsError:  # pragma: no cover
         print("WARNING: Can't find response function file, using approximation instead")
         return approximate_response_function(f, fstar)

legwork/snr.py

@@ -36,13 +36,12 @@ def snr_circ_stationary(m_c, f_orb, dist, t_obs, position=None, polarisation=Non
         Inclination of the source. Must have astropy angular units.
 
     interpolated_g : `function`
-        A function returned by :class:`scipy.interpolate.interp2d` that computes g(n,e) from Peters (1964).
-        The code assumes that the function returns the output sorted as with the interp2d returned functions
-        (and thus unsorts). Default is None and uses exact g(n,e) in this case.
+        A function returned by :class:`scipy.interpolate.RectBivariateSpline` that computes g(n,e)
+        from Peters (1964). Default is None and uses exact g(n,e) in this case.
 
     interpolated_sc : `function`
         A function returned by :class:`scipy.interpolate.interp1d` that computes the LISA sensitivity curve.
-        Default is None and uses exact values. Note: take care to ensure that your interpolated function has
+        Default is None and uses exact values. Note: take care to ensure that yourinterpolated function has
         the same LISA observation time as ``t_obs`` and uses the same instrument.
 
     **kwargs : `various`
@@ -96,9 +95,8 @@ def snr_ecc_stationary(m_c, f_orb, ecc, dist, t_obs, harmonics_required,
         Maximum integer harmonic to compute
 
     interpolated_g : `function`
-        A function returned by :class:`scipy.interpolate.interp2d` that computes g(n,e) from Peters (1964).
-        The code assumes that the function returns the output sorted as with the interp2d returned functions
-        (and thus unsorts). Default is None and uses exact g(n,e) in this case.
+        A function returned by :class:`scipy.interpolate.RectBivariateSpline` that computes g(n,e)
+        from Peters (1964). Default is None and uses exact g(n,e) in this case.
 
     interpolated_sc : `function`
         A function returned by :class:`scipy.interpolate.interp1d` that computes the LISA sensitivity curve.
@@ -186,9 +184,8 @@ def snr_circ_evolving(m_1, m_2, f_orb_i, dist, t_obs, n_step, t_merge=None,
         Time until merger
 
     interpolated_g : `function`
-        A function returned by :class:`scipy.interpolate.interp2d` that computes g(n,e) from Peters (1964).
-        The code assumes that the function returns the output sorted as with the interp2d returned functions
-        (and thus unsorts). Default is None and uses exact g(n,e) in this case.
+        A function returned by :class:`scipy.interpolate.RectBivariateSpline` that computes g(n,e)
+        from Peters (1964). Default is None and uses exact g(n,e) in this case.
 
     interpolated_sc : `function`
         A function returned by :class:`scipy.interpolate.interp1d` that computes the LISA sensitivity curve.
@@ -274,9 +271,8 @@ def snr_ecc_evolving(m_1, m_2, f_orb_i, dist, ecc, harmonics_required, t_obs, n_
         Time until merger
 
     interpolated_g : `function`
-        A function returned by :class:`scipy.interpolate.interp2d` that computes g(n,e) from Peters (1964).
-        The code assumes that the function returns the output sorted as with the interp2d returned functions
-        (and thus unsorts). Default is None and uses exact g(n,e) in this case.
+        A function returned by :class:`scipy.interpolate.RectBivariateSpline` that computes g(n,e)
+        from Peters (1964). Default is None and uses exact g(n,e) in this case.
 
     interpolated_sc : `function`
         A function returned by :class:`scipy.interpolate.interp1d` that computes the LISA sensitivity curve.

legwork/source.py

@@ -3,7 +3,8 @@
 from astropy.coordinates import SkyCoord
 import numpy as np
 from importlib import resources
-from scipy.interpolate import interp1d, interp2d
+from scipy.interpolate import interp1d, RectBivariateSpline
+import os
 
 from legwork import utils, strain, psd, evol
 import legwork.snr as sn
@@ -212,6 +213,12 @@ def __init__(self, m_1, m_2, ecc, dist, n_proc=1, f_orb=None, a=None, position=N
         self.set_g(interpolate_g)
         self.set_sc()
 
+    def __repr__(self):
+        return f"<Source: {self.n_sources} sources>"
+
+    def __len__(self):
+        return self.n_sources
+
     def create_harmonics_functions(self):
         """Create two harmonics related functions as methods for the Source class
 
@@ -225,8 +232,7 @@ def create_harmonics_functions(self):
         the maximum strain for a system with eccentricity `ecc`."""
 
         # open file containing pre-calculated g(n,e) and F(e) values
-        with resources.path(package="legwork", resource="harmonics.npz") as path:
-            lum_info = np.load(path)
+        lum_info = np.load(os.path.join(resources.files("legwork"), "harmonics.npz"), allow_pickle=True)
 
         e_min, e_max, e_len = lum_info["e_lims"]
         e_len = e_len.astype(int)
@@ -309,14 +315,13 @@ def set_g(self, interpolate_g):
         """
         if interpolate_g:
             # open file containing pre-calculated fine g(n,e) grid
-            with resources.path(package="legwork",
-                                resource="peters_g.npy") as path:
-                peters_g = np.load(path)
+            peters_g = np.load(os.path.join(resources.files("legwork"), "peters_g.npy"))
 
             # interpolate grid using scipy
             n_range = np.arange(1, 10000 + 1).astype(int)
             e_range = np.linspace(0, 1, 1000)
-            self.g = interp2d(n_range, e_range, peters_g, kind="cubic")
+            f = RectBivariateSpline(n_range, e_range, peters_g.T)
+            self.g = lambda n, e: f.ev(n, e)
         else:
             self.g = None
 
@@ -1282,6 +1287,9 @@ def get_snr(self, t_obs=None, instrument=None, custom_psd=None, verbose=False):
                                            verbose=verbose)
         return self.snr
 
+    def __repr__(self):
+        return f"<Stationary: {self.n_sources} stationary sources>"
+
 
 class Evolving(Source):
     """Subclass for sources that are evolving"""
@@ -1292,15 +1300,18 @@ def get_snr(self, t_obs=None, instrument=None, custom_psd=None, n_step=100, verb
         return self.snr
 
 
+    def __repr__(self):
+        return f"<Evolving: {self.n_sources} evolving sources>"
+
+
 class VerificationBinaries(Source):
     """Generate a Source class with the LISA verification binaries preloaded.
     Data for the binaries is gathered from Kupfer+18 Table 1 and 2."""
 
     def __init__(self):
         # open file containing verification binary data
-        with resources.path(package="legwork", resource="verification_binaries.npy") as path:
-            vbs = np.load(path, allow_pickle=True)
-            vbs = vbs.item()
+        vbs = np.load(os.path.join(resources.files("legwork"),
+                                   "verification_binaries.npy"), allow_pickle=True).item()
 
         position = SkyCoord(l=vbs["l_gal"], b=vbs["b_gal"], distance=vbs["dist"], frame="galactic")
 
@@ -1313,3 +1324,7 @@ def __init__(self):
         self.labels = vbs["label"]
         self.true_snr = np.array(vbs["snr"])
         self.max_snr_harmonic = np.repeat(2, self.n_sources).astype(int)
+
+
+    def __repr__(self):
+        return f"<VerificationBinaries>"

legwork/strain.py

@@ -90,9 +90,8 @@ def h_0_n(m_c, f_orb, ecc, n, dist, position=None, polarisation=None, inclinatio
         Inclination of the source. Must have astropy angular units.
 
     interpolated_g : `function`
-        A function returned by :class:`scipy.interpolate.interp2d` that computes g(n,e) from Peters (1964).
-        The code assumes that the function returns the output sorted as with the interp2d returned functions
-        (and thus unsorts). Default is None and uses exact g(n,e) in this case.
+        A function returned by :class:`scipy.interpolate.RectBivariateSpline` that computes g(n,e)
+        from Peters (1964). Default is None and uses exact g(n,e) in this case.
 
     Returns
     -------
@@ -117,27 +116,20 @@ def h_0_n(m_c, f_orb, ecc, n, dist, position=None, polarisation=None, inclinatio
     prefac = (2**(28/3) / 5)**(0.5) * c.G**(5/3) / c.c**4
     n_independent_part = prefac * m_c**(5/3) * (np.pi * f_orb)**(2/3) / dist
 
+    # extend harmonic and eccentricity dimensions to full (x, y, z)
+    n = n[np.newaxis, np.newaxis, :]
+    ecc = ecc[..., np.newaxis]
+
     # check whether to interpolate g(n, e)
     if interpolated_g is None:
-        # extend harmonic and eccentricity dimensions to full (x, y, z)
-        n = n[np.newaxis, np.newaxis, :]
-        ecc = ecc[..., np.newaxis]
         n_dependent_part = utils.peters_g(n, ecc)**(1/2) / n
     else:
-        # flatten array to work nicely interp2d
-        g_vals = interpolated_g(n, ecc.flatten())
+        g_vals = interpolated_g(n, ecc)
 
         # set negative values from cubic fit to 0.0
         g_vals[g_vals < 0.0] = 0.0
 
-        # unsort the output array if there is more than one eccentricity
-        if isinstance(ecc, (np.ndarray, list)) and len(ecc) > 1:
-            g_vals = g_vals[np.argsort(ecc.flatten()).argsort()]
-
-        # reshape output to proper dimensions
-        g_vals = g_vals.reshape((*ecc.shape, len(n)))
-
-        n_dependent_part = g_vals**(0.5) / n[np.newaxis, np.newaxis, :]
+        n_dependent_part = g_vals**(0.5) / n
 
     h_0 = n_independent_part[..., np.newaxis] * n_dependent_part
 
@@ -186,9 +178,8 @@ def h_c_n(m_c, f_orb, ecc, n, dist, position=None, polarisation=None, inclinatio
         Inclination of the source. Must have astropy angular units.
 
     interpolated_g : `function`
-        A function returned by :class:`scipy.interpolate.interp2d` that computes g(n,e) from Peters (1964).
-        The code assumes that the function returns the output sorted as with the interp2d returned functions
-        (and thus unsorts). Default is None and uses exact g(n,e) in this case.
+        A function returned by :class:`scipy.interpolate.RectBivariateSpline` that computes g(n,e)
+        from Peters (1964). Default is None and uses exact g(n,e) in this case.
 
     Returns
     -------
@@ -213,27 +204,20 @@ def h_c_n(m_c, f_orb, ecc, n, dist, position=None, polarisation=None, inclinatio
     prefac = (2**(5/3) / (3 * np.pi**(4/3)))**(0.5) * c.G**(5/6) / c.c**(3/2)
     n_independent_part = prefac * m_c**(5/6) / dist * f_orb**(-1/6) / utils.peters_f(ecc)**(0.5)
 
+    # extend harmonic and eccentricity dimensions to full (x, y, z)
+    n = n[np.newaxis, np.newaxis, :]
+    ecc = ecc[..., np.newaxis]
+
     # check whether to interpolate g(n, e)
     if interpolated_g is None:
-        # extend harmonic and eccentricity dimensions to full (x, y, z)
-        n = n[np.newaxis, np.newaxis, :]
-        ecc = ecc[..., np.newaxis]
         n_dependent_part = (utils.peters_g(n, ecc) / n)**(1/2)
     else:
-        # flatten array to work nicely interp2d
-        g_vals = interpolated_g(n, ecc.flatten())
+        g_vals = interpolated_g(n, ecc)
 
         # set negative values from cubic fit to 0.0
         g_vals[g_vals < 0.0] = 0.0
 
-        # unsort the output array if there is more than one eccentricity
-        if isinstance(ecc, (np.ndarray, list)) and len(ecc) > 1:
-            g_vals = g_vals[np.argsort(ecc.flatten()).argsort()]
-
-        # reshape output to proper dimensions
-        g_vals = g_vals.reshape((*ecc.shape, len(n)))
-
-        n_dependent_part = (g_vals / n[np.newaxis, np.newaxis, :])**(0.5)
+        n_dependent_part = (g_vals / n)**(0.5)
 
     h_c = n_independent_part[..., np.newaxis] * n_dependent_part
 

legwork/tests/test_evol.py

@@ -162,7 +162,7 @@ def test_mandel_fit(self):
                                           ecc_i=ecc, exact=False).to(u.yr)
         single_time = evol.get_t_merge_ecc(beta=beta[0], a_i=a_i[0],
                                            ecc_i=ecc[0], exact=False).to(u.yr)
-        self.assertTrue(array_time[0] == single_time)
+        self.assertTrue(np.isclose(array_time[0], single_time))
 
     def test_timestep_creation(self):
         n_values = 10

requirements.txt

@@ -2,11 +2,11 @@
 #   - setup.cfg
 #   - environment.yml
 #   - install.rst
-numba >= 0.50
-numpy >= 1.17
-astropy >= 4.0
-scipy >= 1.5.0
-matplotlib >= 3.3.2
-seaborn >= 0.11.1
-schwimmbad >= 0.3.2
+numba >= 0.58
+numpy >= 1.26
+astropy >= 6.1
+scipy >= 1.13
+matplotlib >= 3.8
+seaborn >= 0.13.2
+schwimmbad >= 0.4.2
 Jinja2 == 2.11