From bd9638897ae32375ddcaae6e339bd2ca2ab93324 Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Sat, 3 Jan 2026 17:57:51 -0600 Subject: [PATCH 01/15] np.in1d --> np.isin --- .../subhalo_based_models/subhalo_selection_kernel.py | 2 +- halotools/utils/array_indexing_manipulations.py | 2 +- halotools/utils/crossmatch.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/halotools/empirical_models/phase_space_models/subhalo_based_models/subhalo_selection_kernel.py b/halotools/empirical_models/phase_space_models/subhalo_based_models/subhalo_selection_kernel.py index 08851c98b..5c403388b 100644 --- a/halotools/empirical_models/phase_space_models/subhalo_based_models/subhalo_selection_kernel.py +++ b/halotools/empirical_models/phase_space_models/subhalo_based_models/subhalo_selection_kernel.py @@ -119,7 +119,7 @@ def calculate_satellite_selection_mask(subhalo_hostids, satellite_occupations, h those rare subhalos with no matching host halo (this situation occurs in <0.1% for typical Rockstar catalogs). - >>> matched_mask = np.in1d(halocat.halo_table['halo_hostid'], halocat.halo_table['halo_id']) + >>> matched_mask = np.isin(halocat.halo_table['halo_hostid'], halocat.halo_table['halo_id']) >>> halos = halocat.halo_table[matched_mask] Now we will sort the catalog by the ``sorting_keys`` list. diff --git a/halotools/utils/array_indexing_manipulations.py b/halotools/utils/array_indexing_manipulations.py index a99162af6..7618d4863 100644 --- a/halotools/utils/array_indexing_manipulations.py +++ b/halotools/utils/array_indexing_manipulations.py @@ -337,7 +337,7 @@ def calculate_entry_multiplicity(sorted_repeated_hostids, unique_possible_hostid unique_appearances_of_hostid, unique_entry_multiplicity = ( np.unique(sorted_repeated_hostids, return_counts=True)) - hostid_has_match = np.in1d(unique_possible_hostids, unique_appearances_of_hostid, + hostid_has_match = np.isin(unique_possible_hostids, unique_appearances_of_hostid, assume_unique=True) entry_multiplicity = np.zeros_like(unique_possible_hostids) diff --git a/halotools/utils/crossmatch.py b/halotools/utils/crossmatch.py index 12da33ee1..f8ebf44b6 100644 --- a/halotools/utils/crossmatch.py +++ b/halotools/utils/crossmatch.py @@ -162,7 +162,7 @@ def crossmatch(x, y, skip_bounds_checking=False): unique_xvals, counts = np.unique(x_sorted, return_counts=True) # Determine which of the unique x values has a match in y - unique_xval_has_match = np.in1d(unique_xvals, y_sorted, assume_unique=True) + unique_xval_has_match = np.isin(unique_xvals, y_sorted, assume_unique=True) # Create a boolean array with True for each value in x with a match, otherwise False idx_x = np.repeat(unique_xval_has_match, counts) From bb55072a5012ed5e094794a0a70416f0681de866 Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 12:36:51 -0600 Subject: [PATCH 02/15] Update changelog --- CHANGES.rst | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index c19190a1b..4c31b86dc 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,24 +1,28 @@ +0.9.4 (unreleased) +------------------ + + 0.9.3 (2025-02-10) ------------------ - + - Compatibility with Numpy 2.0 (https://github.com/astropy/halotools/pull/1114) 0.9.2 (2024-12-02) ------------------ - + - Compatibility with Astropy 7 (https://github.com/astropy/halotools/pull/1107) 0.9.1 (2024-09-11) ------------------ - + - Repackage halotools according to pyproject.toml (https://github.com/astropy/halotools/pull/1093) 0.9.0 (2024-09-06) ------------------ - + - Add models and mock observations of intrinsic alignments (https://github.com/astropy/halotools/pull/1056) From 638ed67045d1359e2a1045a7f0ab91f37c479c20 Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 13:45:52 -0600 Subject: [PATCH 03/15] Fix bug in treatment of period in large_scale_density_spherical_volume.py --- .../large_scale_density_spherical_volume.py | 87 ++++--- ...st_large_scale_density_spherical_volume.py | 58 +++-- .../mock_observables_helpers.py | 91 ++++--- .../pair_counters/rectangular_mesh.py | 239 +++++++++++++----- 4 files changed, 311 insertions(+), 164 deletions(-) diff --git a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py index 80860fe79..f569f1bd5 100644 --- a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py +++ b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py @@ -2,25 +2,32 @@ Module containing functions used to determine various ways of quantifying the large-scale density of a set of points. """ + from __future__ import absolute_import, division, print_function, unicode_literals import numpy as np -from ..pair_counters import npairs_per_object_3d - from ...custom_exceptions import HalotoolsError +from ..mock_observables_helpers import get_period +from ..pair_counters import npairs_per_object_3d +__all__ = ("large_scale_density_spherical_volume",) -__all__ = ('large_scale_density_spherical_volume', ) - -__author__ = ('Andrew Hearin', ) +__author__ = ("Andrew Hearin",) -np.seterr(divide='ignore', invalid='ignore') # ignore divide by zero in e.g. DD/RR +np.seterr(divide="ignore", invalid="ignore") # ignore divide by zero in e.g. DD/RR -def large_scale_density_spherical_volume(sample, tracers, radius, - period=None, sample_volume=None, num_threads=1, approx_cell1_size=None, - norm_by_mean_density=False): +def large_scale_density_spherical_volume( + sample, + tracers, + radius, + period=None, + sample_volume=None, + num_threads=1, + approx_cell1_size=None, + norm_by_mean_density=False, +): """ Calculate the mean density of the input ``sample`` from an input set of tracer particles using a sphere centered on each point @@ -100,50 +107,50 @@ def large_scale_density_spherical_volume(sample, tracers, radius, """ sample, tracers, rbins, period, sample_volume, num_threads, approx_cell1_size = ( _large_scale_density_spherical_volume_process_args( - sample, tracers, radius, period, sample_volume, num_threads, approx_cell1_size) + sample, + tracers, + radius, + period, + sample_volume, + num_threads, + approx_cell1_size, ) + ) - result = npairs_per_object_3d(sample, tracers, rbins, period=period, - num_threads=num_threads, approx_cell1_size=approx_cell1_size)[:, 0] + result = npairs_per_object_3d( + sample, + tracers, + rbins, + period=period, + num_threads=num_threads, + approx_cell1_size=approx_cell1_size, + )[:, 0] - environment_volume = (4/3.)*np.pi*radius**3 - number_density = result/environment_volume + environment_volume = (4 / 3.0) * np.pi * radius**3 + number_density = result / environment_volume if norm_by_mean_density is True: - mean_rho = tracers.shape[0]/sample_volume - return number_density/mean_rho + mean_rho = tracers.shape[0] / sample_volume + return number_density / mean_rho else: return number_density def _large_scale_density_spherical_volume_process_args( - sample, tracers, radius, period, sample_volume, num_threads, approx_cell1_size): - """ - """ + sample, tracers, radius, period, sample_volume, num_threads, approx_cell1_size +): + """ """ sample = np.atleast_1d(sample) tracers = np.atleast_1d(tracers) rbins = np.atleast_1d(radius).astype(float) - rbins = np.append(rbins, rbins[0]+0.0001) - - if period is None: - if sample_volume is None: - msg = ("If period is None, you must pass in ``sample_volume``.") - raise HalotoolsError(msg) - else: - sample_volume = float(sample_volume) + rbins = np.append(rbins, rbins[0] + 0.0001) + + period = get_period(period)[0] + + if sample_volume is None: + sample_volume = period.prod() else: - period = np.atleast_1d(period) - if len(period) == 1: - period = np.array([period, period, period]) - elif len(period) == 3: - pass - else: - msg = ("\nInput ``period`` must either be a float or length-3 sequence") - raise HalotoolsError(msg) - if sample_volume is None: - sample_volume = period.prod() - else: - msg = ("If period is not None, do not pass in sample_volume") - raise HalotoolsError(msg) + msg = "If period is not None, do not pass in sample_volume" + raise HalotoolsError(msg) return sample, tracers, rbins, period, sample_volume, num_threads, approx_cell1_size diff --git a/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_volume.py b/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_volume.py index f92e1238a..6e21fcf60 100644 --- a/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_volume.py +++ b/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_volume.py @@ -1,72 +1,76 @@ -""" -""" -from __future__ import (absolute_import, division, print_function) +""" """ + +from __future__ import absolute_import, division, print_function + import numpy as np import pytest -from ..large_scale_density_spherical_volume import large_scale_density_spherical_volume - -from ...tests.cf_helpers import generate_locus_of_3d_points from ....custom_exceptions import HalotoolsError +from ...tests.cf_helpers import generate_locus_of_3d_points +from ..large_scale_density_spherical_volume import large_scale_density_spherical_volume -__all__ = ('test_large_scale_density_spherical_volume1', ) +__all__ = ("test_large_scale_density_spherical_volume1",) fixed_seed = 43 def test_large_scale_density_spherical_volume_exception_handling(): - """ - """ + """ """ npts1, npts2 = 100, 200 sample = generate_locus_of_3d_points(npts1, xc=0.1, yc=0.1, zc=0.1, seed=fixed_seed) - tracers = generate_locus_of_3d_points(npts2, xc=0.15, yc=0.1, zc=0.1, seed=fixed_seed) + tracers = generate_locus_of_3d_points( + npts2, xc=0.15, yc=0.1, zc=0.1, seed=fixed_seed + ) radius = 0.1 with pytest.raises(HalotoolsError) as err: - result = large_scale_density_spherical_volume( - sample, tracers, radius) + result = large_scale_density_spherical_volume(sample, tracers, radius) substr = "If period is None, you must pass in ``sample_volume``." assert substr in err.value.args[0] with pytest.raises(HalotoolsError) as err: result = large_scale_density_spherical_volume( - sample, tracers, radius, period=[1, 1]) + sample, tracers, radius, period=[1, 1] + ) substr = "Input ``period`` must either be a float or length-3 sequence" assert substr in err.value.args[0] with pytest.raises(HalotoolsError) as err: result = large_scale_density_spherical_volume( - sample, tracers, radius, period=1, sample_volume=0.4) + sample, tracers, radius, period=1, sample_volume=0.4 + ) substr = "If period is not None, do not pass in sample_volume" assert substr in err.value.args[0] def test_large_scale_density_spherical_volume1(): - """ - """ + """ """ npts1, npts2 = 100, 200 sample = generate_locus_of_3d_points(npts1, xc=0.1, yc=0.1, zc=0.1, seed=fixed_seed) - tracers = generate_locus_of_3d_points(npts2, xc=0.15, yc=0.1, zc=0.1, seed=fixed_seed) + tracers = generate_locus_of_3d_points( + npts2, xc=0.15, yc=0.1, zc=0.1, seed=fixed_seed + ) radius = 0.1 - result = large_scale_density_spherical_volume( - sample, tracers, radius, period=1) + result = large_scale_density_spherical_volume(sample, tracers, radius, period=1) - environment_volume = (4/3.)*np.pi*radius**3 - correct_answer = 200/environment_volume + environment_volume = (4 / 3.0) * np.pi * radius**3 + correct_answer = 200 / environment_volume assert np.allclose(result, correct_answer, rtol=0.001) def test_large_scale_density_spherical_volume2(): - """ - """ + """ """ npts1, npts2 = 100, 200 sample = generate_locus_of_3d_points(npts1, xc=0.1, yc=0.1, zc=0.1, seed=fixed_seed) - tracers = generate_locus_of_3d_points(npts2, xc=0.95, yc=0.1, zc=0.1, seed=fixed_seed) + tracers = generate_locus_of_3d_points( + npts2, xc=0.95, yc=0.1, zc=0.1, seed=fixed_seed + ) radius = 0.2 result = large_scale_density_spherical_volume( - sample, tracers, radius, period=[1, 1, 1], norm_by_mean_density=True) + sample, tracers, radius, period=[1, 1, 1], norm_by_mean_density=True + ) mean_density = float(npts2) - environment_volume = (4/3.)*np.pi*radius**3 - correct_answer = 200/environment_volume/mean_density + environment_volume = (4 / 3.0) * np.pi * radius**3 + correct_answer = 200 / environment_volume / mean_density assert np.allclose(result, correct_answer, rtol=0.001) diff --git a/halotools/mock_observables/mock_observables_helpers.py b/halotools/mock_observables/mock_observables_helpers.py index eb2a4754c..3991f9a6f 100644 --- a/halotools/mock_observables/mock_observables_helpers.py +++ b/halotools/mock_observables/mock_observables_helpers.py @@ -1,23 +1,29 @@ -""" Module containing various helper functions used to process the +"""Module containing various helper functions used to process the arguments of functions throughout the `~halotools.mock_observables` package. """ + from __future__ import absolute_import, division, print_function, unicode_literals +import multiprocessing from warnings import warn + import numpy as np -import multiprocessing from ..utils.array_utils import array_is_monotonic - num_available_cores = multiprocessing.cpu_count() -__all__ = ('enforce_sample_respects_pbcs', 'get_num_threads', 'get_period', - 'enforce_sample_has_correct_shape', 'get_separation_bins_array') +__all__ = ( + "enforce_sample_respects_pbcs", + "get_num_threads", + "get_period", + "enforce_sample_has_correct_shape", + "get_separation_bins_array", +) def enforce_sample_respects_pbcs(x, y, z, period): - """ Verify that the input sample is properly bounded in all dimensions by the input period. + """Verify that the input sample is properly bounded in all dimensions by the input period. Parameters ----------- @@ -30,57 +36,67 @@ def enforce_sample_respects_pbcs(x, y, z, period): assert np.all(y >= 0) assert np.all(z >= 0) except: - msg = ("You set periodic boundary conditions to be True by passing in \n" + msg = ( + "You set periodic boundary conditions to be True by passing in \n" "period = (%.2f, %.2f, %.2f), but your input data has negative values,\n" - "indicating that you forgot to apply periodic boundary conditions.\n") + "indicating that you forgot to apply periodic boundary conditions.\n" + ) raise ValueError(msg % (period[0], period[1], period[2])) try: assert np.all(x <= period[0]) except: - msg = ("You set xperiod = %.2f but there are values in the x-dimension \n" - "of the input data that exceed this value") + msg = ( + "You set xperiod = %.2f but there are values in the x-dimension \n" + "of the input data that exceed this value" + ) raise ValueError(msg % period[0]) try: assert np.all(y <= period[1]) except: - msg = ("You set yperiod = %.2f but there are values in the y-dimension \n" - "of the input data that exceed this value") + msg = ( + "You set yperiod = %.2f but there are values in the y-dimension \n" + "of the input data that exceed this value" + ) raise ValueError(msg % period[1]) try: assert np.all(z <= period[2]) except: - msg = ("You set zperiod = %.2f but there are values in the z-dimension \n" - "of the input data that exceed this value") + msg = ( + "You set zperiod = %.2f but there are values in the z-dimension \n" + "of the input data that exceed this value" + ) raise ValueError(msg % period[2]) def get_num_threads(input_num_threads, enforce_max_cores=False): - """ Helper function requires that ``input_num_threads`` either be an + """Helper function requires that ``input_num_threads`` either be an integer or the string ``max``. If ``input_num_threads`` exceeds the number of available cores, a warning will be issued. In this event, ``enforce_max_cores`` is set to True, then ``num_threads`` is automatically set to num_cores. """ - if input_num_threads == 'max': + if input_num_threads == "max": num_threads = num_available_cores else: try: num_threads = int(input_num_threads) assert num_threads == input_num_threads except: - msg = ("Input ``num_threads`` must be an integer") + msg = "Input ``num_threads`` must be an integer" raise ValueError(msg) if num_threads > num_available_cores: - msg = ("Input ``num_threads`` = {0} exceeds the ``num_available_cores`` = {1}.\n") + msg = "Input ``num_threads`` = {0} exceeds the ``num_available_cores`` = {1}.\n" warn(msg.format(num_threads, num_available_cores)) if enforce_max_cores is True: - msg = ("Since ``enforce_max_cores`` is True,\n" - "setting ``num_threads`` to ``num_available_cores``.") + msg = ( + "Since ``enforce_max_cores`` is True,\n" + "setting ``num_threads`` to ``num_available_cores``." + ) warn(msg) num_threads = num_available_cores @@ -89,7 +105,7 @@ def get_num_threads(input_num_threads, enforce_max_cores=False): def get_period(period): - """ Helper function used to process the input ``period`` argument. + """Helper function used to process the input ``period`` argument. If ``period`` is set to None, function returns period, PBCs = (None, False). Otherwise, function returns ([period, period, period], True). """ @@ -101,38 +117,41 @@ def get_period(period): period = np.atleast_1d(period).astype(float) if len(period) == 1: - period = np.array([period[0]]*3).astype(float) + period = np.array([period[0]] * 3).astype(float) try: assert np.all(period < np.inf) assert np.all(period > 0) assert len(period) == 3 except AssertionError: - msg = ("Input ``period`` must be either a scalar or a 3-element sequence.\n" - "All values must bounded positive numbers.\n") + msg = ( + "Input ``period`` must be either a scalar or a 3-element sequence.\n" + "All values must bounded positive numbers.\n" + ) raise ValueError(msg) return period, PBCs def enforce_sample_has_correct_shape(sample, ndim=3): - """ Function inspects the input ``sample`` and enforces that it is of shape (Npts, 3). - """ + """Function inspects the input ``sample`` and enforces that it is of shape (Npts, 3).""" sample = np.atleast_1d(sample) try: input_shape = np.shape(sample) assert len(input_shape) == 2 assert input_shape[1] == ndim except: - msg = ("Input sample of points must be a Numpy ndarray of shape (Npts, {0}).\n" + msg = ( + "Input sample of points must be a Numpy ndarray of shape (Npts, {0}).\n" "To convert a sequence of 1d arrays x, y, z into correct shape expected \n" "throughout the `mock_observables` package:\n\n" - ">>> sample = np.vstack([x, y, z]).T ".format(ndim)) + ">>> sample = np.vstack([x, y, z]).T ".format(ndim) + ) raise TypeError(msg) return sample def get_separation_bins_array(separation_bins): - """ Function verifies that the input ``separation_bins`` is a monotonically increasing + """Function verifies that the input ``separation_bins`` is a monotonically increasing 1d Numpy array with at least two entries, all of which are required to be strictly positive. This helper function can be used equally well with 3d separation bins ``rbins`` or @@ -149,15 +168,17 @@ def get_separation_bins_array(separation_bins): assert array_is_monotonic(separation_bins, strict=True) == 1 assert np.all(separation_bins > 0) except AssertionError: - msg = ("\n Input separation bins must be a monotonically increasing \n" - "1-D array with at least two entries, all of which must be strictly positive.\n") + msg = ( + "\n Input separation bins must be a monotonically increasing \n" + "1-D array with at least two entries, all of which must be strictly positive.\n" + ) raise TypeError(msg) return separation_bins def get_line_of_sight_bins_array(pi_bins): - """ Function verifies that the input ``pi_bins`` is a monotonically increasing + """Function verifies that the input ``pi_bins`` is a monotonically increasing 1d Numpy array with at least two entries. The `get_line_of_sight_bins_array` function differs from the `get_separation_bins_array` function only in that values of zero are permissible. @@ -170,8 +191,10 @@ def get_line_of_sight_bins_array(pi_bins): if len(pi_bins) > 2: assert array_is_monotonic(pi_bins, strict=True) == 1 except AssertionError: - msg = ("\n Input separation bins must be a monotonically increasing \n" - "1-D array with at least two entries.\n") + msg = ( + "\n Input separation bins must be a monotonically increasing \n" + "1-D array with at least two entries.\n" + ) raise TypeError(msg) return pi_bins diff --git a/halotools/mock_observables/pair_counters/rectangular_mesh.py b/halotools/mock_observables/pair_counters/rectangular_mesh.py index 31c05acc6..c228f2a69 100644 --- a/halotools/mock_observables/pair_counters/rectangular_mesh.py +++ b/halotools/mock_observables/pair_counters/rectangular_mesh.py @@ -1,70 +1,89 @@ -""" Module containing `~halotools.mock_observables.RectangularDoubleMesh`, +"""Module containing `~halotools.mock_observables.RectangularDoubleMesh`, the primary data structure used to optimize pairwise calculations throughout the `~halotools.mock_observables` sub-package. """ -import numpy as np + from math import floor -__all__ = ('RectangularDoubleMesh', ) -__author__ = ('Andrew Hearin', ) +import numpy as np + +from ..mock_observables_helpers import get_period + +__all__ = ("RectangularDoubleMesh",) +__author__ = ("Andrew Hearin",) default_max_cells_per_dimension_cell1 = 50 default_max_cells_per_dimension_cell2 = 50 def digitized_position(p, cell_size, num_divs): - """ Function returns a discretized spatial position of input point(s). - """ + """Function returns a discretized spatial position of input point(s).""" ip = np.floor(p // cell_size).astype(int) - return np.where(ip >= num_divs, num_divs-1, ip) + return np.where(ip >= num_divs, num_divs - 1, ip) -def sample1_cell_size(period, search_length, approx_cell_size, - max_cells_per_dimension=default_max_cells_per_dimension_cell1): - """ Function determines the size of the cells of mesh1. +def sample1_cell_size( + period, + search_length, + approx_cell_size, + max_cells_per_dimension=default_max_cells_per_dimension_cell1, +): + """Function determines the size of the cells of mesh1. The conditions that must be met are that the cell size must be less than the search length, must evenly divide the box length, and may not exceed ``max_cells_per_dimension``. """ - if search_length > period/3.: - msg = ("Input ``search_length`` cannot exceed period/3") + period = get_period(period)[0] + assert np.allclose(period, period[0]) + period = float(period[0]) + + if search_length > period / 3.0: + msg = "Input ``search_length`` cannot exceed period/3" raise ValueError(msg) - ndivs = int(floor(period/float(approx_cell_size))) + ndivs = int(floor(period / float(approx_cell_size))) ndivs = max(ndivs, 1) ndivs = min(max_cells_per_dimension, ndivs) - nsearch = int(floor(period/float(search_length))) + nsearch = int(floor(period / float(search_length))) nsearch = max(nsearch, 1) ndivs = min(ndivs, nsearch) ndivs = max(3, ndivs) - cell_size = period/float(ndivs) + cell_size = period / float(ndivs) return cell_size -def sample2_cell_sizes(period, sample1_cell_size, approx_cell_size, - max_cells_per_dimension=default_max_cells_per_dimension_cell2): - """ Function determines the size of the cells of mesh2. +def sample2_cell_sizes( + period, + sample1_cell_size, + approx_cell_size, + max_cells_per_dimension=default_max_cells_per_dimension_cell2, +): + """Function determines the size of the cells of mesh2. The conditions that must be met are that the cell size must be less than the search length, must evenly divide the box length, and may not exceed ``max_cells_per_dimension``. """ + period = get_period(period)[0] + assert np.allclose(period, period[0]) + period = float(period[0]) + num_sample1_cells = int(np.round(period / sample1_cell_size)) - ndivs_sample1_cells = int(np.round(sample1_cell_size/float(approx_cell_size))) + ndivs_sample1_cells = int(np.round(sample1_cell_size / float(approx_cell_size))) ndivs_sample1_cells = max(1, ndivs_sample1_cells) ndivs_sample1_cells = min(max_cells_per_dimension, ndivs_sample1_cells) - num_sample2_cells = num_sample1_cells*ndivs_sample1_cells + num_sample2_cells = num_sample1_cells * ndivs_sample1_cells if num_sample2_cells > max_cells_per_dimension: num2_per_num1 = max_cells_per_dimension // num_sample1_cells - num_sample2_cells = num2_per_num1*num_sample1_cells - cell_size = period/float(num_sample2_cells) + num_sample2_cells = num2_per_num1 * num_sample1_cells + cell_size = period / float(num_sample2_cells) return cell_size class RectangularMesh(object): - """ Underlying mesh structure used to place points into rectangular cells + """Underlying mesh structure used to place points into rectangular cells within a simulation volume. The simulation box is divided into rectanguloid cells whose edges @@ -96,8 +115,18 @@ class RectangularMesh(object): """ - def __init__(self, x1in, y1in, z1in, xperiod, yperiod, zperiod, - approx_xcell_size, approx_ycell_size, approx_zcell_size): + def __init__( + self, + x1in, + y1in, + z1in, + xperiod, + yperiod, + zperiod, + approx_xcell_size, + approx_ycell_size, + approx_zcell_size, + ): """ Parameters ---------- @@ -151,6 +180,21 @@ def __init__(self, x1in, y1in, z1in, xperiod, yperiod, zperiod, self.npts = x1in.shape[0] + try: + xperiod = float(xperiod[0]) + except IndexError: + xperiod = float(xperiod) + + try: + yperiod = float(yperiod[0]) + except IndexError: + yperiod = float(yperiod) + + try: + zperiod = float(zperiod[0]) + except IndexError: + zperiod = float(zperiod) + self.xperiod = xperiod self.yperiod = yperiod self.zperiod = zperiod @@ -158,7 +202,7 @@ def __init__(self, x1in, y1in, z1in, xperiod, yperiod, zperiod, self.num_xdivs = max(int(np.round(xperiod / approx_xcell_size)), 1) self.num_ydivs = max(int(np.round(yperiod / approx_ycell_size)), 1) self.num_zdivs = max(int(np.round(zperiod / approx_zcell_size)), 1) - self.ncells = self.num_xdivs*self.num_ydivs*self.num_zdivs + self.ncells = self.num_xdivs * self.num_ydivs * self.num_zdivs self.xcell_size = self.xperiod / float(self.num_xdivs) self.ycell_size = self.yperiod / float(self.num_ydivs) @@ -171,28 +215,46 @@ def __init__(self, x1in, y1in, z1in, xperiod, yperiod, zperiod, cell_ids = self.cell_id_from_cell_tuple(ix, iy, iz) self.idx_sorted = np.ascontiguousarray(np.argsort(cell_ids)) - cell_id_indices = np.searchsorted(cell_ids, np.arange(self.ncells), - sorter=self.idx_sorted) + cell_id_indices = np.searchsorted( + cell_ids, np.arange(self.ncells), sorter=self.idx_sorted + ) cell_id_indices = np.append(cell_id_indices, self.npts) self.cell_id_indices = np.ascontiguousarray(cell_id_indices) def cell_id_from_cell_tuple(self, ix, iy, iz): - return ix*(self.num_ydivs*self.num_zdivs) + iy*self.num_zdivs + iz + return ix * (self.num_ydivs * self.num_zdivs) + iy * self.num_zdivs + iz class RectangularDoubleMesh(object): - """ Fundamental data structure of the `~halotools.mock_observables` sub-package. + """Fundamental data structure of the `~halotools.mock_observables` sub-package. `~halotools.mock_observables.RectangularDoubleMesh` is built up from two instances of `~halotools.mock_observables.pair_counters.rectangular_mesh.RectangularMesh`. """ - def __init__(self, x1, y1, z1, x2, y2, z2, - approx_x1cell_size, approx_y1cell_size, approx_z1cell_size, - approx_x2cell_size, approx_y2cell_size, approx_z2cell_size, - search_xlength, search_ylength, search_zlength, - xperiod, yperiod, zperiod, PBCs=True, - max_cells_per_dimension_cell1=default_max_cells_per_dimension_cell1, - max_cells_per_dimension_cell2=default_max_cells_per_dimension_cell2): + def __init__( + self, + x1, + y1, + z1, + x2, + y2, + z2, + approx_x1cell_size, + approx_y1cell_size, + approx_z1cell_size, + approx_x2cell_size, + approx_y2cell_size, + approx_z2cell_size, + search_xlength, + search_ylength, + search_zlength, + xperiod, + yperiod, + zperiod, + PBCs=True, + max_cells_per_dimension_cell1=default_max_cells_per_dimension_cell1, + max_cells_per_dimension_cell2=default_max_cells_per_dimension_cell2, + ): """ Parameters ---------- @@ -247,23 +309,65 @@ def __init__(self, x1, y1, z1, x2, y2, z2, self._check_sensible_constructor_inputs() - approx_x1cell_size = sample1_cell_size(xperiod, search_xlength, approx_x1cell_size, - max_cells_per_dimension=max_cells_per_dimension_cell1) - approx_y1cell_size = sample1_cell_size(yperiod, search_ylength, approx_y1cell_size, - max_cells_per_dimension=max_cells_per_dimension_cell1) - approx_z1cell_size = sample1_cell_size(zperiod, search_zlength, approx_z1cell_size, - max_cells_per_dimension=max_cells_per_dimension_cell1) - self.mesh1 = RectangularMesh(x1, y1, z1, xperiod, yperiod, zperiod, - approx_x1cell_size, approx_y1cell_size, approx_z1cell_size) - - approx_x2cell_size = sample2_cell_sizes(xperiod, self.mesh1.xcell_size, approx_x2cell_size, - max_cells_per_dimension=max_cells_per_dimension_cell2) - approx_y2cell_size = sample2_cell_sizes(yperiod, self.mesh1.ycell_size, approx_y2cell_size, - max_cells_per_dimension=max_cells_per_dimension_cell2) - approx_z2cell_size = sample2_cell_sizes(zperiod, self.mesh1.zcell_size, approx_z2cell_size, - max_cells_per_dimension=max_cells_per_dimension_cell2) - self.mesh2 = RectangularMesh(x2, y2, z2, xperiod, yperiod, zperiod, - approx_x2cell_size, approx_y2cell_size, approx_z2cell_size) + approx_x1cell_size = sample1_cell_size( + xperiod, + search_xlength, + approx_x1cell_size, + max_cells_per_dimension=max_cells_per_dimension_cell1, + ) + approx_y1cell_size = sample1_cell_size( + yperiod, + search_ylength, + approx_y1cell_size, + max_cells_per_dimension=max_cells_per_dimension_cell1, + ) + approx_z1cell_size = sample1_cell_size( + zperiod, + search_zlength, + approx_z1cell_size, + max_cells_per_dimension=max_cells_per_dimension_cell1, + ) + self.mesh1 = RectangularMesh( + x1, + y1, + z1, + xperiod, + yperiod, + zperiod, + approx_x1cell_size, + approx_y1cell_size, + approx_z1cell_size, + ) + + approx_x2cell_size = sample2_cell_sizes( + xperiod, + self.mesh1.xcell_size, + approx_x2cell_size, + max_cells_per_dimension=max_cells_per_dimension_cell2, + ) + approx_y2cell_size = sample2_cell_sizes( + yperiod, + self.mesh1.ycell_size, + approx_y2cell_size, + max_cells_per_dimension=max_cells_per_dimension_cell2, + ) + approx_z2cell_size = sample2_cell_sizes( + zperiod, + self.mesh1.zcell_size, + approx_z2cell_size, + max_cells_per_dimension=max_cells_per_dimension_cell2, + ) + self.mesh2 = RectangularMesh( + x2, + y2, + z2, + xperiod, + yperiod, + zperiod, + approx_x2cell_size, + approx_y2cell_size, + approx_z2cell_size, + ) self.num_xcell2_per_xcell1 = self.mesh2.num_xdivs // self.mesh1.num_xdivs self.num_ycell2_per_ycell1 = self.mesh2.num_ydivs // self.mesh1.num_ydivs @@ -271,34 +375,43 @@ def __init__(self, x1, y1, z1, x2, y2, z2, def _check_sensible_constructor_inputs(self): try: - assert self.search_xlength <= self.xperiod/3. + assert self.search_xlength <= self.xperiod / 3.0 except AssertionError: - msg = ("\n The maximum length over which you search for pairs of points \n" + msg = ( + "\n The maximum length over which you search for pairs of points \n" "cannot be larger than Lbox/3 in any dimension. \n" "You tried to search for pairs out to a length of search_xlength = %.2f,\n" "but the size of your box in this dimension is xperiod = %.2f.\n" "If you need to count pairs on these length scales, \n" - "you should use a larger simulation.\n" % (self.search_xlength, self.xperiod)) + "you should use a larger simulation.\n" + % (self.search_xlength, self.xperiod) + ) raise ValueError(msg) try: - assert self.search_ylength <= self.yperiod/3. + assert self.search_ylength <= self.yperiod / 3.0 except AssertionError: - msg = ("\n The maximum length over which you search for pairs of points \n" + msg = ( + "\n The maximum length over which you search for pairs of points \n" "cannot be larger than Lbox/3 in any dimension. \n" "You tried to search for pairs out to a length of search_ylength = %.2f,\n" "but the size of your box in this dimension is yperiod = %.2f.\n" "If you need to count pairs on these length scales, \n" - "you should use a larger simulation.\n" % (self.search_ylength, self.yperiod)) + "you should use a larger simulation.\n" + % (self.search_ylength, self.yperiod) + ) raise ValueError(msg) try: - assert self.search_zlength <= self.zperiod/3. + assert self.search_zlength <= self.zperiod / 3.0 except AssertionError: - msg = ("\n The maximum length over which you search for pairs of points \n" + msg = ( + "\n The maximum length over which you search for pairs of points \n" "cannot be larger than Lbox/3 in any dimension. \n" "You tried to search for pairs out to a length of search_zlength = %.2f,\n" "but the size of your box in this dimension is zperiod = %.2f.\n" "If you need to count pairs on these length scales, \n" - "you should use a larger simulation.\n" % (self.search_zlength, self.zperiod)) + "you should use a larger simulation.\n" + % (self.search_zlength, self.zperiod) + ) raise ValueError(msg) From b81af864f2623f2d0d828e6eeb9768ad59cd0b94 Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 13:54:14 -0600 Subject: [PATCH 04/15] Resolve period bug in large_scale_density_spherical_annulus.py --- .../large_scale_density_spherical_annulus.py | 97 +++++++++++-------- 1 file changed, 55 insertions(+), 42 deletions(-) diff --git a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_annulus.py b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_annulus.py index c753000c4..958694849 100644 --- a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_annulus.py +++ b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_annulus.py @@ -2,25 +2,33 @@ Module containing functions used to determine various ways of quantifying the large-scale density of a set of points. """ + from __future__ import absolute_import, division, print_function, unicode_literals import numpy as np -from ..pair_counters import npairs_per_object_3d - from ...custom_exceptions import HalotoolsError +from ..mock_observables_helpers import get_period +from ..pair_counters import npairs_per_object_3d +__all__ = ("large_scale_density_spherical_annulus",) -__all__ = ('large_scale_density_spherical_annulus', ) - -__author__ = ('Andrew Hearin', ) +__author__ = ("Andrew Hearin",) -np.seterr(divide='ignore', invalid='ignore') # ignore divide by zero in e.g. DD/RR +np.seterr(divide="ignore", invalid="ignore") # ignore divide by zero in e.g. DD/RR -def large_scale_density_spherical_annulus(sample, tracers, inner_radius, outer_radius, - period=None, sample_volume=None, num_threads=1, approx_cell1_size=None, - norm_by_mean_density=False): +def large_scale_density_spherical_annulus( + sample, + tracers, + inner_radius, + outer_radius, + period=None, + sample_volume=None, + num_threads=1, + approx_cell1_size=None, + norm_by_mean_density=False, +): """ Calculate the mean density of the input ``sample`` from an input set of tracer particles using a spherical annulus @@ -101,58 +109,63 @@ def large_scale_density_spherical_annulus(sample, tracers, inner_radius, outer_r """ sample, tracers, rbins, period, sample_volume, num_threads, approx_cell1_size = ( _large_scale_density_spherical_annulus_process_args( - sample, tracers, inner_radius, outer_radius, - period, sample_volume, num_threads, approx_cell1_size) + sample, + tracers, + inner_radius, + outer_radius, + period, + sample_volume, + num_threads, + approx_cell1_size, ) - - result = npairs_per_object_3d(sample, tracers, rbins, period=period, - num_threads=num_threads, approx_cell1_size=approx_cell1_size) + ) + + result = npairs_per_object_3d( + sample, + tracers, + rbins, + period=period, + num_threads=num_threads, + approx_cell1_size=approx_cell1_size, + ) result = np.diff(result, axis=1) - environment_volume = (4/3.)*np.pi*(outer_radius**3 - inner_radius**3) - number_density = result/environment_volume + environment_volume = (4 / 3.0) * np.pi * (outer_radius**3 - inner_radius**3) + number_density = result / environment_volume if norm_by_mean_density is True: - mean_rho = tracers.shape[0]/sample_volume - return number_density/mean_rho + mean_rho = tracers.shape[0] / sample_volume + return number_density / mean_rho else: return number_density def _large_scale_density_spherical_annulus_process_args( - sample, tracers, inner_radius, outer_radius, - period, sample_volume, num_threads, approx_cell1_size): - """ - """ + sample, + tracers, + inner_radius, + outer_radius, + period, + sample_volume, + num_threads, + approx_cell1_size, +): + """ """ sample = np.atleast_1d(sample) tracers = np.atleast_1d(tracers) + period = get_period(period)[0] try: assert outer_radius > inner_radius except AssertionError: - msg = ("Input ``outer_radius`` must be larger than input ``inner_radius``") + msg = "Input ``outer_radius`` must be larger than input ``inner_radius``" raise HalotoolsError(msg) rbins = np.array([inner_radius, outer_radius]) - if period is None: - if sample_volume is None: - msg = ("If period is None, you must pass in ``sample_volume``.") - raise HalotoolsError(msg) - else: - sample_volume = float(sample_volume) + if sample_volume is None: + sample_volume = period.prod() else: - period = np.atleast_1d(period) - if len(period) == 1: - period = np.array([period, period, period]) - elif len(period) == 3: - pass - else: - msg = ("\nInput ``period`` must either be a float or length-3 sequence") - raise HalotoolsError(msg) - if sample_volume is None: - sample_volume = period.prod() - else: - msg = ("If period is not None, do not pass in sample_volume") - raise HalotoolsError(msg) + msg = "If period is not None, do not pass in sample_volume" + raise HalotoolsError(msg) return sample, tracers, rbins, period, sample_volume, num_threads, approx_cell1_size From 4ee028bcb8ed7156646b797337998fd56911c965 Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 14:10:57 -0600 Subject: [PATCH 05/15] Fix error handling in large_scale_density_spherical_volume.py --- .../large_scale_density_spherical_volume.py | 13 +++++++------ .../test_large_scale_density_spherical_volume.py | 9 +-------- 2 files changed, 8 insertions(+), 14 deletions(-) diff --git a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py index f569f1bd5..cc4c9683a 100644 --- a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py +++ b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py @@ -145,12 +145,13 @@ def _large_scale_density_spherical_volume_process_args( rbins = np.atleast_1d(radius).astype(float) rbins = np.append(rbins, rbins[0] + 0.0001) - period = get_period(period)[0] - if sample_volume is None: - sample_volume = period.prod() - else: - msg = "If period is not None, do not pass in sample_volume" - raise HalotoolsError(msg) + if period is None: + raise HalotoolsError("If sample_volume is None, must pass period") + if period is not None: + raise HalotoolsError("If period is not None, do not pass in sample_volume") + + period = get_period(period)[0] + sample_volume = period.prod() return sample, tracers, rbins, period, sample_volume, num_threads, approx_cell1_size diff --git a/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_volume.py b/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_volume.py index 6e21fcf60..ed9459dc3 100644 --- a/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_volume.py +++ b/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_volume.py @@ -25,14 +25,7 @@ def test_large_scale_density_spherical_volume_exception_handling(): with pytest.raises(HalotoolsError) as err: result = large_scale_density_spherical_volume(sample, tracers, radius) - substr = "If period is None, you must pass in ``sample_volume``." - assert substr in err.value.args[0] - - with pytest.raises(HalotoolsError) as err: - result = large_scale_density_spherical_volume( - sample, tracers, radius, period=[1, 1] - ) - substr = "Input ``period`` must either be a float or length-3 sequence" + substr = "must pass" assert substr in err.value.args[0] with pytest.raises(HalotoolsError) as err: From 81c60cc1433d470761d859f97d66e124dfddfac9 Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 14:17:31 -0600 Subject: [PATCH 06/15] Fix error handling in void_prob_func.py --- .../mock_observables/pair_counters/pairs.py | 113 ++++++++++-------- .../test_marked_npairs_3d.py | 13 +- .../tests/test_underdensity_prob_func.py | 8 +- .../tests/test_void_prob_func.py | 25 +--- .../void_statistics/void_prob_func.py | 9 +- 5 files changed, 83 insertions(+), 85 deletions(-) diff --git a/halotools/mock_observables/pair_counters/pairs.py b/halotools/mock_observables/pair_counters/pairs.py index 745f681d8..32082b8ae 100755 --- a/halotools/mock_observables/pair_counters/pairs.py +++ b/halotools/mock_observables/pair_counters/pairs.py @@ -5,10 +5,11 @@ """ from __future__ import absolute_import, division, print_function, unicode_literals + import numpy as np -__all__ = ['npairs', 'wnpairs', 'xy_z_npairs', 'xy_z_wnpairs', 's_mu_npairs'] -__author__ = ['Duncan Campbell'] +__all__ = ["npairs", "wnpairs", "xy_z_npairs", "xy_z_wnpairs", "s_mu_npairs"] +__author__ = ["Duncan Campbell"] def npairs(sample1, sample2, rbins, period=None): @@ -51,22 +52,24 @@ def npairs(sample1, sample2, rbins, period=None): # Process period entry and check for consistency. if period is None: - period = np.array([np.inf]*np.shape(sample1)[-1]) + period = np.array([np.inf] * np.shape(sample1)[-1]) else: period = np.asarray(period).astype("float64") if np.shape(period) == (): - period = np.array([period]*np.shape(sample1)[-1]) + period = np.array([period] * np.shape(sample1)[-1]) elif np.shape(period)[0] != np.shape(sample1)[-1]: raise ValueError("period should have len == dimension of points") return None N1 = len(sample1) N2 = len(sample2) - dd = np.zeros((N1*N2,)) # store radial pair separations - for i in range(0, N1): # calculate distance between every point and every other point + dd = np.zeros((N1 * N2,)) # store radial pair separations + for i in range( + 0, N1 + ): # calculate distance between every point and every other point x1 = sample1[i, :] x2 = sample2 - dd[i*N2:i*N2+N2] = distance(x1, x2, period) + dd[i * N2 : i * N2 + N2] = distance(x1, x2, period) # sort results dd.sort() @@ -127,23 +130,25 @@ def xy_z_npairs(sample1, sample2, rp_bins, pi_bins, period=None): # Process period entry and check for consistency. if period is None: - period = np.array([np.inf]*np.shape(sample1)[-1]) + period = np.array([np.inf] * np.shape(sample1)[-1]) else: period = np.asarray(period).astype("float64") if np.shape(period) == (): - period = np.array([period]*np.shape(sample1)[-1]) + period = np.array([period] * np.shape(sample1)[-1]) elif np.shape(period)[0] != np.shape(sample1)[-1]: raise ValueError("period should have len == dimension of points") return None N1 = len(sample1) N2 = len(sample2) - dd = np.zeros((N1*N2, 2)) # store pair separations - for i in range(0, N1): # calculate distance between every point and every other point + dd = np.zeros((N1 * N2, 2)) # store pair separations + for i in range( + 0, N1 + ): # calculate distance between every point and every other point x1 = sample1[i, :] x2 = sample2 - dd[i*N2:i*N2+N2, 1] = parallel_distance(x1, x2, period) - dd[i*N2:i*N2+N2, 0] = perpendicular_distance(x1, x2, period) + dd[i * N2 : i * N2 + N2, 1] = parallel_distance(x1, x2, period) + dd[i * N2 : i * N2 + N2, 0] = perpendicular_distance(x1, x2, period) # count number less than r n = np.zeros((rp_bins.size, pi_bins.size), dtype=np.int64) @@ -200,18 +205,18 @@ def wnpairs(sample1, sample2, r, period=None, weights1=None, weights2=None): # Process period entry and check for consistency. if period is None: - period = np.array([np.inf]*np.shape(sample1)[-1]) + period = np.array([np.inf] * np.shape(sample1)[-1]) else: period = np.asarray(period).astype("float64") if np.shape(period) == (): - period = np.array([period]*np.shape(sample1)[-1]) + period = np.array([period] * np.shape(sample1)[-1]) if np.shape(period)[0] != np.shape(sample1)[-1]: raise ValueError("period should have len == dimension of points") return None # Process weights1 entry and check for consistency. if weights1 is None: - weights1 = np.array([1.0]*np.shape(sample1)[0], dtype=np.float64) + weights1 = np.array([1.0] * np.shape(sample1)[0], dtype=np.float64) else: weights1 = np.asarray(weights1).astype("float64") if np.shape(weights1)[0] != np.shape(sample1)[0]: @@ -219,7 +224,7 @@ def wnpairs(sample1, sample2, r, period=None, weights1=None, weights2=None): return None # Process weights2 entry and check for consistency. if weights2 is None: - weights2 = np.array([1.0]*np.shape(sample2)[0], dtype=np.float64) + weights2 = np.array([1.0] * np.shape(sample2)[0], dtype=np.float64) else: weights2 = np.asarray(weights2).astype("float64") if np.shape(weights2)[0] != np.shape(sample2)[0]: @@ -229,7 +234,9 @@ def wnpairs(sample1, sample2, r, period=None, weights1=None, weights2=None): N1 = len(sample1) N2 = len(sample2) dd = np.zeros((N1, N2), dtype=np.float64) # store radial pair separations - for i in range(0, N1): # calculate distance between every point and every other point + for i in range( + 0, N1 + ): # calculate distance between every point and every other point x1 = sample1[i, :] x2 = sample2 dd[i, :] = distance(x1, x2, period) @@ -238,12 +245,14 @@ def wnpairs(sample1, sample2, r, period=None, weights1=None, weights2=None): n = np.zeros((r.size,), dtype=np.float64) for i in range(r.size): for j in range(N1): - n[i] += np.sum(np.extract(dd[j, :] <= r[i], weights2))*weights1[j] + n[i] += np.sum(np.extract(dd[j, :] <= r[i], weights2)) * weights1[j] return n -def xy_z_wnpairs(sample1, sample2, rp_bins, pi_bins, period=None, weights1=None, weights2=None): +def xy_z_wnpairs( + sample1, sample2, rp_bins, pi_bins, period=None, weights1=None, weights2=None +): r""" Calculate the number of weighted pairs with parellal separations less than or equal to pi_bins[i], and perpendicular separations less than or equal to rp_bins[i]. @@ -298,18 +307,18 @@ def xy_z_wnpairs(sample1, sample2, rp_bins, pi_bins, period=None, weights1=None, # Process period entry and check for consistency. if period is None: - period = np.array([np.inf]*np.shape(sample1)[-1]) + period = np.array([np.inf] * np.shape(sample1)[-1]) else: period = np.asarray(period).astype("float64") if np.shape(period) == (): - period = np.array([period]*np.shape(sample1)[-1]) + period = np.array([period] * np.shape(sample1)[-1]) elif np.shape(period)[0] != np.shape(sample1)[-1]: raise ValueError("period should have len == dimension of points") return None # Process weights1 entry and check for consistency. if weights1 is None: - weights1 = np.array([1.0]*np.shape(sample1)[0], dtype=np.float64) + weights1 = np.array([1.0] * np.shape(sample1)[0], dtype=np.float64) else: weights1 = np.asarray(weights1).astype("float64") if np.shape(weights1)[0] != np.shape(sample1)[0]: @@ -317,7 +326,7 @@ def xy_z_wnpairs(sample1, sample2, rp_bins, pi_bins, period=None, weights1=None, return None # Process weights2 entry and check for consistency. if weights2 is None: - weights2 = np.array([1.0]*np.shape(sample2)[0], dtype=np.float64) + weights2 = np.array([1.0] * np.shape(sample2)[0], dtype=np.float64) else: weights2 = np.asarray(weights2).astype("float64") if np.shape(weights2)[0] != np.shape(sample2)[0]: @@ -326,20 +335,24 @@ def xy_z_wnpairs(sample1, sample2, rp_bins, pi_bins, period=None, weights1=None, N1 = len(sample1) N2 = len(sample2) - dd = np.zeros((N1*N2, 2)) # store pair separations - ww = np.zeros((N1*N2, 1)) # store pair separations - for i in range(0, N1): # calculate distance between every point and every other point + dd = np.zeros((N1 * N2, 2)) # store pair separations + ww = np.zeros((N1 * N2, 1)) # store pair separations + for i in range( + 0, N1 + ): # calculate distance between every point and every other point x1 = sample1[i, :] x2 = sample2 - dd[i*N2:i*N2+N2, 1] = parallel_distance(x1, x2, period) - dd[i*N2:i*N2+N2, 0] = perpendicular_distance(x1, x2, period) - ww[i*N2:i*N2+N2] = weights1[i]*weights2 + dd[i * N2 : i * N2 + N2, 1] = parallel_distance(x1, x2, period) + dd[i * N2 : i * N2 + N2, 0] = perpendicular_distance(x1, x2, period) + ww[i * N2 : i * N2 + N2] = weights1[i] * weights2 # count number less than r n = np.zeros((rp_bins.size, pi_bins.size), dtype=np.float64) for i in range(rp_bins.size): for j in range(pi_bins.size): - n[i, j] += np.sum(np.extract((dd[:, 0] <= rp_bins[i]) & (dd[:, 1] <= pi_bins[j]), ww)) + n[i, j] += np.sum( + np.extract((dd[:, 0] <= rp_bins[i]) & (dd[:, 1] <= pi_bins[j]), ww) + ) return n @@ -400,11 +413,11 @@ def s_mu_npairs(sample1, sample2, s_bins, mu_bins, period=None): # Process period entry and check for consistency. if period is None: - period = np.array([np.inf]*np.shape(sample1)[-1]) + period = np.array([np.inf] * np.shape(sample1)[-1]) else: period = np.asarray(period).astype("float64") if np.shape(period) == (): - period = np.array([period]*np.shape(sample1)[-1]) + period = np.array([period] * np.shape(sample1)[-1]) elif np.shape(period)[0] != np.shape(sample1)[-1]: raise ValueError("period should have len == dimension of points") return None @@ -413,14 +426,14 @@ def s_mu_npairs(sample1, sample2, s_bins, mu_bins, period=None): # note that this array can be very large for large N1 and N2 N1 = len(sample1) N2 = len(sample2) - dd = np.zeros((N1*N2, 2)) + dd = np.zeros((N1 * N2, 2)) # calculate distance between every point and every other point for i in range(0, N1): x1 = sample1[i, :] x2 = sample2 - dd[i*N2:i*N2+N2, 0] = distance(x1, x2, period) - dd[i*N2:i*N2+N2, 1] = np.cos(theta_LOS(x1, x2, period)) + dd[i * N2 : i * N2 + N2, 0] = distance(x1, x2, period) + dd[i * N2 : i * N2 + N2, 1] = np.cos(theta_LOS(x1, x2, period)) # put mu bins in increasing theta_LOS order mu_bins = np.sort(mu_bins)[::-1] @@ -458,7 +471,7 @@ def distance(x1, x2, period=None): x1 = np.atleast_2d(x1) x2 = np.atleast_2d(x2) if period is None: - period = np.array([np.inf]*np.shape(x1)[-1]) + period = np.array([np.inf] * np.shape(x1)[-1]) # check for consistency if np.shape(x1)[-1] != np.shape(x2)[-1]: @@ -469,7 +482,7 @@ def distance(x1, x2, period=None): raise ValueError("period must have length equal to the dimension of x1 and x2.") m = np.minimum(np.fabs(x1 - x2), period - np.fabs(x1 - x2)) - distance = np.sqrt(np.sum(m*m, axis=len(np.shape(m))-1)) + distance = np.sqrt(np.sum(m * m, axis=len(np.shape(m)) - 1)) return distance @@ -500,7 +513,7 @@ def parallel_distance(x1, x2, period=None): x1 = np.atleast_2d(x1) x2 = np.atleast_2d(x2) if period is None: - period = np.array([np.inf]*np.shape(x1)[-1]) + period = np.array([np.inf] * np.shape(x1)[-1]) # check for consistency if np.shape(x1)[-1] != np.shape(x2)[-1]: @@ -510,8 +523,10 @@ def parallel_distance(x1, x2, period=None): if np.shape(period)[0] != np.shape(x1)[-1]: raise ValueError("period must have length equal to the dimension of x1 and x2.") - m = np.minimum(np.fabs(x1[:, -1] - x2[:, -1]), period[-1] - np.fabs(x1[:, -1] - x2[:, -1])) - distance = np.sqrt(m*m) + m = np.minimum( + np.fabs(x1[:, -1] - x2[:, -1]), period[-1] - np.fabs(x1[:, -1] - x2[:, -1]) + ) + distance = np.sqrt(m * m) return distance @@ -542,7 +557,7 @@ def perpendicular_distance(x1, x2, period=None): x1 = np.atleast_2d(x1) x2 = np.atleast_2d(x2) if period is None: - period = np.array([np.inf]*np.shape(x1)[-1]) + period = np.array([np.inf] * np.shape(x1)[-1]) # check for consistency if np.shape(x1)[-1] != np.shape(x2)[-1]: @@ -552,8 +567,10 @@ def perpendicular_distance(x1, x2, period=None): if np.shape(period)[0] != np.shape(x1)[-1]: raise ValueError("period must have length equal to the dimension of x1 and x2.") - m = np.minimum(np.fabs(x1[:, :-1] - x2[:, :-1]), period[:-1] - np.fabs(x1[:, :-1] - x2[:, :-1])) - distance = np.sqrt(np.sum(m*m, axis=len(np.shape(m))-1)) + m = np.minimum( + np.fabs(x1[:, :-1] - x2[:, :-1]), period[:-1] - np.fabs(x1[:, :-1] - x2[:, :-1]) + ) + distance = np.sqrt(np.sum(m * m, axis=len(np.shape(m)) - 1)) return distance @@ -589,7 +606,7 @@ def theta_LOS(x1, x2, period=None): x1 = np.atleast_2d(x1) x2 = np.atleast_2d(x2) if period is None: - period = np.array([np.inf]*np.shape(x1)[-1]) + period = np.array([np.inf] * np.shape(x1)[-1]) # check for consistency if np.shape(x1)[-1] != np.shape(x2)[-1]: @@ -604,9 +621,9 @@ def theta_LOS(x1, x2, period=None): # deal with zero separation r = np.sqrt(r_perp**2 + r_parallel**2) - mask = (r>0.0) + mask = r > 0.0 - theta = np.zeros(len(r)) # set to zero if r==0 - theta[mask] = np.pi/2.0 - np.arctan2(r_parallel[mask],r_perp[mask]) + theta = np.zeros(len(r)) # set to zero if r==0 + theta[mask] = np.pi / 2.0 - np.arctan2(r_parallel[mask], r_perp[mask]) return theta diff --git a/halotools/mock_observables/pair_counters/test_pair_counters/test_marked_npairs_3d.py b/halotools/mock_observables/pair_counters/test_pair_counters/test_marked_npairs_3d.py index 744f409a3..1f8eba745 100644 --- a/halotools/mock_observables/pair_counters/test_pair_counters/test_marked_npairs_3d.py +++ b/halotools/mock_observables/pair_counters/test_pair_counters/test_marked_npairs_3d.py @@ -1,17 +1,16 @@ -""" -""" +""" """ + from __future__ import absolute_import, division, print_function, unicode_literals +from pathlib import Path + import numpy as np import pytest from astropy.utils.misc import NumpyRNGContext -from pathlib import Path - -from ..pairs import wnpairs as pure_python_weighted_pairs -from ..marked_npairs_3d import marked_npairs_3d, _func_signature_int_from_wfunc from ....custom_exceptions import HalotoolsError - +from ..marked_npairs_3d import _func_signature_int_from_wfunc, marked_npairs_3d +from ..pairs import wnpairs as pure_python_weighted_pairs error_msg = ( "\nThe `test_marked_npairs_wfuncs_behavior` function performs \n" diff --git a/halotools/mock_observables/void_statistics/tests/test_underdensity_prob_func.py b/halotools/mock_observables/void_statistics/tests/test_underdensity_prob_func.py index 0cf511e8b..d463dd443 100644 --- a/halotools/mock_observables/void_statistics/tests/test_underdensity_prob_func.py +++ b/halotools/mock_observables/void_statistics/tests/test_underdensity_prob_func.py @@ -1,18 +1,18 @@ -""" Module providing unit-testing for the +"""Module providing unit-testing for the `~halotools.mock_observables.underdensity_prob_func function. """ + from __future__ import absolute_import, division, print_function import numpy as np import pytest from astropy.utils.misc import NumpyRNGContext +from ....custom_exceptions import HalotoolsError +from ...tests.cf_helpers import generate_locus_of_3d_points from ..underdensity_prob_func import underdensity_prob_func from ..void_prob_func import void_prob_func -from ...tests.cf_helpers import generate_locus_of_3d_points -from ....custom_exceptions import HalotoolsError - __all__ = ("test_upf1", "test_upf2", "test_upf3", "test_upf4") fixed_seed = 43 diff --git a/halotools/mock_observables/void_statistics/tests/test_void_prob_func.py b/halotools/mock_observables/void_statistics/tests/test_void_prob_func.py index 72d8d6245..c51acc0d2 100644 --- a/halotools/mock_observables/void_statistics/tests/test_void_prob_func.py +++ b/halotools/mock_observables/void_statistics/tests/test_void_prob_func.py @@ -1,16 +1,15 @@ -""" Module providing unit-testing for the +"""Module providing unit-testing for the `~halotools.mock_observables.void_prob_func` function. """ + from __future__ import absolute_import, division, print_function + import numpy as np import pytest from astropy.utils.misc import NumpyRNGContext -from ..void_prob_func import void_prob_func - -from ...tests.cf_helpers import generate_locus_of_3d_points - from ....custom_exceptions import HalotoolsError +from ..void_prob_func import void_prob_func __all__ = ("test_vpf1", "test_vpf2", "test_vpf3") @@ -80,22 +79,6 @@ def test_vpf3(): assert np.allclose(vpf, vpf2, rtol=0.01) -def test_vpf_process_args1(): - Npts = 1000 - Lbox = 1 - period = np.array([Lbox, Lbox, Lbox]) - with NumpyRNGContext(fixed_seed): - sample1 = np.random.random((Npts, 3)) - random_sphere_centers = np.random.random((Npts, 3)) - n_ran = 1000 - rbins = np.logspace(-1.5, -1, 5) - - with pytest.raises(HalotoolsError) as err: - __ = void_prob_func(sample1, rbins, n_ran=n_ran, period=[Lbox, Lbox]) - substr = "Input ``period`` must either be a float or length-3 sequence" - assert substr in err.value.args[0] - - def test_vpf_process_args2(): Npts = 1000 Lbox = 1 diff --git a/halotools/mock_observables/void_statistics/void_prob_func.py b/halotools/mock_observables/void_statistics/void_prob_func.py index 2fb31c942..77c642e31 100644 --- a/halotools/mock_observables/void_statistics/void_prob_func.py +++ b/halotools/mock_observables/void_statistics/void_prob_func.py @@ -6,14 +6,12 @@ from __future__ import absolute_import, division, print_function, unicode_literals import numpy as np - from astropy.utils.misc import NumpyRNGContext -from ..pair_counters import npairs_per_object_3d - -from ...utils.array_utils import array_is_monotonic from ...custom_exceptions import HalotoolsError - +from ...utils.array_utils import array_is_monotonic +from ..mock_observables_helpers import get_period +from ..pair_counters import npairs_per_object_3d __all__ = ("void_prob_func",) __author__ = ["Duncan Campbell", "Andrew Hearin"] @@ -207,6 +205,7 @@ def _void_prob_func_process_args( ) raise HalotoolsError(msg) + period = get_period(period)[0] if period is None: xmin, xmax = np.min(sample1), np.max(sample1) ymin, ymax = np.min(sample1), np.max(sample1) From c900877b5d51046d2429e56eb42ea57eafba1f35 Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 14:25:36 -0600 Subject: [PATCH 07/15] Fix bug in large_scale_density_spherical_volume.py --- .../large_scale_density_spherical_volume.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py index cc4c9683a..e54b891f5 100644 --- a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py +++ b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py @@ -149,7 +149,8 @@ def _large_scale_density_spherical_volume_process_args( if period is None: raise HalotoolsError("If sample_volume is None, must pass period") if period is not None: - raise HalotoolsError("If period is not None, do not pass in sample_volume") + if sample_volume is not None: + raise HalotoolsError("If period is not None, do not pass in sample_volume") period = get_period(period)[0] sample_volume = period.prod() From e5002847503dc681a18eae7ab250d41e6aa942d4 Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 14:33:01 -0600 Subject: [PATCH 08/15] Fix error handling in large_scale_density_spherical_annulus.py --- .../large_scale_density_spherical_annulus.py | 12 ++- .../large_scale_density_spherical_volume.py | 3 +- ...t_large_scale_density_spherical_annulus.py | 76 ++++++++++--------- 3 files changed, 52 insertions(+), 39 deletions(-) diff --git a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_annulus.py b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_annulus.py index 958694849..7dc19ec2f 100644 --- a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_annulus.py +++ b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_annulus.py @@ -153,7 +153,6 @@ def _large_scale_density_spherical_annulus_process_args( """ """ sample = np.atleast_1d(sample) tracers = np.atleast_1d(tracers) - period = get_period(period)[0] try: assert outer_radius > inner_radius @@ -162,10 +161,15 @@ def _large_scale_density_spherical_annulus_process_args( raise HalotoolsError(msg) rbins = np.array([inner_radius, outer_radius]) + if sample_volume is None: + if period is None: + raise HalotoolsError("If sample_volume is None, must pass period") + if period is not None: + if sample_volume is not None: + raise HalotoolsError("If period is not None, do not pass in sample_volume") + + period = get_period(period)[0] if sample_volume is None: sample_volume = period.prod() - else: - msg = "If period is not None, do not pass in sample_volume" - raise HalotoolsError(msg) return sample, tracers, rbins, period, sample_volume, num_threads, approx_cell1_size diff --git a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py index e54b891f5..48e5ba5c6 100644 --- a/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py +++ b/halotools/mock_observables/large_scale_density/large_scale_density_spherical_volume.py @@ -153,6 +153,7 @@ def _large_scale_density_spherical_volume_process_args( raise HalotoolsError("If period is not None, do not pass in sample_volume") period = get_period(period)[0] - sample_volume = period.prod() + if sample_volume is None: + sample_volume = period.prod() return sample, tracers, rbins, period, sample_volume, num_threads, approx_cell1_size diff --git a/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_annulus.py b/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_annulus.py index 69848f037..2cc5e1dd8 100644 --- a/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_annulus.py +++ b/halotools/mock_observables/large_scale_density/tests/test_large_scale_density_spherical_annulus.py @@ -1,81 +1,89 @@ -""" -""" -from __future__ import (absolute_import, division, print_function) +""" """ + +from __future__ import absolute_import, division, print_function + import numpy as np import pytest -from ..large_scale_density_spherical_annulus import large_scale_density_spherical_annulus - -from ...tests.cf_helpers import generate_locus_of_3d_points from ....custom_exceptions import HalotoolsError +from ...tests.cf_helpers import generate_locus_of_3d_points +from ..large_scale_density_spherical_annulus import ( + large_scale_density_spherical_annulus, +) -__all__ = ('test_large_scale_density_spherical_annulus_exception_handling', ) +__all__ = ("test_large_scale_density_spherical_annulus_exception_handling",) fixed_seed = 43 def test_large_scale_density_spherical_annulus_exception_handling(): - """ - """ + """ """ npts1, npts2 = 100, 200 sample = generate_locus_of_3d_points(npts1, xc=0.1, yc=0.1, zc=0.1, seed=fixed_seed) - tracers = generate_locus_of_3d_points(npts2, xc=0.15, yc=0.1, zc=0.1, seed=fixed_seed) + tracers = generate_locus_of_3d_points( + npts2, xc=0.15, yc=0.1, zc=0.1, seed=fixed_seed + ) inner_radius, outer_radius = 0.1, 0.2 with pytest.raises(HalotoolsError) as err: result = large_scale_density_spherical_annulus( - sample, tracers, inner_radius, outer_radius) - substr = "If period is None, you must pass in ``sample_volume``." + sample, tracers, inner_radius, outer_radius + ) + substr = "If sample_volume is None, must pass period" assert substr in err.value.args[0] with pytest.raises(HalotoolsError) as err: result = large_scale_density_spherical_annulus( - sample, tracers, inner_radius, outer_radius, period=[1, 1]) - substr = "Input ``period`` must either be a float or length-3 sequence" - assert substr in err.value.args[0] - - with pytest.raises(HalotoolsError) as err: - result = large_scale_density_spherical_annulus( - sample, tracers, inner_radius, outer_radius, period=1, sample_volume=0.4) + sample, tracers, inner_radius, outer_radius, period=1, sample_volume=0.4 + ) substr = "If period is not None, do not pass in sample_volume" assert substr in err.value.args[0] with pytest.raises(HalotoolsError) as err: result = large_scale_density_spherical_annulus( - sample, tracers, 0.5, outer_radius, period=1, sample_volume=0.4) + sample, tracers, 0.5, outer_radius, period=1, sample_volume=0.4 + ) substr = "Input ``outer_radius`` must be larger than input ``inner_radius``" assert substr in err.value.args[0] def test_large_scale_density_spherical_annulus1(): - """ - """ + """ """ npts1, npts2 = 10, 200 sample = generate_locus_of_3d_points(npts1, xc=0.1, yc=0.1, zc=0.1, seed=fixed_seed) - tracers = generate_locus_of_3d_points(npts2, xc=0.15, yc=0.1, zc=0.1, seed=fixed_seed) + tracers = generate_locus_of_3d_points( + npts2, xc=0.15, yc=0.1, zc=0.1, seed=fixed_seed + ) inner_radius, outer_radius = 0.04, 0.1 result = large_scale_density_spherical_annulus( - sample, tracers, inner_radius, outer_radius, period=1) + sample, tracers, inner_radius, outer_radius, period=1 + ) - environment_volume = (4/3.)*np.pi*(outer_radius**3 - inner_radius**3) - correct_answer = 200/environment_volume + environment_volume = (4 / 3.0) * np.pi * (outer_radius**3 - inner_radius**3) + correct_answer = 200 / environment_volume print(result) print(correct_answer) assert np.allclose(result, correct_answer, rtol=0.001) def test_large_scale_density_spherical_annulus2(): - """ - """ + """ """ npts1, npts2 = 100, 200 sample = generate_locus_of_3d_points(npts1, xc=0.1, yc=0.1, zc=0.1, seed=fixed_seed) - tracers = generate_locus_of_3d_points(npts2, xc=0.95, yc=0.1, zc=0.1, seed=fixed_seed) + tracers = generate_locus_of_3d_points( + npts2, xc=0.95, yc=0.1, zc=0.1, seed=fixed_seed + ) inner_radius, outer_radius = 0.1, 0.2 result = large_scale_density_spherical_annulus( - sample, tracers, inner_radius, outer_radius, - period=[1, 1, 1], norm_by_mean_density=True) - - environment_volume = (4/3.)*np.pi*(outer_radius**3 - inner_radius**3) + sample, + tracers, + inner_radius, + outer_radius, + period=[1, 1, 1], + norm_by_mean_density=True, + ) + + environment_volume = (4 / 3.0) * np.pi * (outer_radius**3 - inner_radius**3) mean_density = float(npts2) - correct_answer = 200/environment_volume/mean_density + correct_answer = 200 / environment_volume / mean_density assert np.allclose(result, correct_answer, rtol=0.001) From 0ec263d1c776507d9e3e51f5905801745543bf64 Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 14:35:11 -0600 Subject: [PATCH 09/15] Fix exception handling in rectangular_mesh.py --- .../mock_observables/pair_counters/rectangular_mesh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/halotools/mock_observables/pair_counters/rectangular_mesh.py b/halotools/mock_observables/pair_counters/rectangular_mesh.py index c228f2a69..e4143fd43 100644 --- a/halotools/mock_observables/pair_counters/rectangular_mesh.py +++ b/halotools/mock_observables/pair_counters/rectangular_mesh.py @@ -182,17 +182,17 @@ def __init__( try: xperiod = float(xperiod[0]) - except IndexError: + except (IndexError, TypeError): xperiod = float(xperiod) try: yperiod = float(yperiod[0]) - except IndexError: + except (IndexError, TypeError): yperiod = float(yperiod) try: zperiod = float(zperiod[0]) - except IndexError: + except (IndexError, TypeError): zperiod = float(zperiod) self.xperiod = xperiod From 2e74c5403f01f82a05cd491ce37c4f6a8f31d4fc Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 14:50:00 -0600 Subject: [PATCH 10/15] Use get_period within wnpairs function --- halotools/mock_observables/pair_counters/pairs.py | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/halotools/mock_observables/pair_counters/pairs.py b/halotools/mock_observables/pair_counters/pairs.py index 32082b8ae..eb0dd866f 100755 --- a/halotools/mock_observables/pair_counters/pairs.py +++ b/halotools/mock_observables/pair_counters/pairs.py @@ -8,6 +8,8 @@ import numpy as np +from ..mock_observables_helpers import get_period + __all__ = ["npairs", "wnpairs", "xy_z_npairs", "xy_z_wnpairs", "s_mu_npairs"] __author__ = ["Duncan Campbell"] @@ -204,15 +206,7 @@ def wnpairs(sample1, sample2, r, period=None, weights1=None, weights2=None): return None # Process period entry and check for consistency. - if period is None: - period = np.array([np.inf] * np.shape(sample1)[-1]) - else: - period = np.asarray(period).astype("float64") - if np.shape(period) == (): - period = np.array([period] * np.shape(sample1)[-1]) - if np.shape(period)[0] != np.shape(sample1)[-1]: - raise ValueError("period should have len == dimension of points") - return None + period = get_period(period)[0] # Process weights1 entry and check for consistency. if weights1 is None: @@ -245,7 +239,8 @@ def wnpairs(sample1, sample2, r, period=None, weights1=None, weights2=None): n = np.zeros((r.size,), dtype=np.float64) for i in range(r.size): for j in range(N1): - n[i] += np.sum(np.extract(dd[j, :] <= r[i], weights2)) * weights1[j] + res = np.sum(np.extract(dd[j, :] <= r[i], weights2)) * weights1[j] + n[i] += res return n From 4083f159bd2a81642784323a947d34eb48a5a18b Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 15:00:04 -0600 Subject: [PATCH 11/15] Fix bugs in test_marked_npairs_3d.py --- .../test_pair_counters/test_marked_npairs_3d.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/halotools/mock_observables/pair_counters/test_pair_counters/test_marked_npairs_3d.py b/halotools/mock_observables/pair_counters/test_pair_counters/test_marked_npairs_3d.py index 1f8eba745..ddec2ef67 100644 --- a/halotools/mock_observables/pair_counters/test_pair_counters/test_marked_npairs_3d.py +++ b/halotools/mock_observables/pair_counters/test_pair_counters/test_marked_npairs_3d.py @@ -82,8 +82,8 @@ def test_marked_npairs_3d_periodic(): random_sample, rbins, period=period, - weights1=ran_weights1, - weights2=ran_weights1, + weights1=ran_weights1.flatten(), + weights2=ran_weights1.flatten(), ) assert np.allclose(test_result, result, rtol=1e-05), "pair counts are incorrect" @@ -171,8 +171,8 @@ def test_marked_npairs_nonperiodic(): random_sample, rbins, period=None, - weights1=ran_weights1, - weights2=ran_weights1, + weights1=ran_weights1.flatten(), + weights2=ran_weights1.flatten(), ) assert np.allclose(test_result, result, rtol=1e-05), "pair counts are incorrect" From f97d4008d720780809e7f8276955bc41bdd48693 Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 15:07:14 -0600 Subject: [PATCH 12/15] Remove obsolete test --- .../test_wnpairs_pure_python.py | 41 ++++++------------- 1 file changed, 12 insertions(+), 29 deletions(-) diff --git a/halotools/mock_observables/pair_counters/test_pair_counters/test_wnpairs_pure_python.py b/halotools/mock_observables/pair_counters/test_pair_counters/test_wnpairs_pure_python.py index 94e0f2abd..96e4487f2 100644 --- a/halotools/mock_observables/pair_counters/test_pair_counters/test_wnpairs_pure_python.py +++ b/halotools/mock_observables/pair_counters/test_pair_counters/test_wnpairs_pure_python.py @@ -1,5 +1,5 @@ -""" -""" +""" """ + from __future__ import absolute_import, division, print_function, unicode_literals import numpy as np @@ -9,8 +9,7 @@ from ..pairs import wnpairs as pure_python_brute_force_wnpairs_3d from ..pairs import xy_z_wnpairs - -__all__ = ('test_wnpairs_pure_python1', ) +__all__ = ("test_wnpairs_pure_python1",) fixed_seed = 43 @@ -28,33 +27,19 @@ def test_wnpairs_pure_python1(): assert substr in err.value.args[0] -def test_wnpairs_pure_python2(): - """ - """ - npts = 10 - with NumpyRNGContext(fixed_seed): - sample1 = np.random.random((npts, 2)) - sample2 = np.random.random((npts, 2)) - rbins = np.linspace(0.01, 0.1, 5) - - with pytest.raises(ValueError) as err: - __ = pure_python_brute_force_wnpairs_3d(sample1, sample2, rbins, period=[1, 1, 1]) - substr = "period should have len == dimension of points" - assert substr in err.value.args[0] - - def test_wnpairs_pure_python3(): npts = 10 with NumpyRNGContext(fixed_seed): sample1 = np.random.random((npts, 3)) sample2 = np.random.random((npts, 3)) - weights1 = np.random.rand(npts-1) + weights1 = np.random.rand(npts - 1) weights2 = np.random.rand(npts) rbins = np.linspace(0.01, 0.1, 5) with pytest.raises(ValueError) as err: - __ = pure_python_brute_force_wnpairs_3d(sample1, sample2, rbins, - weights1=weights1, weights2=weights2) + __ = pure_python_brute_force_wnpairs_3d( + sample1, sample2, rbins, weights1=weights1, weights2=weights2 + ) substr = "weights1 should have same len as sample1" assert substr in err.value.args[0] @@ -65,12 +50,13 @@ def test_wnpairs_pure_python4(): sample1 = np.random.random((npts, 3)) sample2 = np.random.random((npts, 3)) weights1 = np.random.rand(npts) - weights2 = np.random.rand(npts-1) + weights2 = np.random.rand(npts - 1) rbins = np.linspace(0.01, 0.1, 5) with pytest.raises(ValueError) as err: - __ = pure_python_brute_force_wnpairs_3d(sample1, sample2, rbins, - weights1=weights1, weights2=weights2) + __ = pure_python_brute_force_wnpairs_3d( + sample1, sample2, rbins, weights1=weights1, weights2=weights2 + ) substr = "weights2 should have same len as sample2" assert substr in err.value.args[0] @@ -90,8 +76,7 @@ def test_xy_z_wnpairs_pure_python1(): def test_xy_z_wnpairs_pure_python2(): - """ - """ + """ """ npts = 10 with NumpyRNGContext(fixed_seed): sample1 = np.random.random((npts, 2)) @@ -103,5 +88,3 @@ def test_xy_z_wnpairs_pure_python2(): __ = xy_z_wnpairs(sample1, sample2, rp_bins, pi_bins, period=[1, 1, 1]) substr = "period should have len == dimension of points" assert substr in err.value.args[0] - - From 1af6d09913637ecbd8737c023218c926d1e0108b Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 15:15:43 -0600 Subject: [PATCH 13/15] python black --- .../tests/test_biased_isotropic_velocity.py | 21 +++++++++++-------- .../kernels/unbiased_isotropic_velocity.py | 5 ++--- 2 files changed, 14 insertions(+), 12 deletions(-) diff --git a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/tests/test_biased_isotropic_velocity.py b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/tests/test_biased_isotropic_velocity.py index 0e84e4032..e0c25ee50 100644 --- a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/tests/test_biased_isotropic_velocity.py +++ b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/tests/test_biased_isotropic_velocity.py @@ -1,13 +1,16 @@ -""" -""" +""" """ + import numpy as np from astropy.utils.data import get_pkg_data_filename -from ..unbiased_isotropic_velocity import dimensionless_radial_velocity_dispersion as unbiased_dimless_vel_rad_disp -from ..biased_isotropic_velocity import dimensionless_radial_velocity_dispersion as biased_dimless_vel_rad_disp - +from ..biased_isotropic_velocity import ( + dimensionless_radial_velocity_dispersion as biased_dimless_vel_rad_disp, +) +from ..unbiased_isotropic_velocity import ( + dimensionless_radial_velocity_dispersion as unbiased_dimless_vel_rad_disp, +) -__all__ = ('test_unbiased_vel_rad_disp1', ) +__all__ = ("test_unbiased_vel_rad_disp1",) def test_unbiased_vel_rad_disp1(): @@ -25,14 +28,14 @@ def test_unbiased_vel_rad_disp2(): halo_conc = 5 unbiased_result = unbiased_dimless_vel_rad_disp(scaled_radius, halo_conc) - gal_conc = 2*halo_conc + gal_conc = 2 * halo_conc biased_result = biased_dimless_vel_rad_disp(scaled_radius, halo_conc, gal_conc) assert not np.allclose(unbiased_result, biased_result) def test_unbiased_vel_rad_disp_external_ch10_cg5(): halo_conc, gal_conc = 10, 5 - fname = get_pkg_data_filename('data/van_den_bosch_nfw_vr_disp_ch10_cg5.dat') + fname = get_pkg_data_filename("data/van_den_bosch_nfw_vr_disp_ch10_cg5.dat") x = np.loadtxt(fname) frank_r_by_Rvir = x[:, 1] frank_dimless_sigma_rad = x[:, 2] @@ -42,7 +45,7 @@ def test_unbiased_vel_rad_disp_external_ch10_cg5(): def test_unbiased_vel_rad_disp_external_ch5_cg10(): halo_conc, gal_conc = 5, 10 - fname = get_pkg_data_filename('data/van_den_bosch_nfw_vr_disp_ch5_cg10.dat') + fname = get_pkg_data_filename("data/van_den_bosch_nfw_vr_disp_ch5_cg10.dat") x = np.loadtxt(fname) frank_r_by_Rvir = x[:, 1] frank_dimless_sigma_rad = x[:, 2] diff --git a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/unbiased_isotropic_velocity.py b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/unbiased_isotropic_velocity.py index 57ec1d134..81a76fce7 100644 --- a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/unbiased_isotropic_velocity.py +++ b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/unbiased_isotropic_velocity.py @@ -1,11 +1,10 @@ -""" -""" +""" """ + import numpy as np from scipy.integrate import quad as quad_integration from .mass_profile import _g_integral - __all__ = ("dimensionless_radial_velocity_dispersion",) From 90003a32f52c0ece04aa2924921f940255c9bb38 Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 15:37:20 -0600 Subject: [PATCH 14/15] fix most array shape bugs in nfw models --- .../analytic_models/profile_model_template.py | 10 +- .../nfw/kernels/biased_isotropic_velocity.py | 17 +- .../satellites/nfw/kernels/mass_profile.py | 4 +- .../satellites/nfw/nfw_phase_space.py | 4 +- .../satellites/nfw/nfw_profile.py | 22 +-- .../test_nfw_profile/test_nfw_profile.py | 147 +++++++++--------- 6 files changed, 107 insertions(+), 97 deletions(-) diff --git a/halotools/empirical_models/phase_space_models/analytic_models/profile_model_template.py b/halotools/empirical_models/phase_space_models/analytic_models/profile_model_template.py index 04002e928..b53a735ad 100644 --- a/halotools/empirical_models/phase_space_models/analytic_models/profile_model_template.py +++ b/halotools/empirical_models/phase_space_models/analytic_models/profile_model_template.py @@ -3,17 +3,17 @@ a container class for the distribution of mass and/or galaxies within dark matter halos. """ -from __future__ import division, print_function, absolute_import, unicode_literals + +from __future__ import absolute_import, division, print_function, unicode_literals import numpy as np -from scipy.integrate import quad as quad_integration -from scipy.optimize import minimize as scipy_minimize from astropy import units as u from astropy.constants import G - -from . import halo_boundary_functions +from scipy.integrate import quad as quad_integration +from scipy.optimize import minimize as scipy_minimize from ... import model_defaults +from . import halo_boundary_functions newtonG = G.to(u.km * u.km * u.Mpc / (u.Msun * u.s * u.s)) diff --git a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/biased_isotropic_velocity.py b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/biased_isotropic_velocity.py index ba31219f5..1bc6dd5b5 100644 --- a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/biased_isotropic_velocity.py +++ b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/biased_isotropic_velocity.py @@ -1,11 +1,10 @@ -""" -""" +""" """ + import numpy as np from scipy.integrate import quad as quad_integration from .mass_profile import _g_integral - __all__ = ("dimensionless_radial_velocity_dispersion",) @@ -58,21 +57,25 @@ def dimensionless_radial_velocity_dispersion( gal_conc * gal_conc * x * (1.0 + gal_conc * x) ** 2 / _g_integral(halo_conc) ) extra_args = halo_conc / np.atleast_1d(gal_conc).astype("f4") + assert extra_args.shape == (1,) + extra_args = float(extra_args[0]) lower_limit = gal_conc * x upper_limit = float("inf") for i in range(len(x)): + a = lower_limit[i] + b = upper_limit term1, _ = quad_integration( _jeans_integrand_term1, - lower_limit[i], - upper_limit, + a, + b, epsrel=profile_integration_tol, args=extra_args, ) term2, _ = quad_integration( _jeans_integrand_term2, - lower_limit[i], - upper_limit, + a, + b, epsrel=profile_integration_tol, args=extra_args, ) diff --git a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/mass_profile.py b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/mass_profile.py index a4ff427cb..013396523 100644 --- a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/mass_profile.py +++ b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/kernels/mass_profile.py @@ -1,5 +1,5 @@ -""" -""" +""" """ + import numpy as np __all__ = ("cumulative_mass_PDF", "dimensionless_mass_density") diff --git a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/nfw_phase_space.py b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/nfw_phase_space.py index 3602eb887..cbfdc927a 100644 --- a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/nfw_phase_space.py +++ b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/nfw_phase_space.py @@ -1,5 +1,5 @@ -""" -""" +""" """ + from __future__ import absolute_import, division, print_function import numpy as np diff --git a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/nfw_profile.py b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/nfw_profile.py index 4b442fef7..97572e145 100644 --- a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/nfw_profile.py +++ b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/nfw_profile.py @@ -5,20 +5,20 @@ Navarry, Frenk and White (1995), `arXiv:9508025 `_. a sub-class of `~halotools.empirical_models.AnalyticDensityProf`. """ -from __future__ import division, print_function, absolute_import, unicode_literals -import numpy as np - -from .conc_mass import direct_from_halo_catalog, dutton_maccio14 -from .kernels import nfw_dimensionless_mass_density, nfw_cumulative_mass_PDF -from .kernels import standalone_mc_generate_nfw_radial_positions +from __future__ import absolute_import, division, print_function, unicode_literals -from ...profile_model_template import AnalyticDensityProf - -from ..... import model_defaults +import numpy as np from ......sim_manager import sim_defaults - +from ..... import model_defaults +from ...profile_model_template import AnalyticDensityProf +from .conc_mass import direct_from_halo_catalog, dutton_maccio14 +from .kernels import ( + nfw_cumulative_mass_PDF, + nfw_dimensionless_mass_density, + standalone_mc_generate_nfw_radial_positions, +) __all__ = ("NFWProfile",) __author__ = ("Andrew Hearin", "Benedikt Diemer") @@ -44,7 +44,7 @@ def __init__( conc_mass_model=model_defaults.conc_mass_model, concentration_key=model_defaults.concentration_key, halo_boundary_key=None, - **kwargs + **kwargs, ): r""" Parameters diff --git a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/tests/test_nfw_profile/test_nfw_profile.py b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/tests/test_nfw_profile/test_nfw_profile.py index baa32fc2c..bb3f3d3d9 100644 --- a/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/tests/test_nfw_profile/test_nfw_profile.py +++ b/halotools/empirical_models/phase_space_models/analytic_models/satellites/nfw/tests/test_nfw_profile/test_nfw_profile.py @@ -1,53 +1,53 @@ -""" Module providing unit-testing for `~halotools.empirical_models.NFWProfile` class -""" +"""Module providing unit-testing for `~halotools.empirical_models.NFWProfile` class""" + from __future__ import absolute_import, division, print_function, unicode_literals import numpy as np - import pytest - -from astropy.cosmology import WMAP9, FLRW - -from ...nfw_profile import NFWProfile +from astropy.cosmology import FLRW, WMAP9 from ........utils.array_utils import array_is_monotonic +from ...nfw_profile import NFWProfile -__all__ = ('test_instance_attrs', 'analytic_nfw_density_outer_shell_normalization', - 'monte_carlo_density_outer_shell_normalization') +__all__ = ( + "test_instance_attrs", + "analytic_nfw_density_outer_shell_normalization", + "monte_carlo_density_outer_shell_normalization", +) def test_attr_inheritance(): - r""" Test that `~halotools.empirical_models.NFWProfile` + r"""Test that `~halotools.empirical_models.NFWProfile` possesses the necessary attributes and methods. """ - model_instance = NFWProfile(cosmology=WMAP9, redshift=2, mdef='vir') + model_instance = NFWProfile(cosmology=WMAP9, redshift=2, mdef="vir") - assert hasattr(model_instance, 'cosmology') + assert hasattr(model_instance, "cosmology") assert isinstance(model_instance.cosmology, FLRW) - assert hasattr(model_instance, 'redshift') + assert hasattr(model_instance, "redshift") assert model_instance.redshift == 2 - assert hasattr(model_instance, 'mdef') - assert model_instance.mdef == 'vir' + assert hasattr(model_instance, "mdef") + assert model_instance.mdef == "vir" - assert hasattr(model_instance, 'halo_boundary_key') - assert model_instance.halo_boundary_key == 'halo_rvir' + assert hasattr(model_instance, "halo_boundary_key") + assert model_instance.halo_boundary_key == "halo_rvir" - assert hasattr(model_instance, 'prim_haloprop_key') - assert model_instance.prim_haloprop_key == 'halo_mvir' + assert hasattr(model_instance, "prim_haloprop_key") + assert model_instance.prim_haloprop_key == "halo_mvir" - assert hasattr(model_instance, 'param_dict') - assert hasattr(model_instance, 'publications') - assert hasattr(model_instance, 'halo_prof_param_keys') + assert hasattr(model_instance, "param_dict") + assert hasattr(model_instance, "publications") + assert hasattr(model_instance, "halo_prof_param_keys") - assert hasattr(model_instance, 'virial_velocity') + assert hasattr(model_instance, "virial_velocity") vvir = model_instance.virial_velocity(total_mass=1e12) def analytic_nfw_density_outer_shell_normalization(radii, conc): - r""" Density of an NFW profile normalized by the density evaluated at the outermost value of the input ``radii`` array. + r"""Density of an NFW profile normalized by the density evaluated at the outermost value of the input ``radii`` array. For an NFW profile we have the following analytical relation: @@ -73,13 +73,13 @@ def analytic_nfw_density_outer_shell_normalization(radii, conc): """ outer_radius = radii[-1] - numerator = outer_radius*(1 + conc*outer_radius)**2 - denominator = radii*(1 + conc*radii)**2 - return numerator/denominator + numerator = outer_radius * (1 + conc * outer_radius) ** 2 + denominator = radii * (1 + conc * radii) ** 2 + return numerator / denominator def monte_carlo_density_outer_shell_normalization(rbins, radial_positions): - r""" Density of a Monte Carlo realization of a spherically symmetric profile normalized by the density evaluated at the midpoint of the outermost bin of the input ``rbins`` array. + r"""Density of a Monte Carlo realization of a spherically symmetric profile normalized by the density evaluated at the midpoint of the outermost bin of the input ``rbins`` array. Parameters ------------ @@ -98,35 +98,36 @@ def monte_carlo_density_outer_shell_normalization(rbins, radial_positions): Ratio of the profile number density scaled by the number density at the midpoint of the outermost bin of the input ``rbins`` array. """ - rbin_midpoints = 0.5*(rbins[:-1] + rbins[1:]) + rbin_midpoints = 0.5 * (rbins[:-1] + rbins[1:]) counts = np.histogram(radial_positions, bins=rbins)[0].astype(np.float64) outer_radius = rbin_midpoints[-1] outer_counts = counts[-1] - return rbin_midpoints, (counts/rbin_midpoints**2)/(outer_counts/outer_radius**2) + return rbin_midpoints, (counts / rbin_midpoints**2) / ( + outer_counts / outer_radius**2 + ) def test_instance_attrs(): - r""" Require that all model variants have ``cosmology``, ``redshift`` and ``mdef`` attributes. - """ + r"""Require that all model variants have ``cosmology``, ``redshift`` and ``mdef`` attributes.""" default_nfw = NFWProfile(concentration_bins=np.array((5, 10, 15))) wmap9_nfw = NFWProfile(cosmology=WMAP9, concentration_bins=np.array((5, 10, 15))) - m200_nfw = NFWProfile(mdef='200m', concentration_bins=np.array((5, 10, 15))) + m200_nfw = NFWProfile(mdef="200m", concentration_bins=np.array((5, 10, 15))) - assert hasattr(default_nfw, 'cosmology') - assert hasattr(wmap9_nfw, 'cosmology') - assert hasattr(m200_nfw, 'cosmology') + assert hasattr(default_nfw, "cosmology") + assert hasattr(wmap9_nfw, "cosmology") + assert hasattr(m200_nfw, "cosmology") - assert hasattr(default_nfw, 'redshift') - assert hasattr(wmap9_nfw, 'redshift') - assert hasattr(m200_nfw, 'redshift') + assert hasattr(default_nfw, "redshift") + assert hasattr(wmap9_nfw, "redshift") + assert hasattr(m200_nfw, "redshift") - assert hasattr(default_nfw, 'mdef') - assert hasattr(wmap9_nfw, 'mdef') - assert hasattr(m200_nfw, 'mdef') + assert hasattr(default_nfw, "mdef") + assert hasattr(wmap9_nfw, "mdef") + assert hasattr(m200_nfw, "mdef") def test_mass_density(): - r""" Require the returned value of the + r"""Require the returned value of the `~halotools.empirical_models.NFWProfile.mass_density` function to be self-consistent with the returned value of the `~halotools.empirical_models.NFWProfile.dimensionless_mass_density` function. @@ -138,23 +139,22 @@ def test_mass_density(): default_nfw = NFWProfile(concentration_bins=np.array((5, 10, 15))) wmap9_nfw = NFWProfile(cosmology=WMAP9, concentration_bins=np.array((5, 10, 15))) - m200_nfw = NFWProfile(mdef='200m', concentration_bins=np.array((5, 10, 15))) + m200_nfw = NFWProfile(mdef="200m", concentration_bins=np.array((5, 10, 15))) model_list = [default_nfw, wmap9_nfw, m200_nfw] for model in model_list: result = model.mass_density(radius, mass, conc) halo_radius = model.halo_mass_to_halo_radius(mass) - scaled_radius = radius/halo_radius - derived_result = ( - model.density_threshold * - model.dimensionless_mass_density(scaled_radius, conc) - ) + scaled_radius = radius / halo_radius + derived_result = model.density_threshold * model.dimensionless_mass_density( + scaled_radius, conc + ) assert np.allclose(derived_result, result, rtol=1e-4) def test_cumulative_mass_PDF(): - r""" Require the `~halotools.empirical_models.NFWProfile.cumulative_mass_PDF` + r"""Require the `~halotools.empirical_models.NFWProfile.cumulative_mass_PDF` method in all model variants to respect a number of consistency conditions. 1. Returned value is a strictly monotonically increasing array between 0 and 1. @@ -185,10 +185,10 @@ def test_cumulative_mass_PDF(): Npts = 100 total_mass = np.zeros(Npts) + 1e12 scaled_radius = np.logspace(-2, -0.01, Npts) - conc = 5 + conc = 5.0 default_nfw = NFWProfile(concentration_bins=np.array((5, 10, 15))) - m200_nfw = NFWProfile(mdef='200m', concentration_bins=np.array((5, 10, 15))) + m200_nfw = NFWProfile(mdef="200m", concentration_bins=np.array((5, 10, 15))) model_list = [default_nfw, m200_nfw] for model in model_list: @@ -203,20 +203,21 @@ def test_cumulative_mass_PDF(): # and the direct numerical integral of the analytical expression for # dimensionless_mass_density super_class_result = super(NFWProfile, model).cumulative_mass_PDF( - scaled_radius, conc) + scaled_radius, conc + ) assert np.allclose(super_class_result, result, rtol=1e-4) # Verify that we get a self-consistent result between # enclosed_mass and cumulative_mass_PDF halo_radius = model.halo_mass_to_halo_radius(total_mass) - radius = scaled_radius*halo_radius + radius = scaled_radius * halo_radius enclosed_mass = model.enclosed_mass(radius, total_mass, conc) - derived_enclosed_mass = result*total_mass + derived_enclosed_mass = result * total_mass assert np.allclose(enclosed_mass, derived_enclosed_mass, rtol=1e-4) def test_vmax(): - r""" Require that the analytic approximation used to estimate the NFW :math:`V_{\rm max}` + r"""Require that the analytic approximation used to estimate the NFW :math:`V_{\rm max}` by the `~halotools.empirical_models.NFWProfile.vmax` method agrees with the maximum value of :math:`V_{\rm circ}(r)` computed over the entire profile of the halo computed using the super-class method @@ -228,7 +229,7 @@ def test_vmax(): radius_array = np.logspace(-2, 0, npts) default_nfw = NFWProfile() - m200_nfw = NFWProfile(mdef='200m') + m200_nfw = NFWProfile(mdef="200m") model_list = [default_nfw, m200_nfw] for model in model_list: @@ -240,7 +241,7 @@ def test_vmax(): @pytest.mark.slow def test_mc_generate_nfw_radial_positions(): - r""" Require that the points returned by the + r"""Require that the points returned by the `~halotools.empirical_models.NFWProfile.mc_generate_nfw_radial_positions` function do indeed trace an NFW profile. @@ -280,15 +281,18 @@ def test_mc_generate_nfw_radial_positions(): for conc in conc_to_test: radial_positions = default_nfw.mc_generate_nfw_radial_positions( - halo_radius=halo_radius, conc=conc, num_pts=num_pts, seed=43) + halo_radius=halo_radius, conc=conc, num_pts=num_pts, seed=43 + ) radial_positions /= halo_radius rbin_midpoints, monte_carlo_ratio = ( - monte_carlo_density_outer_shell_normalization(rbins, radial_positions)) + monte_carlo_density_outer_shell_normalization(rbins, radial_positions) + ) - analytical_ratio = ( - analytic_nfw_density_outer_shell_normalization(rbin_midpoints, conc)) + analytical_ratio = analytic_nfw_density_outer_shell_normalization( + rbin_midpoints, conc + ) assert np.allclose(monte_carlo_ratio, analytical_ratio, 0.05) @@ -303,29 +307,32 @@ def test_mc_generate_nfw_radial_positions_stochasticity(): model = NFWProfile() r1 = model.mc_generate_nfw_radial_positions( - halo_radius=halo_radius, conc=conc, num_pts=num_pts, seed=43) + halo_radius=halo_radius, conc=conc, num_pts=num_pts, seed=43 + ) r2 = model.mc_generate_nfw_radial_positions( - halo_radius=halo_radius, conc=conc, num_pts=num_pts, seed=43) + halo_radius=halo_radius, conc=conc, num_pts=num_pts, seed=43 + ) r3 = model.mc_generate_nfw_radial_positions( - halo_radius=halo_radius, conc=conc, num_pts=num_pts, seed=44) + halo_radius=halo_radius, conc=conc, num_pts=num_pts, seed=44 + ) assert np.allclose(r1, r2, rtol=0.001) assert not np.allclose(r1, r3, rtol=0.001) def test_user_defined_halo_radius(): - """Regression test for #940 - """ + """Regression test for #940""" M = 1e14 model = NFWProfile() virial_radius = model.halo_mass_to_halo_radius(M) user_defined_radius = 7 r = model.mc_generate_nfw_radial_positions( - halo_radius=user_defined_radius, conc=5, num_pts=int(1000), seed=43) + halo_radius=user_defined_radius, conc=5, num_pts=int(1000), seed=43 + ) assert np.all(r <= user_defined_radius) assert np.any(r > virial_radius) r = model.mc_generate_nfw_radial_positions( - halo_mass=M, conc=5, num_pts=int(1000), seed=43) + halo_mass=M, conc=5, num_pts=int(1000), seed=43 + ) assert np.all(r <= virial_radius) - From 124aa8f83813c4a9fc567f73884ddaa9fe07838b Mon Sep 17 00:00:00 2001 From: Andrew Hearin Date: Mon, 26 Jan 2026 15:46:09 -0600 Subject: [PATCH 15/15] Fix bug in NFW cumulative_mass_PDF --- .../analytic_models/profile_model_template.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/halotools/empirical_models/phase_space_models/analytic_models/profile_model_template.py b/halotools/empirical_models/phase_space_models/analytic_models/profile_model_template.py index b53a735ad..c732ece33 100644 --- a/halotools/empirical_models/phase_space_models/analytic_models/profile_model_template.py +++ b/halotools/empirical_models/phase_space_models/analytic_models/profile_model_template.py @@ -223,9 +223,12 @@ def cumulative_mass_PDF(self, scaled_radius, *prof_params): x = np.atleast_1d(scaled_radius).astype(np.float64) enclosed_mass = np.zeros_like(x) + def _scalar_integrand(*args, **kwargs): + return self._enclosed_dimensionless_mass_integrand(*args, **kwargs)[0] + for i in range(len(x)): enclosed_mass[i], _ = quad_integration( - self._enclosed_dimensionless_mass_integrand, + _scalar_integrand, 0.0, x[i], epsrel=1e-5, @@ -233,7 +236,7 @@ def cumulative_mass_PDF(self, scaled_radius, *prof_params): ) total, _ = quad_integration( - self._enclosed_dimensionless_mass_integrand, + _scalar_integrand, 0.0, 1.0, epsrel=1e-5,