MAINT: np.trapz -> scipy trapezoid (#974)

* MAINT: np.trapz -> np.trapezoid

* REFACTOR: switch to scipy trapezoid

Author: ColmTalbot

bilby/bilby_mcmc/sampler.py

@@ -6,6 +6,7 @@
 
 import numpy as np
 import pandas as pd
+from scipy.integrate import trapezoid
 from scipy.optimize import differential_evolution
 
 from ..core.result import rejection_sample
@@ -1090,8 +1091,8 @@ def _calculate_stepping_stone(betas, ln_likes):
 
     @staticmethod
     def _compute_evidence_from_mean_lnlikes(betas, mean_lnlikes):
-        lnZ = np.trapz(mean_lnlikes, betas)
-        z2 = np.trapz(mean_lnlikes[::-1][::2][::-1], betas[::-1][::2][::-1])
+        lnZ = trapezoid(mean_lnlikes, betas)
+        z2 = trapezoid(mean_lnlikes[::-1][::2][::-1], betas[::-1][::2][::-1])
         lnZerr = np.abs(lnZ - z2)
         return lnZ, lnZerr
 

bilby/core/prior/interpolated.py

@@ -1,4 +1,5 @@
 import numpy as np
+from scipy.integrate import trapezoid
 from scipy.interpolate import interp1d
 
 from .base import Prior
@@ -163,9 +164,9 @@ def _update_instance(self):
 
     def _initialize_attributes(self):
         from scipy.integrate import cumulative_trapezoid
-        if np.trapz(self._yy, self.xx) != 1:
+        if trapezoid(self._yy, self.xx) != 1:
             logger.debug('Supplied PDF for {} is not normalised, normalising.'.format(self.name))
-        self._yy /= np.trapz(self._yy, self.xx)
+        self._yy /= trapezoid(self._yy, self.xx)
         self.YY = cumulative_trapezoid(self._yy, self.xx, initial=0)
         # Need last element of cumulative distribution to be exactly one.
         self.YY[-1] = 1

bilby/core/sampler/ptemcee.py

@@ -7,6 +7,7 @@
 
 import numpy as np
 import pandas as pd
+from scipy.integrate import trapezoid
 
 from ..utils import check_directory_exists_and_if_not_mkdir, logger, safe_file_dump
 from .base_sampler import (
@@ -1396,9 +1397,9 @@ def compute_evidence(
         mean_lnlikes = mean_lnlikes[~idxs]
         betas = betas[~idxs]
 
-    lnZ = np.trapz(mean_lnlikes, betas)
-    z1 = np.trapz(mean_lnlikes, betas)
-    z2 = np.trapz(mean_lnlikes[::-1][::2][::-1], betas[::-1][::2][::-1])
+    lnZ = trapezoid(mean_lnlikes, betas)
+    z1 = trapezoid(mean_lnlikes, betas)
+    z2 = trapezoid(mean_lnlikes[::-1][::2][::-1], betas[::-1][::2][::-1])
     lnZerr = np.abs(z1 - z2)
 
     if make_plots:
@@ -1410,7 +1411,7 @@ def compute_evidence(
         evidence = []
         for i in range(int(len(betas) / 2.0)):
             min_betas.append(betas[i])
-            evidence.append(np.trapz(mean_lnlikes[i:], betas[i:]))
+            evidence.append(trapezoid(mean_lnlikes[i:], betas[i:]))
 
         ax2.semilogx(min_betas, evidence, "-o")
         ax2.set_ylabel(

bilby/gw/prior.py

@@ -2,7 +2,7 @@
 import copy
 
 import numpy as np
-from scipy.integrate import quad
+from scipy.integrate import quad, trapezoid
 from scipy.interpolate import InterpolatedUnivariateSpline, interp1d
 from scipy.special import hyp2f1
 from scipy.stats import norm
@@ -1520,7 +1520,7 @@ def _build_attributes(self):
         """
         from scipy.integrate import cumulative_trapezoid
         yy = self._all_interped(self.pix_xx)
-        yy /= np.trapz(yy, self.pix_xx)
+        yy /= trapezoid(yy, self.pix_xx)
         YY = cumulative_trapezoid(yy, self.pix_xx, initial=0)
         YY[-1] = 1
         self.inverse_cdf = interp1d(x=YY, y=self.pix_xx, bounds_error=True)

examples/core_examples/multivariate_gaussian_prior.py

@@ -9,6 +9,7 @@
 import numpy as np
 from bilby.core.likelihood import GaussianLikelihood
 from scipy import linalg, stats
+from scipy.integrate import trapezoid
 
 # A few simple setup steps
 label = "multivariate_gaussian_prior"
@@ -75,7 +76,7 @@ def model(time, m, c):
     gp = np.zeros(len(x))
     for i in range(nmodes):  # loop over modes
         gp += weights[i] * stats.norm.pdf(x, loc=mus[i][j], scale=mvg.sigmas[i][j])
-    gp = gp / np.trapz(gp, x)  # renormalise
+    gp = gp / trapezoid(gp, x)  # renormalise
 
     axs[aidx[j]].plot(x, gp, "k--", lw=2)
 

test/core/prior/prior_test.py

@@ -3,6 +3,7 @@
 import numpy as np
 import os
 import scipy.stats as ss
+from scipy.integrate import trapezoid
 
 
 class TestPriorClasses(unittest.TestCase):
@@ -573,7 +574,7 @@ def test_normalized(self):
                 domain = np.linspace(prior.minimum, prior.maximum, 10000)
             else:
                 domain = np.linspace(prior.minimum, prior.maximum, 1000)
-            self.assertAlmostEqual(np.trapz(prior.prob(domain), domain), 1, 3)
+            self.assertAlmostEqual(trapezoid(prior.prob(domain), domain), 1, 3)
 
     def test_accuracy(self):
         """Test that each of the priors' functions is calculated accurately, as compared to scipy's calculations"""

test/gw/likelihood/marginalization_test.py

@@ -9,6 +9,7 @@
 import numpy as np
 import bilby
 from bilby.gw.detector import calibration
+from scipy.integrate import trapezoid
 from scipy.special import logsumexp
 
 
@@ -380,7 +381,7 @@ def _template(self, marginalized, non_marginalized, key, prior=None, values=None
             ln_likes[ii] = non_marginalized.log_likelihood_ratio()
         like = np.exp(ln_likes - max(ln_likes))
 
-        marg_like = np.log(np.trapz(like * prior_values, values)) + max(ln_likes)
+        marg_like = np.log(trapezoid(like * prior_values, values)) + max(ln_likes)
         self.assertAlmostEqual(
             marg_like, marginalized.log_likelihood_ratio(), delta=0.5
         )