fix: roll back default change of kbins

src/powerbox/tools.py

@@ -14,21 +14,47 @@
 from . import dft
 
 
-def _getbins(bins, coords, log):
-    if not np.iterable(bins):
+def _getbins(
+    bins: np.ndarray | int,
+    coord_mags: np.ndarray,
+    log: bool,
+    bins_upto_boxlen: bool | None = None,
+):
+    if np.iterable(bins):
+        return bins
+
+    if bins_upto_boxlen is None:
+        warnings.warn(
+            (
+                "In the future, bins will be generated by default up to the smallest "
+                "length over any dimension, instead of the largest magnitude for the box."
+                "Set bins_upto_boxlen to silence this warning."
+            ),
+            stacklevel=2,
+            category=FutureWarning,
+        )
+        bins_upto_boxlen = False
+
+    if bins_upto_boxlen:
         try:
             # Fails if coords is not a cube / inhomogeneous.
-            max_radius = np.min([np.max(coords, axis=i) for i in range(coords.ndim)])
+            max_radius = np.min(
+                [np.max(coord_mags, axis=i) for i in range(coord_mags.ndim)]
+            )
         except ValueError:
-            maxs = [np.max(coords, axis=i) for i in range(coords.ndim)]
+            maxs = [np.max(coord_mags, axis=i) for i in range(coord_mags.ndim)]
             maxs_flat = []
             [maxs_flat.extend(m.ravel()) for m in maxs]
             max_radius = np.min(maxs_flat)
-        if not log:
-            bins = np.linspace(coords.min(), max_radius, bins + 1)
-        else:
-            mn = coords[coords > 0].min()
-            bins = np.logspace(np.log10(mn), np.log10(max_radius), bins + 1)
+    else:
+        max_radius = coord_mags.max()
+
+    if not log:
+        bins = np.linspace(coord_mags.min(), max_radius, bins + 1)
+    else:
+        mn = coord_mags[coord_mags > 0].min()
+        bins = np.logspace(np.log10(mn), np.log10(max_radius), bins + 1)
+
     return bins
 
 
@@ -44,6 +70,7 @@
     interpolation_method=None,
     interp_points_generator=None,
     return_sumweights=False,
+    bins_upto_boxlen: bool | None = None,
 ):
     r"""
     Average a given field within radial bins.
@@ -107,6 +134,12 @@
         (which can be adjusted by supplying a different interp_points_generator
         function with a different angular resolution).
 
+    bins_upto_boxlen : bool, optional
+        If set to True and the bins are determined automatically, calculate bins only
+        up to the maximum k along any dimension. Otherwise, calculate bins up to the
+        maximum magnitude of k (i.e. a factor of sqrt(ndim) higher). Default is False
+        for backwards compatibility.
+
     Returns
     -------
     field_1d : 1D-array
@@ -149,18 +182,25 @@
         raise ValueError("Only linear interpolation is supported.")
     if len(coords) == len(field.shape):
         # coords are a segmented list of dimensional co-ordinates
-        coords_grid = _magnitude_grid(coords)
+        coord_mags = _magnitude_grid(coords)
     elif interpolation_method is not None:
         raise ValueError(
             "Must supply a list of len(field.shape) of 1D coordinate arrays for coords when interpolating!"
         )
     else:
         # coords are the magnitude of the co-ordinates
         # since we are not interpolating, then we can just use the magnitude of the co-ordinates
-        coords_grid = coords
+        coord_mags = coords
+
     if interpolation_method is None:
         indx, bins, sumweights = _get_binweights(
-            coords_grid, weights, bins, average, bin_ave=bin_ave, log_bins=log_bins
+            coord_mags,
+            weights,
+            bins,
+            average,
+            bin_ave=bin_ave,
+            log_bins=log_bins,
+            bins_upto_boxlen=bins_upto_boxlen,
         )
 
         if np.any(sumweights == 0):
@@ -169,7 +209,8 @@
             )
         res = _field_average(indx, field, weights, sumweights)
     else:
-        bins = _getbins(bins, coords_grid, log_bins)
+        bins = _getbins(bins, coord_mags, log_bins, bins_upto_boxlen)
+
         if bin_ave:
             if log_bins:
                 bins = np.exp((np.log(bins[1:]) + np.log(bins[:-1])) / 2)
@@ -207,18 +248,26 @@
         return np.sqrt(np.sum(np.meshgrid(*([X**2 for X in x]), indexing="ij"), axis=0))
 
 
-def _get_binweights(coords, weights, bins, average=True, bin_ave=True, log_bins=False):
+def _get_binweights(
+    coord_mags,
+    weights,
+    bins,
+    average=True,
+    bin_ave=True,
+    log_bins=False,
+    bins_upto_boxlen: bool | None = None,
+):
     # Get a vector of bin edges
-    bins = _getbins(bins, coords, log_bins)
+    bins = _getbins(bins, coord_mags, log_bins, bins_upto_boxlen=bins_upto_boxlen)
 
-    indx = np.digitize(coords.flatten(), bins)
+    indx = np.digitize(coord_mags.flatten(), bins)
 
     if average or bin_ave:
         if not np.isscalar(weights):
-            if coords.shape != weights.shape:
+            if coord_mags.shape != weights.shape:
                 raise ValueError(
                     "coords and weights must have the same shape!",
-                    coords.shape,
+                    coord_mags.shape,
                     weights.shape,
                 )
             sumweights = np.bincount(
@@ -236,7 +285,9 @@
         if bin_ave:
             bins = (
                 np.bincount(
-                    indx, weights=(weights * coords).flatten(), minlength=len(bins) + 1
+                    indx,
+                    weights=(weights * coord_mags).flatten(),
+                    minlength=len(bins) + 1,
                 )[1:-1]
                 / binweight
             )
@@ -557,6 +608,7 @@
     interpolation_method=None,
     interp_points_generator=None,
     return_sumweights=False,
+    bins_upto_boxlen: bool | None = None,
 ):
     """
     Average the first n dimensions of a given field within radial bins.
@@ -621,6 +673,12 @@
         (which can be adjusted by supplying a different interp_points_generator
         function with a different angular resolution).
 
+    bins_upto_boxlen : bool, optional
+        If set to True and the bins are determined automatically, calculate bins only
+        up to the maximum k along any dimension. Otherwise, calculate bins up to the
+        maximum magnitude of k (i.e. a factor of sqrt(ndim) higher). Default is False
+        for backwards compatibility.
+
     Returns
     -------
     field : (m-n+1)-array
@@ -683,30 +741,37 @@
             interpolation_method=interpolation_method,
             interp_points_generator=interp_points_generator,
             return_sumweights=return_sumweights,
+            bins_upto_boxlen=bins_upto_boxlen,
         )
 
     if len(coords) == len(field.shape):
         # coords are a segmented list of dimensional co-ordinates
-        coords_grid = _magnitude_grid([c for i, c in enumerate(coords) if i < n])
+        coord_mags = _magnitude_grid([c for i, c in enumerate(coords) if i < n])
     elif interpolation_method is not None:
         raise ValueError(
             "Must supply a list of len(field.shape) of 1D coordinate arrays for coords when interpolating!"
         )
     else:
         # coords are the magnitude of the co-ordinates
         # since we are not interpolating, then we can just use the magnitude of the co-ordinates
-        coords_grid = coords
+        coord_mags = coords
 
-    coords_grid = _magnitude_grid([c for i, c in enumerate(coords) if i < n])
+    coord_mags = _magnitude_grid([c for i, c in enumerate(coords) if i < n])
     n1 = np.prod(field.shape[:n])
     n2 = np.prod(field.shape[n:])
     if interpolation_method is None:
         indx, bins, sumweights = _get_binweights(
-            coords_grid, weights, bins, average, bin_ave=bin_ave, log_bins=log_bins
+            coord_mags,
+            weights,
+            bins,
+            average,
+            bin_ave=bin_ave,
+            log_bins=log_bins,
+            bins_upto_boxlen=bins_upto_boxlen,
         )
         res = np.zeros((len(sumweights), n2), dtype=field.dtype)
     if interpolation_method is not None:
-        bins = _getbins(bins, coords_grid, log_bins)
+        bins = _getbins(bins, coord_mags, log_bins, bins_upto_boxlen)
         if bin_ave:
             if log_bins:
                 bins = np.exp((np.log(bins[1:]) + np.log(bins[:-1])) / 2)
@@ -977,6 +1042,7 @@
     interpolation_method=None,
     interp_points_generator=None,
     return_sumweights=False,
+    bins_upto_boxlen: bool | None = None,
 ):
     r"""
     Calculate isotropic power spectrum of a field, or cross-power of two similar fields.
@@ -1065,6 +1131,11 @@
         Note that for the linear interpolation case,
         this corresponds to the number of samples averaged over
         (which can be adjusted with the angular_resolution parameter).
+    bins_upto_boxlen : bool, optional
+        If set to True and the bins are determined automatically, calculate bins only
+        up to the maximum k along any dimension. Otherwise, calculate bins up to the
+        maximum magnitude of k (i.e. a factor of sqrt(ndim) higher). Default is False
+        for backwards compatibility.
 
     Returns
     -------
@@ -1183,6 +1254,7 @@
         interpolation_method=interpolation_method,
         interp_points_generator=interp_points_generator,
         return_sumweights=return_sumweights,
+        bins_upto_boxlen=bins_upto_boxlen,
     )
     res = list(res)
     # Remove shot-noise

tests/test_discrete.py

@@ -1,9 +1,12 @@
 import pytest
 
 import numpy as np
+from functools import partial
 
 from powerbox import LogNormalPowerBox, PowerBox, get_power
 
+get_power = partial(get_power, bins_upto_boxlen=True)
+
 
 def test_discrete_power_gaussian():
     pb = PowerBox(

tests/test_lognormal.py

@@ -1,7 +1,10 @@
 import numpy as np
+from functools import partial
 
 from powerbox import LogNormalPowerBox, PowerBox, get_power
 
+get_power = partial(get_power, bins_upto_boxlen=True)
+
 
 def test_ln_vs_straight():
     # Set up two boxes with exactly the same parameters

tests/test_power.py

@@ -1,8 +1,11 @@
 import numpy as np
 import warnings
+from functools import partial
 
 from powerbox import PowerBox, get_power, ignore_zero_absk, ignore_zero_ki, power2delta
 
+get_power = partial(get_power, bins_upto_boxlen=True)
+
 
 def test_power1d():
     p = [0] * 40