maint: add warnings tests ci

.github/workflows/warnings-tests.yaml

@@ -0,0 +1,53 @@
+name: Warnings Tests
+
+# Test on all pushes, except when the push is literally just a tag (because we
+# tag automatically via CI, and therefore there's no extra code in that push).
+# Also, only test on pull requests into master/dev.
+on:
+  push:
+    tags-ignore:
+      - 'v*'
+  pull_request:
+    branches:
+      - 'master'
+      - 'dev'
+
+jobs:
+  tests:
+    if: "!contains(github.event.pull_request.labels.*.name, 'auto-pr')"
+    env:
+      ENV_NAME: tests
+      PYTHON: ${{ matrix.python-version }}
+      OS: ${{ matrix.os }}
+    name: Testing
+    #    needs: [linter]
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest]
+        python-version: ["3.11"]
+    steps:
+      - uses: actions/checkout@master
+        with:
+          fetch-depth: 1
+      - uses: actions/setup-python@v4
+        with:
+          python-version: ${{ matrix.python-version }}
+      - uses: fortran-lang/setup-fortran@v1
+        id: setup-fortran
+        with:
+          compiler: gcc
+          version: 13
+      - name: Install
+        run:  pip install .[dev]
+      - name: Run Tests
+        run: |
+          python -m pytest -Werror --cov=hmf --cov-config=.coveragerc --cov-report xml:./coverage.xml --durations=25
+      - uses: codecov/codecov-action@v3
+        if: matrix.os == 'ubuntu-latest' && success() && !contains(github.event.pull_request.labels.*.name, 'auto-pr')
+        with:
+          file: ./coverage.xml # optional
+          fail_ci_if_error: true
+          verbose: true
+          token: ${{ secrets.CODECOV_TOKEN }}

development/growth_factor_integration.ipynb

@@ -93,10 +93,10 @@
     "    integrand = 1.0 / (np.exp(lna) * cosmo.efunc(_zvec)) ** 3\n",
     "\n",
     "    if not getvec:\n",
-    "        integral = intg.simpson(np.exp(lna) * integrand, lna,even=\"avg\")\n",
+    "        integral = intg.simpson(np.exp(lna) * integrand, x=lna)\n",
     "        dplus = 5.0 * cosmo.Om0 * cosmo.efunc(z) * integral / 2.0\n",
     "    else:\n",
-    "        integral = intg.cumtrapz(np.exp(lna) * integrand, dx=dlna, initial=0.0)\n",
+    "        integral = intg.cumulative_trapezoid(np.exp(lna) * integrand, dx=dlna, initial=0.0)\n",
     "        dplus = 5.0 * cosmo.Om0 * cosmo.efunc(_zvec) * integral / 2.0\n",
     "\n",
     "    return dplus"
@@ -245,10 +245,10 @@
     "    integrand = 1.0 / (np.exp(lna) * cosmo.efunc(_zvec)) ** 3\n",
     "\n",
     "    if not getvec:\n",
-    "        integral = intg.simpson(integrand, np.exp(lna),even=\"avg\")\n",
+    "        integral = intg.simpson(integrand, x=np.exp(lna))\n",
     "        dplus = 5.0 * cosmo.Om0 * cosmo.efunc(z) * integral / 2.0\n",
     "    else:\n",
-    "        integral = intg.cumtrapz(np.exp(lna) * integrand, dx=dlna, initial=0.0)\n",
+    "        integral = intg.cumulative_trapezoid(np.exp(lna) * integrand, dx=dlna, initial=0.0)\n",
     "        dplus = 5.0 * cosmo.Om0 * cosmo.efunc(_zvec) * integral / 2.0\n",
     "\n",
     "    return dplus"

src/hmf/density_field/halofit.py

@@ -48,7 +48,7 @@ def _get_spec(
     def get_log_sigma2(lnr):
         R = np.exp(lnr)
         integrand = delta_k * np.exp(-((k * R) ** 2))
-        return np.log(_simps(integrand, np.log(k)))
+        return np.log(_simps(integrand, x=np.log(k)))
 
     def get_sigma_abs(lnr):
         return np.abs(get_log_sigma2(lnr))
@@ -77,7 +77,7 @@ def get_sigma_abs(lnr):
     return knl, n_eff, n_curv
 
 
-def halofit(k, delta_k, sigma_8=None, z=0, cosmo=None, takahashi=True):
+def halofit(k, delta_k, *, sigma_8=None, z=0, cosmo=None, takahashi=True):
     """
     Implementation of HALOFIT (Smith+2003).
 

src/hmf/density_field/transfer.py

@@ -323,5 +323,5 @@ def nonlinear_delta_k(self):
         .. math:: \Delta_k = \frac{k^3 P_{\rm nl}(k)}{2\pi^2}
         """
         return _hfit(
-            self.k, self.delta_k, self.sigma_8, self.z, self.cosmo, self.takahashi
+            self.k, self.delta_k, z=self.z, cosmo=self.cosmo, takahashi=self.takahashi
         )

src/hmf/halos/mass_definitions.py

@@ -180,6 +180,8 @@ def change_definition(
 
         if not hasattr(rhos, "__len__"):
             rhos = [rhos]
+
+        if not hasattr(c, "__len__"):
             c = [c]
 
         c_new = np.array(
@@ -356,7 +358,7 @@ def _find_new_concentration(rho_s, halo_density, h=None, x_guess=5.0):
     # provide lower and upper limits for the root finder. To balance stability and
     # performance, we do so iteratively: if there is no result within relatively
     # aggressive limits, we try again with more conservative limits.
-    args = rho_s, halo_density
+    target = halo_density / rho_s
     x = None
     i = 0
     XDELTA_GUESS_FACTORS = [5.0, 10.0, 20.0, 100.0, 10000.0]
@@ -366,22 +368,23 @@ def _find_new_concentration(rho_s, halo_density, h=None, x_guess=5.0):
         def h(x):
             return np.log(1.0 + x) - x / (1.0 + x)
 
-    fnc = (
-        lambda x, rhos, density_threshold: rhos * h(x) * 3.0 / x**3 - density_threshold
-    )
+    def fnc(x):
+        return h(x) * 3.0 / x**3 - target
 
     while x is None and i < len(XDELTA_GUESS_FACTORS):
         try:
             xmin = x_guess / XDELTA_GUESS_FACTORS[i]
             xmax = x_guess * XDELTA_GUESS_FACTORS[i]
-            x = sp.optimize.brentq(fnc, xmin, xmax, args)
-        except Exception:
+            x = sp.optimize.brentq(fnc, xmin, xmax)
+        except Exception as e:
+            warnings.warn(f"raised following error: {e}")
             i += 1
 
     if x is None:
         raise OptimizationException(
-            "Could not determine x where the density threshold %.2f is satisfied."
-            % halo_density
+            "Could not determine x where the density threshold "
+            f"{halo_density / rho_s:.2f} is satisfied. Largest x-range tried was "
+            f"x={xmin} -- {xmax} which brackets {fnc(xmin)} -- {fnc(xmax)}"
         )
 
     return x

src/hmf/helpers/sample.py

@@ -13,20 +13,24 @@
 def _prepare_mf(log_mmin, **mf_kwargs):
     h = hmf.MassFunction(Mmin=log_mmin, **mf_kwargs)
     mask = h.ngtm > 0
-    icdf = _spline((h.ngtm[mask] / h.ngtm[0])[::-1], np.log10(h.m[mask][::-1]), k=3)
+    _icdf = np.log10(h.ngtm[mask] / h.ngtm[0])
+    icdf = _spline(_icdf[::-1], np.log10(h.m[mask][::-1]), k=3)
 
     return icdf, h
 
 
-def _choose_halo_masses_num(N, icdf, xmin=0):
+def _choose_halo_masses_num(N, icdf, xmin=0, rng=None):
     # Generate random variates from 0 to maxcum
-    x = np.random.uniform(low=xmin, high=1.0, size=int(N))
+    if rng is None:
+        rng = np.random.default_rng()
 
+    x = rng.uniform(low=xmin, high=1.0, size=int(N))
+    logm = icdf(np.log10(x))
     # Generate halo masses from mf distribution
-    return 10 ** icdf(x)
+    return 10**logm
 
 
-def sample_mf(N, log_mmin, sort=False, **mf_kwargs):
+def sample_mf(N, log_mmin, sort=False, rng=None, **mf_kwargs):
     """
     Create a sample of halo masses from a theoretical mass function.
 
@@ -60,8 +64,7 @@ def sample_mf(N, log_mmin, sort=False, **mf_kwargs):
     >>> m,hmf = sample_mf(1e6,10.0,hmf_model="PS",Mmax=17)
     """
     icdf, h = _prepare_mf(log_mmin, **mf_kwargs)
-
-    m = _choose_halo_masses_num(N, icdf, xmin=h.ngtm.min() / h.ngtm[0])
+    m = _choose_halo_masses_num(N, icdf, xmin=h.ngtm.min() / h.ngtm[0], rng=rng)
 
     if sort:
         m.sort()

src/hmf/mass_function/integrate_hmf.py

@@ -13,7 +13,7 @@ class NaNException(Exception):
     pass
 
 
-def hmf_integral_gtm(M, dndm, mass_density=False):
+def hmf_integral_gtm(m, dndm, mass_density=False):
     """
     Cumulatively integrate dn/dm.
 
@@ -52,15 +52,17 @@ def hmf_integral_gtm(M, dndm, mass_density=False):
     True
 
     """
+    n = len(m)
+
     # Eliminate NaN's
-    m = M[np.logical_not(np.isnan(dndm))]
-    dndm = dndm[np.logical_not(np.isnan(dndm))]
+    mask = np.isfinite(dndm)
+    m = m[mask]
+    dndm = dndm[mask]
     dndlnm = m * dndm
 
     if len(m) < 4:
         raise NaNException(
-            "There are too few real numbers in dndm: len(dndm) = %s, #NaN's = %s"
-            % (len(M), len(M) - len(dndm))
+            f"There are too few real numbers in dndm: len(dndm) = {n}, #NaN's = {n - len(m)}"
         )
 
     # Calculate the mass function (and its integral) from the highest M up to 10**18
@@ -70,30 +72,26 @@ def hmf_integral_gtm(M, dndm, mass_density=False):
         mf = mf_func(m_upper)
 
         if not mass_density:
-            int_upper = intg.simpson(
-                np.exp(mf), dx=m_upper[2] - m_upper[1], even="first"
-            )
+            int_upper = intg.simpson(np.exp(mf), dx=m_upper[2] - m_upper[1])
         else:
-            int_upper = intg.simpson(
-                np.exp(m_upper + mf), dx=m_upper[2] - m_upper[1], even="first"
-            )
+            int_upper = intg.simpson(np.exp(m_upper + mf), dx=m_upper[2] - m_upper[1])
     else:
         int_upper = 0
 
     # Calculate the cumulative integral (backwards) of [m*]dndlnm
     if not mass_density:
         ngtm = np.concatenate(
             (
-                intg.cumtrapz(dndlnm[::-1], dx=np.log(m[1]) - np.log(m[0]))[::-1],
+                intg.cumulative_trapezoid(dndlnm[::-1], dx=np.log(m[1] / m[0]))[::-1],
                 np.zeros(1),
             )
         )
     else:
         ngtm = np.concatenate(
             (
-                intg.cumtrapz(m[::-1] * dndlnm[::-1], dx=np.log(m[1]) - np.log(m[0]))[
-                    ::-1
-                ],
+                intg.cumulative_trapezoid(
+                    m[::-1] * dndlnm[::-1], dx=np.log(m[1]) - np.log(m[0])
+                )[::-1],
                 np.zeros(1),
             )
         )

tests/test_fcoll.py

@@ -80,12 +80,8 @@ def test_ranges_cut(self, peacock, Mmin, Mmax):
 
         anl = fcoll_Peacock(np.sqrt(peacock.nu))
         num = peacock.rho_gtm / peacock.mean_density0
-        err = np.abs((num - anl) / anl)[
-            np.logical_and(peacock.m > 10**10, peacock.m < 10**15)
-        ]
-        err = err[np.logical_not(np.isnan(err))]
-        print(np.max(err))
-        assert np.max(err) < 0.4
+        mask = np.logical_and(peacock.m > 10**10, peacock.m < 10**15)
+        np.testing.assert_allclose(num[mask], anl[mask], rtol=0.4)
 
     @pytest.fixture(scope="class")
     def ps(self):

tests/test_filters.py

@@ -137,12 +137,10 @@ def test_sigma_Rhalf(self, cls):
         t = 2 + 2 + 1
         true = 1.0 / (2 * pi**2 * t * thisr**t)
 
-        with warnings.catch_warnings(record=True) as w:
+        with pytest.warns(UserWarning):
             # should also raise a warning
             R = 0.5
             s2 = cls.sigma(R)[0] ** 2
-            assert w
-        print(s2, true)
         assert np.isclose(s2, true)
 
     def test_sigma1_Rhalf(self, cls):
@@ -151,13 +149,10 @@ def test_sigma1_Rhalf(self, cls):
         t = 4 + 2 + 1
         true = 1.0 / (2 * pi**2 * t * thisr**t)
 
-        with warnings.catch_warnings(record=True) as w:
+        with pytest.warns(UserWarning):
             # should also raise a warning
             R = 0.5
             s2 = cls.sigma(R, 1)[0] ** 2
-            assert w
-
-        print(s2, true)
         assert np.isclose(s2, true)
 
     def test_dlnssdlnr_Rhalf(self, cls):

tests/test_hmf.py

@@ -1,9 +1,9 @@
 """Tests of HMF."""
 
+import pytest
 from pytest import raises
 
 import numpy as np
-import warnings
 
 from hmf import MassFunction
 
@@ -47,6 +47,5 @@ def test_str_filter():
 
 def test_mass_nonlinear_outside_range():
     h = MassFunction(Mmin=8, Mmax=9)
-    with warnings.catch_warnings(record=True) as w:
+    with pytest.warns(UserWarning):
         assert h.mass_nonlinear > 0
-        assert len(w)

tests/test_integrate_hmf.py

@@ -48,9 +48,10 @@ def test_high_z(self):
         m = np.logspace(10, 18, 500)
         dndm = self.tggd(m, 9.0, -1.93, 0.4)
         ngtm = self.anl_int(m, 9.0, -1.93, 0.4)
-
-        print(ngtm / hmf_integral_gtm(m, dndm))
-        assert np.allclose(ngtm, hmf_integral_gtm(m, dndm), rtol=0.03)
+        mask = m < 10**15
+        np.testing.assert_allclose(
+            ngtm[mask], hmf_integral_gtm(m, dndm)[mask], rtol=0.03, atol=1e-8
+        )
 
     # def test_low_mmax_z0(self):
     #     m = np.logspace(10,15,500)

tests/test_sample.py

@@ -1,14 +1,12 @@
-from pytest import raises
-
 import numpy as np
 from scipy.interpolate import InterpolatedUnivariateSpline as spline
 
 from hmf.helpers.sample import dndm_from_sample, sample_mf
 
 
 def test_circular():
-    np.random.seed(1234)
-    m, h = sample_mf(1e5, 11, transfer_model="EH")
+    rng = np.random.default_rng(1234)
+    m, h = sample_mf(1e5, 11, transfer_model="EH", rng=rng)
     centres, hist = dndm_from_sample(m, 1e5 / h.ngtm[0])
 
     s = spline(np.log10(h.m), np.log10(h.dndm))
@@ -18,12 +16,9 @@ def test_circular():
 
 
 def test_mmax_big():
-    # raises ValueError because ngtm=0 exactly at m=18
-    # due to hard limit of integration in integrate_hmf.
-    np.random.seed(12345)
+    rng = np.random.default_rng(12345)
 
-    m, h = sample_mf(1e5, 11, transfer_model="EH", Mmax=18)
+    m, h = sample_mf(1e5, 11, transfer_model="EH", Mmax=18, rng=rng)
     # centres,hist = dndm_from_sample(m,1e5/h.ngtm[0])
     # print centres,hist
-    with raises(ValueError):
-        dndm_from_sample(m, 1e5 / h.ngtm[0])
+    dndm_from_sample(m, 1e5 / h.ngtm[0])