Coverage for src/galsbi/sham_ot/histograms.py: 100%
227 statements
« prev ^ index » next coverage.py v7.17.2, created at 2026-11-09 18:05 +0000
1# Copyright (C) 2025 ETH Zurich
2# Institute for Particle Physics and Astrophysics
3# Author: Silvan Fischbacher
4# created: Tue Jun 24 2025
7import pickle
8import re
10import numpy as np
11from cosmic_toolbox.logger import get_logger
12from scipy.special import erf
14from galsbi.sham_ot import utils
16LOGGER = get_logger(__file__)
19class HaloHistogram:
20 """
21 Class to handle histograms of halos for optimal use in the SHAM-OT framework.
22 The class saves the following histograms:
24 - 2D histogram of redshift vs halo mass for all halos,
25 shape (n_z_bins, n_mass_bins)
26 - 2D histogram of redshift vs halo mass for central halos,
27 shape (n_z_bins, n_mass_bins)
28 - 3D histogram of redshift vs halo mass vs time since last merger for subhalos,
29 shape (n_z_bins, n_mass_bins, n_time_bins)
31 The indexing follows: hist[z_bin_index, mass_bin_index, time_bin_index],
32 i.e. histogram[0] gives the 0th redshift bin as a histogram over log masses.
33 """
35 def __init__(
36 self,
37 z_min=None,
38 z_max=None,
39 min_mass=None,
40 max_mass=6e15,
41 min_time_since_merger=0,
42 max_time_since_merger=6,
43 n_bins_mass=500,
44 n_bins_z=500,
45 n_bins_time=100,
46 ):
47 """
48 Initializes the HaloHistogram class with specified mass and redshift ranges
49 and the number of bins for each.
50 Defaults for redshift and minimum mass are taken directly from the simulation
51 if not specified.
53 :param z_min: Minimum redshift of the histogram
54 :param z_max: Maximum redshift of the histogram
55 :param min_mass: Minimum mass of the histogram
56 :param max_mass: Maximum mass of the histogram (default: 6e15)
57 :param min_time_since_merger: Minimum time since last merger (default: 0)
58 :param max_time_since_merger: Maximum time since last merger (default: 10)
59 :param n_bins_mass: Number of bins in the mass dimension (default: 500)
60 :param n_bins_z: Number of bins in the redshift dimension (default: 500)
61 :param n_bins_time: Number of bins in the time since last merger dimension
62 (default: 100)
63 """
64 self.z_min = z_min
65 self.z_max = z_max
66 self.min_mass = min_mass
67 self.max_mass = max_mass
68 self.min_time_since_merger = min_time_since_merger
69 self.max_time_since_merger = max_time_since_merger
70 self.n_bins_mass = n_bins_mass
71 self.n_bins_z = n_bins_z
72 self.n_bins_time = n_bins_time
74 # Flag to enforce initialization of bins and histograms the first time
75 # a halo catalog is added
76 self.no_catalogs_added = True
78 @property
79 def min_log_mass(self):
80 """
81 Returns the minimum log mass of the histogram.
82 If min_mass is not set, it returns None.
83 """
84 if self.min_mass is not None:
85 return np.log10(self.min_mass)
86 return None
88 @property
89 def max_log_mass(self):
90 """
91 Returns the maximum log mass of the histogram.
92 If max_mass is not set, it returns None.
93 """
94 return np.log10(self.max_mass)
96 @property
97 def log_mass_bins(self):
98 """
99 Returns the logarithmic mass bins for the histogram.
100 If min_log_mass or max_log_mass is not set, it returns None.
101 """
102 if hasattr(self, "mass_bins"):
103 return np.log10(self.mass_bins)
104 return None
106 def add_halo_catalog(
107 self,
108 halo_catalog,
109 mass_field="m",
110 redshift_field="z_true",
111 is_central_field="is_halo",
112 time_since_merger_field="delta_t",
113 pinocchio_output_filename=None,
114 ):
115 """
116 Adds a halo catalog to the histograms. The catalog should be a recarray or
117 dictionary-like object containing the necessary fields:
118 - Halo mass
119 - Redshift
120 - Is central halo
121 - Time since last merger
123 :param halo_catalog: A recarray or dictionary-like object containing halo data.
124 :param mass_field: Field name for halo mass
125 (default: 'm')
126 :param redshift_field: Field name for redshift
127 (default: 'z_true')
128 :param is_central_field: Field name for central halo indicator
129 (default: 'is_halo')
130 :param time_since_merger_field: Field name for time since last merger
131 (default: 'delta_t')
132 """
134 if self.no_catalogs_added:
135 if self.z_min is None or self.z_max is None or self.min_mass is None:
136 m_part, min_num_part, z_min, z_max = read_pinocchio_config_details(
137 pinocchio_output_filename
138 )
139 if self.z_min is None:
140 self.z_min = z_min
141 if self.z_max is None:
142 self.z_max = z_max
143 if self.min_mass is None:
144 self.min_mass = 0.999 * m_part * min_num_part
146 log_mass_bins = np.linspace(
147 self.min_log_mass, self.max_log_mass, self.n_bins_mass + 1
148 )
149 self.mass_bins = np.power(10, log_mass_bins)
150 self.z_bins = np.linspace(self.z_min, self.z_max, self.n_bins_z + 1)
151 self.time_bins = np.linspace(
152 self.min_time_since_merger,
153 self.max_time_since_merger,
154 self.n_bins_time + 1,
155 )
157 self.histogram_all = np.zeros((self.n_bins_z, self.n_bins_mass), dtype=int)
158 self.histogram_central = np.zeros(
159 (self.n_bins_z, self.n_bins_mass), dtype=int
160 )
161 self.histogram_subhalo = np.zeros(
162 (self.n_bins_z, self.n_bins_mass, self.n_bins_time), dtype=int
163 )
164 self.no_catalogs_added = False
166 is_central = halo_catalog[is_central_field].astype(bool)
168 self.histogram_all += np.histogram2d(
169 halo_catalog[redshift_field],
170 np.log10(halo_catalog[mass_field]),
171 bins=(self.z_bins, self.log_mass_bins),
172 )[0].astype(int)
174 self.histogram_central += np.histogram2d(
175 halo_catalog[redshift_field][is_central],
176 np.log10(halo_catalog[mass_field][is_central]),
177 bins=(self.z_bins, self.log_mass_bins),
178 )[0].astype(int)
180 self.histogram_subhalo += np.histogramdd(
181 (
182 halo_catalog[redshift_field][~is_central],
183 np.log10(halo_catalog[mass_field][~is_central]),
184 halo_catalog[time_since_merger_field][~is_central],
185 ),
186 bins=(self.z_bins, self.log_mass_bins, self.time_bins),
187 )[0].astype(int)
189 if np.any(halo_catalog[mass_field] > self.max_mass):
190 too_high = halo_catalog[mass_field] > self.max_mass
191 LOGGER.warning(
192 f"{sum(too_high)} halos have mass > {self.max_mass}. "
193 "They are added to the heaviest mass bin"
194 )
195 # Add overflow halos to the correct redshift bins in the last mass bin
196 overflow_hist_all = np.histogram(
197 halo_catalog[redshift_field][too_high], bins=self.z_bins
198 )[0].astype(int)
199 overflow_hist_central = np.histogram(
200 halo_catalog[redshift_field][too_high & is_central], bins=self.z_bins
201 )[0].astype(int)
202 overflow_hist_subhalo = np.histogramdd(
203 (
204 halo_catalog[redshift_field][too_high & ~is_central],
205 halo_catalog[time_since_merger_field][too_high & ~is_central],
206 ),
207 bins=(self.z_bins, self.time_bins),
208 )[0].astype(int)
210 self.histogram_all[:, -1] += overflow_hist_all
211 self.histogram_central[:, -1] += overflow_hist_central
212 self.histogram_subhalo[:, -1, :] += overflow_hist_subhalo
214 if np.any(halo_catalog[time_since_merger_field] > self.max_time_since_merger):
215 too_late = (
216 halo_catalog[time_since_merger_field] > self.max_time_since_merger
217 )
218 LOGGER.warning(
219 f"{sum(too_late)} subhalos have time since last merger > "
220 "{self.max_time_since_merger}. They are added to the latest time bin."
221 )
222 subhalos_too_late = too_late & ~is_central
223 if np.any(subhalos_too_late):
224 overflow_hist_subhalo = np.histogram2d(
225 halo_catalog[redshift_field][subhalos_too_late],
226 np.log10(halo_catalog[mass_field][subhalos_too_late]),
227 bins=(self.z_bins, self.log_mass_bins),
228 )[0].astype(int)
229 self.histogram_subhalo[:, :, -1] += overflow_hist_subhalo
231 def add_halo_histogram(self, halo_histogram):
232 """
233 Adds another HaloHistogram instance to the current one.
235 :param halo_histogram: Another HaloHistogram instance to add.
236 """
237 if not isinstance(halo_histogram, HaloHistogram):
238 raise TypeError("halo_histogram must be an instance of HaloHistogram")
240 if (
241 not np.array_equal(self.mass_bins, halo_histogram.mass_bins)
242 or not np.array_equal(self.z_bins, halo_histogram.z_bins)
243 or not np.array_equal(self.time_bins, halo_histogram.time_bins)
244 ):
245 raise ValueError("Mass, redshift, or time bins do not match.")
247 self.histogram_all += halo_histogram.histogram_all
248 self.histogram_central += halo_histogram.histogram_central
249 self.histogram_subhalo += halo_histogram.histogram_subhalo
251 def get_histogram(self, galaxy_type="blue", model="Berner+24", **kwargs):
252 """
253 Returns the histogram for red or blue galaxies for the specified SHAM model.
254 Currently, the following models are supported:
256 - "Berner+24": free parameters m_limit and quenching time, see arXiv:2310.15223
257 - "Peng+10": quenching fractions based on mass, see arXiv:1003.4747
259 :param galaxy_type: Type of galaxy ("blue" or "red")
260 :param model: Model name (default: "Berner+24")
261 :param m_limit: Mass limit for central halos in Berner+24 (default: 9e12)
262 :param quenching_time: Quenching time for subhalos in Berner+24 (default: 2)
263 :param mlim_sigma: Scatter in mass limit for extended Berner+24 (default: 0.3)
264 :param initial_quenched_fraction: Initial quenched fraction for subhalos in
265 extended Berner+24 (default: 0.1)
266 :param quenching_time_sigma: Scatter in quenching time for extended Berner+24
267 (default: 0.5)
268 :param M_mu: Characteristic mass for Peng+10 (default: 11.935)
269 :param sigma: Mass scatter for Peng+10 (default: 0.51)
270 :param b: Satellite quenching efficiency for Peng+10 (default: 0.41)
271 :return: Histogram for the specified galaxy type
272 """
273 if model == "Berner+24":
274 kwargs.setdefault("m_limit", 8e12)
275 kwargs.setdefault("quenching_time", 2.0)
276 return self._get_histogram_berner(
277 galaxy_type, kwargs["m_limit"], kwargs["quenching_time"]
278 )
279 elif model == "Berner+24_ext":
280 kwargs.setdefault("m_limit", 8e12)
281 kwargs.setdefault("quenching_time", 2.0)
282 kwargs.setdefault("mlim_sigma", 0.3)
283 kwargs.setdefault("initial_quenched_fraction", 0.1)
284 kwargs.setdefault("quenching_time_sigma", 0.5)
285 return self._get_histogram_berner_ext(
286 galaxy_type,
287 kwargs["m_limit"],
288 kwargs["quenching_time"],
289 mlim_sigma=kwargs["mlim_sigma"],
290 initial_quenched_fraction=kwargs["initial_quenched_fraction"],
291 quenching_time_sigma=kwargs["quenching_time_sigma"],
292 )
293 elif model == "Peng+10":
294 kwargs.setdefault("M_mu", 11.935)
295 kwargs.setdefault("sigma", 0.51)
296 kwargs.setdefault("b", 0.41)
297 return self._get_histogram_peng(
298 galaxy_type, kwargs["M_mu"], kwargs["sigma"], kwargs["b"]
299 )
300 else:
301 raise ValueError(
302 f"Model '{model}' is not supported. Currently only the following "
303 "models are supported: 'Berner+24', 'Berner+24_ext, 'Peng+10'."
304 )
306 def _get_histogram_berner(self, galaxy_type, m_limit, quenching_time):
307 """
308 Get histogram using Berner+24 model with hard cuts.
309 """
310 histogram = np.zeros_like(self.histogram_all)
312 # Create fractional masks for mass and time bins
313 mass_bin_centers = self.get_mass_bins()[1]
314 time_bin_centers = self.get_time_bins()[1]
316 # For mass bins: fraction of bin that is <= m_limit
317 mass_blue_mask = np.zeros(len(mass_bin_centers))
318 for i, (bin_start, bin_end) in enumerate(
319 zip(self.mass_bins[:-1], self.mass_bins[1:], strict=True)
320 ):
321 if bin_end <= m_limit:
322 mass_blue_mask[i] = 1.0 # Entire bin is blue
323 elif bin_start > m_limit:
324 mass_blue_mask[i] = 0.0 # Entire bin is red
325 else:
326 # Limit is within bin, calculate fraction
327 mass_blue_mask[i] = (m_limit - bin_start) / (bin_end - bin_start)
329 # For time bins: fraction of bin that is <= quenching_time
330 time_blue_mask = np.zeros(len(time_bin_centers))
331 for i, (bin_start, bin_end) in enumerate(
332 zip(self.time_bins[:-1], self.time_bins[1:], strict=True)
333 ):
334 if bin_end <= quenching_time:
335 time_blue_mask[i] = 1.0 # Entire bin is blue
336 elif bin_start > quenching_time:
337 time_blue_mask[i] = 0.0 # Entire bin is red
338 else:
339 # Limit is within bin, calculate fraction
340 time_blue_mask[i] = (quenching_time - bin_start) / (bin_end - bin_start)
342 if galaxy_type == "blue":
343 # Blue galaxies:
344 # central halos with mass <= m_limit
345 # subhalos with delta_t <= quenching_time
346 histogram += np.round(
347 self.histogram_central * mass_blue_mask[None, :]
348 ).astype(int)
349 histogram += (
350 np.round(self.histogram_subhalo * time_blue_mask[None, None, :])
351 .sum(axis=2)
352 .astype(int)
353 )
354 elif galaxy_type == "red":
355 # Red galaxies:
356 # central halos with mass > m_limit
357 # subhalos with delta_t > quenching_time
358 histogram += np.round(
359 self.histogram_central * (1 - mass_blue_mask)[None, :]
360 ).astype(int)
361 histogram += (
362 np.round(self.histogram_subhalo * (1 - time_blue_mask)[None, None, :])
363 .sum(axis=2)
364 .astype(int)
365 )
366 else:
367 raise ValueError(
368 f"Galaxy type '{galaxy_type}' not supported. Use 'blue' or 'red'."
369 )
371 return histogram
373 def _get_histogram_berner_ext(
374 self,
375 galaxy_type,
376 m_limit,
377 quenching_time,
378 mlim_sigma,
379 initial_quenched_fraction,
380 quenching_time_sigma,
381 ):
382 """
383 Get histogram using extended Berner+24 model with mass-dependent quenching
384 fractions.
385 """
387 histogram = np.zeros_like(self.histogram_all, dtype=float)
389 mass_centers = self.get_mass_bins()[1]
390 time_centers = self.get_time_bins()[1]
392 red_centrals = utils.adapted_err_function(
393 np.log10(mass_centers), np.log10(m_limit), mlim_sigma, min_val=0, max_val=1
394 )
395 blue_centrals = 1 - red_centrals
397 red_subhalos = utils.adapted_err_function(
398 time_centers,
399 quenching_time,
400 quenching_time_sigma,
401 min_val=initial_quenched_fraction,
402 max_val=1,
403 )
404 blue_subhalos = 1 - red_subhalos
406 histogram_central = self.histogram_central.astype(float)
407 histogram_subhalo = self.histogram_subhalo.astype(float)
409 if galaxy_type == "blue":
410 # Blue galaxies = (1 - quenched fraction) * histogram
411 histogram += histogram_central * blue_centrals[None, :]
412 # subhalo histogram is 3D, the selection is over the 3rd dimension
413 histogram += (histogram_subhalo * blue_subhalos[None, None, :]).sum(axis=2)
414 elif galaxy_type == "red":
415 # Red galaxies = quenched fraction * histogram
416 histogram += histogram_central * red_centrals[None, :]
417 # subhalo histogram is 3D, the selection is over the 3rd dimension
418 histogram += (histogram_subhalo * red_subhalos[None, None, :]).sum(axis=2)
419 else:
420 raise ValueError(
421 f"Galaxy type '{galaxy_type}' not supported. Use 'blue' or 'red'."
422 )
424 return histogram
426 def _get_histogram_peng(self, galaxy_type, M_mu, sigma, b):
427 """
428 Get histogram using Peng+10 model with mass-dependent quenching fractions.
429 Attention: M_mu is in log10 mass units.
430 """
432 histogram = np.zeros_like(self.histogram_all, dtype=float)
434 mass_centers = self.get_mass_bins()[1]
436 def _r0(M_mu, sigma, masses):
437 return np.log10(masses / (10**M_mu)) / (np.sqrt(2) * sigma)
439 def _quenched_fraction_central(M_mu, sigma, masses):
440 return (1 / 2) * (1 + erf(_r0(M_mu, sigma, masses)))
442 def _quenched_fraction_subhalo(M_mu, sigma, b, masses):
443 return ((1 - b) / 2) * (1 + erf(_r0(M_mu, sigma, masses))) + b
445 # Calculate quenching fractions for each mass bin
446 quenched_frac_central = _quenched_fraction_central(M_mu, sigma, mass_centers)
447 quenched_frac_subhalo = 1 - _quenched_fraction_subhalo(
448 M_mu, sigma, b, mass_centers
449 )
451 histogram_central = self.histogram_central.astype(float)
452 histogram_subhalo = self.histogram_subhalo.astype(float).sum(axis=2)
454 if galaxy_type == "blue":
455 # Blue galaxies = (1 - quenched fraction) * histogram
456 histogram += histogram_central * (1 - quenched_frac_central)[None, :]
457 histogram += histogram_subhalo * (1 - quenched_frac_subhalo)[None, :]
458 elif galaxy_type == "red":
459 # Red galaxies = quenched fraction * histogram
460 histogram += histogram_central * quenched_frac_central[None, :]
461 histogram += histogram_subhalo * quenched_frac_subhalo[None, :]
462 else:
463 raise ValueError(
464 f"Galaxy type '{galaxy_type}' not supported. Use 'blue' or 'red'."
465 )
467 return histogram
469 def save(self, filename):
470 """
471 Saves the class instance to a file using pickle.
473 :param filename: Name of the file to save the instance to.
474 """
476 with open(filename, "wb") as f:
477 pickle.dump(self, f)
478 LOGGER.info(f"Saved HaloHistogram to {filename}")
480 def get_total_counts(self):
481 """
482 Returns the total counts of all histograms.
484 :return: Tuple of total counts for all, central, and subhalo histograms.
485 """
486 total_all = np.sum(self.histogram_all)
487 total_central = np.sum(self.histogram_central)
488 total_subhalo = np.sum(self.histogram_subhalo)
490 return total_all, total_central, total_subhalo
492 def get_z_bins(self):
493 """
494 Returns the redshift bin edges and centers.
495 :return: Array of redshift bin edges and centers.
496 """
497 return self.z_bins, (self.z_bins[:-1] + self.z_bins[1:]) / 2
499 def get_mass_bins(self):
500 """
501 Returns the mass bin edges and log centers.
502 :return: Array of mass bin edges and log centers.
503 """
504 centers_log = (self.log_mass_bins[:-1] + self.log_mass_bins[1:]) / 2
505 return self.mass_bins, np.power(10, centers_log)
507 def get_time_bins(self):
508 """
509 Returns the time since last merger bin edges and centers.
510 :return: Array of time since last merger bin edges and centers.
511 """
512 return self.time_bins, (self.time_bins[:-1] + self.time_bins[1:]) / 2
515def load_halo_histogram(filename):
516 """
517 Loads a HaloHistogram instance from a file.
519 :param filename: Name of the file to load the instance from.
520 :return: Loaded HaloHistogram instance.
521 """
522 with open(filename, "rb") as f:
523 halo_hist = pickle.load(f)
524 LOGGER.info(f"Loaded HaloHistogram from {filename}")
525 return halo_hist
528def read_pinocchio_config_details(pinocchio_output_filename):
529 """
530 Reads the Pinocchio output file to extract details about the simulation,
531 such as particle mass, minimum number of particles, and redshift range.
532 :param pinocchio_output_filename: Path to the Pinocchio output file.
533 :return: Tuple containing particle mass, minimum number of particles,
534 minimum redshift, and maximum redshift.
535 """
536 with open(pinocchio_output_filename) as f:
537 text = f.read()
539 m_part = float(
540 re.search("Particle Mass \\(Msun/h\\)\\s+(.+)\n", text).groups()[0],
541 )
542 min_num_part = int(
543 re.search("MinHaloMass \\(particles\\)\\s+(\\d+)\n", text).groups()[0],
544 )
545 z_max, z_min = (
546 float(i)
547 for i in re.search(
548 "The Past Light Cone will be reconstruct from z=(.+) to z=(.+)\n", text
549 ).groups()
550 )
552 return (
553 m_part,
554 min_num_part,
555 z_min,
556 z_max,
557 )
560def rebin_halo_histogram(
561 halo_hist, new_n_bins_mass=None, new_n_bins_z=None, new_n_bins_time=None
562):
563 """
564 Rebins a HaloHistogram to a new number of bins by combining neighboring bins.
565 The new number of bins must be a divisor of the current number of bins.
567 :param halo_hist: HaloHistogram instance to rebin
568 :param new_n_bins_mass: New number of mass bins (must divide current n_bins_mass)
569 :param new_n_bins_z: New number of redshift bins (must divide current n_bins_z)
570 :param new_n_bins_time: New number of time bins (must divide current n_bins_time)
571 :return: New HaloHistogram instance with rebinned histograms
572 """
573 if halo_hist.no_catalogs_added:
574 raise ValueError("Cannot rebin histogram before any catalogs have been added")
576 # Use current values if new values are not specified
577 if new_n_bins_mass is None:
578 new_n_bins_mass = halo_hist.n_bins_mass
579 if new_n_bins_z is None:
580 new_n_bins_z = halo_hist.n_bins_z
581 if new_n_bins_time is None:
582 new_n_bins_time = halo_hist.n_bins_time
584 # Check if new bin numbers are valid divisors
585 if halo_hist.n_bins_mass % new_n_bins_mass != 0:
586 raise ValueError(
587 f"new_n_bins_mass ({new_n_bins_mass}) must be a divisor of "
588 f"current n_bins_mass ({halo_hist.n_bins_mass})"
589 )
590 if halo_hist.n_bins_z % new_n_bins_z != 0:
591 raise ValueError(
592 f"new_n_bins_z ({new_n_bins_z}) must be a divisor of "
593 f"current n_bins_z ({halo_hist.n_bins_z})"
594 )
595 if halo_hist.n_bins_time % new_n_bins_time != 0:
596 raise ValueError(
597 f"new_n_bins_time ({new_n_bins_time}) must be a divisor of "
598 f"current n_bins_time ({halo_hist.n_bins_time})"
599 )
601 # Calculate rebin factors
602 mass_factor = halo_hist.n_bins_mass // new_n_bins_mass
603 z_factor = halo_hist.n_bins_z // new_n_bins_z
604 time_factor = halo_hist.n_bins_time // new_n_bins_time
606 # Create new histogram instance
607 new_hist = HaloHistogram(
608 z_min=halo_hist.z_min,
609 z_max=halo_hist.z_max,
610 min_mass=halo_hist.min_mass,
611 max_mass=halo_hist.max_mass,
612 min_time_since_merger=halo_hist.min_time_since_merger,
613 max_time_since_merger=halo_hist.max_time_since_merger,
614 n_bins_mass=new_n_bins_mass,
615 n_bins_z=new_n_bins_z,
616 n_bins_time=new_n_bins_time,
617 )
619 # Manually initialize the bins (same as what would happen in add_halo_catalog)
620 log_mass_bins = np.linspace(
621 new_hist.min_log_mass, new_hist.max_log_mass, new_n_bins_mass + 1
622 )
623 new_hist.mass_bins = np.power(10, log_mass_bins)
624 new_hist.z_bins = np.linspace(new_hist.z_min, new_hist.z_max, new_n_bins_z + 1)
625 new_hist.time_bins = np.linspace(
626 new_hist.min_time_since_merger,
627 new_hist.max_time_since_merger,
628 new_n_bins_time + 1,
629 )
630 new_hist.no_catalogs_added = False
632 # Rebin the histograms
633 new_hist.histogram_all = _rebin_2d(halo_hist.histogram_all, z_factor, mass_factor)
634 new_hist.histogram_central = _rebin_2d(
635 halo_hist.histogram_central, z_factor, mass_factor
636 )
637 new_hist.histogram_subhalo = _rebin_3d(
638 halo_hist.histogram_subhalo, z_factor, mass_factor, time_factor
639 )
641 return new_hist
644def _rebin_2d(histogram, factor_axis0, factor_axis1):
645 """
646 Helper function to rebin a 2D histogram by combining neighboring bins.
647 """
648 old_shape = histogram.shape
649 new_shape = (old_shape[0] // factor_axis0, old_shape[1] // factor_axis1)
651 # Reshape and sum over the dimensions to be combined
652 reshaped = histogram.reshape(new_shape[0], factor_axis0, new_shape[1], factor_axis1)
653 return reshaped.sum(axis=(1, 3))
656def _rebin_3d(histogram, factor_axis0, factor_axis1, factor_axis2):
657 """
658 Helper function to rebin a 3D histogram by combining neighboring bins.
659 """
660 old_shape = histogram.shape
661 new_shape = (
662 old_shape[0] // factor_axis0,
663 old_shape[1] // factor_axis1,
664 old_shape[2] // factor_axis2,
665 )
667 # Reshape and sum over the dimensions to be combined
668 reshaped = histogram.reshape(
669 new_shape[0],
670 factor_axis0,
671 new_shape[1],
672 factor_axis1,
673 new_shape[2],
674 factor_axis2,
675 )
676 return reshaped.sum(axis=(1, 3, 5))