Build and Inspect SFC Fields#
SFCProjection maps a particle catalogue onto scaling-function coefficients (SFCs). The resulting SFCField is the reusable field object consumed by PyHermes windows and correlation tasks.
This notebook starts from the same public Quijote halo catalogue and YAML as quick_start.ipynb. It then builds the complete example field set in one top-to-bottom run: the standard and higher-resolution number fields, a halo-mass field, redshift-space fields, and a matching random field. It also reconstructs the represented density field and compares it with direct halo binning. For catalogue formats and conversion, see particle_io.ipynb.
[1]:
from pathlib import Path
import os
import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable
from matplotlib.colors import LinearSegmentedColormap, Normalize, PowerNorm
import numpy as np
from scipy.ndimage import gaussian_filter
from pyhermes.base import SFCProjection
from pyhermes.io import read_particle_data
from pyhermes.param.parambase import read_param
from pyhermes.utils.redshift_space import hubble_at_redshift, redshift_space_positions
from pyhermes.utils.sampling import random_box_positions
cwd = Path.cwd().resolve()
if cwd.name == "notebooks" and cwd.parent.name == "examples":
examples_dir = cwd.parent
elif cwd.name == "examples":
examples_dir = cwd
elif (cwd / "examples").is_dir():
examples_dir = cwd / "examples"
else:
raise RuntimeError("Run this notebook from the project root, examples/, or examples/notebooks/.")
os.chdir(examples_dir)
Path("./output").mkdir(exist_ok=True)
Path("./figs").mkdir(exist_ok=True)
print(f"Working directory: {Path.cwd()}")
Working directory: /Users/xutengpeng/xutp/PycharmProjects/PyHermes/examples
1. Project the configuration-defined catalogue#
The YAML is the source of truth for the catalogue URL and cache, box, multiresolution basis, normalization, companion particle data, and output paths. The equivalent command-line run is:
python scripts/run_sfc_projection.py configs/param_sfc_projection.yaml
[2]:
projection_params = read_param("./configs/param_sfc_projection.yaml")
projection_config = projection_params["SFCProjection"]
projection_task = SFCProjection(projection_params)
field_j8 = projection_task.run(overwrite=True)
11:37:13 - INFO - pyhermes.param.parambase:ParamBase - Reading configure file: './configs/param_sfc_projection.yaml'
11:37:13 - INFO - pyhermes.param.parambase:ParamBase - Input parameter file format is YAML
11:37:13 - INFO - pyhermes.param.parambase:ParamBase - Set default parameters of module <base> ...
11:37:13 - INFO - pyhermes.param.parambase:ParamBase - Default 'SFCProjection.fin.path' from 'empty' to 'https://pyhermes.astroslacker.com/downloads/group_tab_004.0'
11:37:13 - INFO - pyhermes.param.parambase:ParamBase - Default 'SFCProjection.fin.format' from 'bin' to 'fof'
11:37:13 - INFO - pyhermes.param.parambase:ParamBase - Default 'SFCProjection.fin.download.cache_path' from 'empty' to './data/quijote_halos/8000/groups_004/group_tab_004.0'
11:37:13 - INFO - pyhermes.param.parambase:ParamBase - Default 'SFCProjection.fin.download.sha256' from 'empty' to '4a1c6ca4f6747a70e9e552685226ecf5d678c6c97551e2caa7cc3883502eac85'
11:37:13 - INFO - pyhermes.param.parambase:ParamBase - Default 'SFCProjection.threads' from '1' to '2'
11:37:13 - INFO - pyhermes.param.parambase:ParamBase - Default 'SFCProjection.save_particle_data' from 'False' to 'True'
11:37:13 - INFO - pyhermes.param.parambase:ParamBase - Default 'SFCProjection.particle_data_path' from 'empty' to './output/quijote8000_snap004_particles.npz'
11:37:13 - INFO - pyhermes.param.parambase:ParamBase - Default 'SFCProjection.fout_path' from 'empty' to './output/quijote8000_snap004_sfc.pkl'
11:37:13 - INFO - pyhermes.pipeline.pipeline:SFCProjection - SFCProjection runtime configuration: running on 1 MPI ranks with 2 threads per rank
11:37:13 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Preparing SFCProjection input fields ...
11:37:13 - INFO - pyhermes.pipeline.pipeline:SFCProjection - J=8, L=256, box_size=1000.0, phi_resolution=1024, wavelet_mode=db2, wavelet_level=10
11:37:13 - INFO - pyhermes.io.resources:resolve_data_path - Using cached particle data: data/quijote_halos/8000/groups_004/group_tab_004.0
11:37:13 - INFO - pyhermes.io.readers:read_particle_data - Selected input particle format: fof
11:37:13 - INFO - pyhermes.io.readers:read_fof - Reading Quijote FoF halo data from ---> data/quijote_halos/8000/groups_004/group_tab_004.0 <---
11:37:13 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Input particles ready | source=file=path=https://pyhermes.astroslacker.com/downloads/group_tab_004.0 format=fof | particle_count=406728 | catalog_weight_key=None | field_value_key=None | weight_normalization=catalog | catalog_weight_sum=4.067280e+05 | field_integral=1.000000e+00
11:37:13 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Saved particle positions, catalogue weights, and field values to ./output/quijote8000_snap004_particles.npz
11:37:13 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Single process mode
11:37:15 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for scaling function: 1.7707 sec
11:37:15 - WARNING - pyhermes.io.funcs:check_fout - Output file './output/quijote8000_snap004_sfc.pkl' already exists and will be overwritten (overwrite=True).
11:37:15 - INFO - pyhermes.io.base:SFCField - Writing SFCField data to ---> ./output/quijote8000_snap004_sfc.pkl <---
11:37:15 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for task: 2.1628 sec
[3]:
print(f"epsilon shape: {field_j8.epsilon.shape}")
print(f"J={field_j8.J}, cells per axis={field_j8.L}, cell spacing={field_j8.box_size / field_j8.L:.3f} Mpc/h")
print(f"field integral: {field_j8.field_integral:.6g}")
print(f"normalization: {field_j8.weight_normalization}")
epsilon shape: (256, 256, 256)
J=8, cells per axis=256, cell spacing=3.906 Mpc/h
field integral: 1
normalization: catalog
3. Reconstruct the projected density field#
The stored coefficients are expansion coefficients, not density samples. At multiresolution level \(j\), PyHermes reconstructs the represented field as
\[F_j(\mathbf{x}) = \sum_{\boldsymbol{\ell}} \epsilon^j_{\boldsymbol{\ell}}\,\phi^j_{\boldsymbol{\ell}}(\mathbf{x}),
\qquad
\phi^j_{\boldsymbol{\ell}}(\mathbf{x}) = \prod_{d=1}^{3}\phi^j_{\ell_d}(x_d).\]
For a slab centred on \(z_0\) with thickness \(\Delta z\), the plotted projected field is
\[\Sigma_j(x,y) = \int_{z_0-\Delta z/2}^{z_0+\Delta z/2} F_j(x,y,z)\,\mathrm{d}z.\]
Panels A-C evaluate this expression at \(J=7,8,9\). Panel D directly bins the same halo subset in the same slab. All panels use the same catalogue normalization, display smoothing, colormap, and color normalization, so their structures can be compared directly.
[7]:
field_j7 = project_arrays(pos, J=7)
box_size = float(projection_config["box_size"])
x_range = (300.0, 700.0)
y_range = (300.0, 700.0)
slice_z = 500.0
slice_thickness = 50.0
def periodic_slab_mask(values, center, thickness, period):
periodic_delta = (values - center + 0.5 * period) % period - 0.5 * period
return np.abs(periodic_delta) <= 0.5 * thickness
halo_mask = (
periodic_slab_mask(pos[:, 2], slice_z, slice_thickness, box_size)
& (pos[:, 0] >= x_range[0])
& (pos[:, 0] <= x_range[1])
& (pos[:, 1] >= y_range[0])
& (pos[:, 1] <= y_range[1])
)
halo_cut = pos[halo_mask]
print(f"Haloes in the displayed slab: {halo_cut.shape[0]:,}")
11:37:19 - INFO - pyhermes.param.parambase:ParamBase - Set default parameters of module <base> ...
11:37:19 - INFO - pyhermes.pipeline.pipeline:SFCProjection - SFCProjection runtime configuration: running on 1 MPI ranks with 2 threads per rank
11:37:19 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Preparing SFCProjection input fields ...
11:37:19 - INFO - pyhermes.pipeline.pipeline:SFCProjection - J=7, L=128, box_size=1000.0, phi_resolution=1024, wavelet_mode=db2, wavelet_level=10
11:37:19 - INFO - pyhermes.pipeline.pipeline:SFCProjection - No catalog_weight provided; using unit catalogue weights for 406728 particles.
11:37:19 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Input particles ready | source=custom particle_pos array | particle_count=406728 | catalog_weight_key=None | field_value_key=None | weight_normalization=catalog | catalog_weight_sum=4.067280e+05 | field_integral=1.000000e+00
11:37:19 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Single process mode
11:37:19 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for scaling function: 0.0796 sec
11:37:19 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for task: 0.0895 sec
Haloes in the displayed slab: 3,351
[8]:
def reconstructed_projected_density(field, nx=256, nz=16, chunk_size=262_144):
x_edges = np.linspace(*x_range, nx + 1)
y_edges = np.linspace(*y_range, nx + 1)
x = 0.5 * (x_edges[:-1] + x_edges[1:])
y = 0.5 * (y_edges[:-1] + y_edges[1:])
dz = slice_thickness / nz
z = (slice_z - 0.5 * slice_thickness + (np.arange(nz) + 0.5) * dz) % box_size
xx, yy = np.meshgrid(x, y, indexing="ij")
xy = np.column_stack([xx.ravel(), yy.ravel()])
projected = np.zeros(xx.size, dtype=np.float64)
for z_value in z:
sample_pos = np.column_stack([xy, np.full(xy.shape[0], z_value)])
for start in range(0, sample_pos.shape[0], chunk_size):
stop = min(start + chunk_size, sample_pos.shape[0])
projected[start:stop] += (
field.field_density_at_pos(sample_pos[start:stop], value_unit="physical") * dz
)
return projected.reshape(nx, nx)
def halo_projected_density_reference(halo_positions, particle_count, nx=256):
pixel_area = ((x_range[1] - x_range[0]) / nx) * ((y_range[1] - y_range[0]) / nx)
weights = np.full(halo_positions.shape[0], 1.0 / particle_count)
histogram, _, _ = np.histogram2d(
halo_positions[:, 0],
halo_positions[:, 1],
bins=nx,
range=[x_range, y_range],
weights=weights,
)
return histogram / pixel_area
def smooth_for_display(image, nx, smooth_mpc):
if smooth_mpc <= 0.0:
return image
sigma_pixels = smooth_mpc * nx / (x_range[1] - x_range[0])
return gaussian_filter(image, sigma=sigma_pixels, mode="nearest")
[9]:
reconstruct_nx = 256
reconstruct_nz = 16
display_smooth_mpc = 1.0
use_power_norm = True
gamma = 0.5
reconstructed_maps = [
reconstructed_projected_density(field, nx=reconstruct_nx, nz=reconstruct_nz)
for field in (field_j7, field_j8, field_j9)
]
halo_reference = halo_projected_density_reference(
halo_cut,
field_j8.particle_count,
nx=reconstruct_nx,
)
display_maps = [
smooth_for_display(np.clip(image, 0.0, None), reconstruct_nx, display_smooth_mpc)
for image in (*reconstructed_maps, halo_reference)
]
positive_values = np.concatenate([image[image > 0.0] for image in display_maps])
vmax = np.percentile(positive_values, 99.7)
density_norm = (
PowerNorm(gamma=gamma, vmin=0.0, vmax=vmax)
if use_power_norm
else Normalize(vmin=0.0, vmax=vmax)
)
density_cmap = LinearSegmentedColormap.from_list(
"reconstructed_density",
["#ffffff", "#d7f0f0", "#7cc4c4", "#2f6f9f", "#f2b35e", "#b8322a"],
)
fig, axes = plt.subplots(2, 2, figsize=(8.8, 8.0), sharex=True, sharey=True)
extent = (*x_range, *y_range)
panels = [
(axes[0, 0], display_maps[0], r"A: PyHermes $J=7$", "green"),
(axes[0, 1], display_maps[1], r"B: PyHermes $J=8$", "green"),
(axes[1, 0], display_maps[2], r"C: PyHermes $J=9$", "green"),
(axes[1, 1], display_maps[3], "D: direct halo binning", "red"),
]
for ax, image, label, edgecolor in panels:
ax.imshow(
image.T,
origin="lower",
extent=extent,
cmap=density_cmap,
norm=density_norm,
interpolation="bilinear",
)
ax.text(
0.95,
0.95,
label,
transform=ax.transAxes,
ha="right",
va="top",
fontsize=11,
bbox=dict(boxstyle="round,pad=0.25", facecolor="white", edgecolor=edgecolor, linewidth=1.2),
)
for ax in axes.ravel():
ax.set_xlim(*x_range)
ax.set_ylim(*y_range)
ax.set_aspect("equal")
ax.tick_params(labelsize=10)
ax.label_outer()
fig.supxlabel(r"$x\ [h^{-1}\mathrm{Mpc}]$", fontsize=14)
fig.supylabel(r"$y\ [h^{-1}\mathrm{Mpc}]$", fontsize=14)
fig.subplots_adjust(left=0.12, right=0.84, bottom=0.10, top=0.96, wspace=0.12, hspace=0.12)
mappable = ScalarMappable(norm=density_norm, cmap=density_cmap)
mappable.set_array([])
colorbar_ax = fig.add_axes([0.91, 0.20, 0.025, 0.60])
colorbar = fig.colorbar(mappable, cax=colorbar_ax)
colorbar.set_label(r"projected normalized density [$(h^{-1}\mathrm{Mpc})^{-2}$]", fontsize=11)
colorbar.ax.tick_params(labelsize=9)
fig.savefig(
"./figs/sfc_field_reconstructed_density_j7_j8_j9_halo_binned.png",
dpi=200,
bbox_inches="tight",
)
plt.show()
4. Transform coordinates, then project once#
Redshift-space distortions enter through particle coordinates. The next cell constructs both the common distant-observer line of sight along z and a diagonal line of sight, then projects number and mass fields with the same SFC basis. The shifted particle positions are saved as companion NPZ files for later particle-centred statistics.
[10]:
redshift = 0.0
hubble = hubble_at_redshift(redshift)
rsd_pos_z = redshift_space_positions(pos, velocity, box_size, hubble, redshift, los="z")
rsd_pos_diag = redshift_space_positions(pos, velocity, box_size, hubble, redshift, los=[1, 1, 1])
rsd_field = project_arrays(
rsd_pos_z,
fout_path="./output/quijote8000_snap004_rsd_sfc.pkl",
particle_data_path="./data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd.npz",
)
rsd_mass_field = project_arrays(
rsd_pos_z,
field_value=mass,
fout_path="./output/quijote8000_snap004_rsd_sfc_massweight.pkl",
particle_data_path="./data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd_massweight.npz",
)
rsd_diag_field = project_arrays(
rsd_pos_diag,
fout_path="./output/quijote8000_snap004_rsd_diag_sfc.pkl",
particle_data_path="./data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd_diag.npz",
)
11:37:22 - INFO - pyhermes.param.parambase:ParamBase - Set default parameters of module <base> ...
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - SFCProjection runtime configuration: running on 1 MPI ranks with 2 threads per rank
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Preparing SFCProjection input fields ...
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - J=8, L=256, box_size=1000.0, phi_resolution=1024, wavelet_mode=db2, wavelet_level=10
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - No catalog_weight provided; using unit catalogue weights for 406728 particles.
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Input particles ready | source=custom particle_pos array | particle_count=406728 | catalog_weight_key=None | field_value_key=None | weight_normalization=catalog | catalog_weight_sum=4.067280e+05 | field_integral=1.000000e+00
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Saved particle positions, catalogue weights, and field values to ./data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd.npz
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Single process mode
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for scaling function: 0.1199 sec
11:37:22 - WARNING - pyhermes.io.funcs:check_fout - Output file './output/quijote8000_snap004_rsd_sfc.pkl' already exists and will be overwritten (overwrite=True).
11:37:22 - INFO - pyhermes.io.base:SFCField - Writing SFCField data to ---> ./output/quijote8000_snap004_rsd_sfc.pkl <---
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for task: 0.2533 sec
11:37:22 - INFO - pyhermes.param.parambase:ParamBase - Set default parameters of module <base> ...
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - SFCProjection runtime configuration: running on 1 MPI ranks with 2 threads per rank
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Preparing SFCProjection input fields ...
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - J=8, L=256, box_size=1000.0, phi_resolution=1024, wavelet_mode=db2, wavelet_level=10
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - No catalog_weight provided; using unit catalogue weights for 406728 particles.
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Input particles ready | source=custom particle_pos array | particle_count=406728 | catalog_weight_key=None | field_value_key=custom | weight_normalization=catalog | catalog_weight_sum=4.067280e+05 | field_integral=4.830485e+13
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Saved particle positions, catalogue weights, and field values to ./data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd_massweight.npz
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Single process mode
11:37:22 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for scaling function: 0.1259 sec
11:37:22 - WARNING - pyhermes.io.funcs:check_fout - Output file './output/quijote8000_snap004_rsd_sfc_massweight.pkl' already exists and will be overwritten (overwrite=True).
11:37:22 - INFO - pyhermes.io.base:SFCField - Writing SFCField data to ---> ./output/quijote8000_snap004_rsd_sfc_massweight.pkl <---
11:37:23 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for task: 0.2346 sec
11:37:23 - INFO - pyhermes.param.parambase:ParamBase - Set default parameters of module <base> ...
11:37:23 - INFO - pyhermes.pipeline.pipeline:SFCProjection - SFCProjection runtime configuration: running on 1 MPI ranks with 2 threads per rank
11:37:23 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Preparing SFCProjection input fields ...
11:37:23 - INFO - pyhermes.pipeline.pipeline:SFCProjection - J=8, L=256, box_size=1000.0, phi_resolution=1024, wavelet_mode=db2, wavelet_level=10
11:37:23 - INFO - pyhermes.pipeline.pipeline:SFCProjection - No catalog_weight provided; using unit catalogue weights for 406728 particles.
11:37:23 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Input particles ready | source=custom particle_pos array | particle_count=406728 | catalog_weight_key=None | field_value_key=None | weight_normalization=catalog | catalog_weight_sum=4.067280e+05 | field_integral=1.000000e+00
11:37:23 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Saved particle positions, catalogue weights, and field values to ./data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd_diag.npz
11:37:23 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Single process mode
11:37:23 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for scaling function: 0.1131 sec
11:37:23 - WARNING - pyhermes.io.funcs:check_fout - Output file './output/quijote8000_snap004_rsd_diag_sfc.pkl' already exists and will be overwritten (overwrite=True).
11:37:23 - INFO - pyhermes.io.base:SFCField - Writing SFCField data to ---> ./output/quijote8000_snap004_rsd_diag_sfc.pkl <---
11:37:23 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for task: 0.2646 sec
[11]:
slice_resolution = 160
z_slice = 0.5 * box_size
slice_axis = (np.arange(slice_resolution) + 0.5) * box_size / slice_resolution
slice_x, slice_y = np.meshgrid(slice_axis, slice_axis, indexing="xy")
slice_pos = np.column_stack([
slice_x.ravel(),
slice_y.ravel(),
np.full(slice_x.size, z_slice),
])
def plane_density_contrast(field):
density = field.field_density_at_pos(slice_pos, value_unit="physical")
mean_density = field.field_mean_density(value_unit="physical")
return (density / mean_density - 1.0).reshape(slice_x.shape)
delta_real = plane_density_contrast(field_j8)
delta_rsd = plane_density_contrast(rsd_field)
[12]:
vmax = np.percentile(np.abs(np.concatenate([delta_real.ravel(), delta_rsd.ravel()])), 99.0)
fig, axes = plt.subplots(1, 2, figsize=(11.2, 4.8), sharex=True, sharey=True, constrained_layout=True)
for ax, delta, title in zip(axes, [delta_real, delta_rsd], ["Real space", "Redshift space"]):
image = ax.imshow(
delta,
origin="lower",
extent=(0.0, box_size, 0.0, box_size),
cmap="coolwarm",
vmin=-vmax,
vmax=vmax,
interpolation="nearest",
)
ax.set_title(title)
ax.set_xlabel(r"$x$ [Mpc/$h$]")
axes[0].set_ylabel(r"$y$ [Mpc/$h$]")
fig.colorbar(image, ax=axes, label=r"$\delta_h$", shrink=0.88)
fig.savefig("./figs/sfc_projection_real_rsd_slice.png", dpi=180, bbox_inches="tight")
plt.show()
5. Build the matching random field#
Correlation estimators need a random catalogue with the same selection function as the data. For this periodic-box example, the selection is uniform. The fixed seed makes the 10-million-point random field exactly reproducible, and its companion NPZ stores the generated positions.
[13]:
random_count = 10_000_000
random_pos = random_box_positions(count=random_count, box_size=box_size, seed=42).astype(
np.float32,
copy=False,
)
random_field = project_arrays(
random_pos,
fout_path="./output/random_sfc.pkl",
particle_data_path="./data/random_1e7.npz",
)
del random_pos
11:37:24 - INFO - pyhermes.param.parambase:ParamBase - Set default parameters of module <base> ...
11:37:24 - INFO - pyhermes.pipeline.pipeline:SFCProjection - SFCProjection runtime configuration: running on 1 MPI ranks with 2 threads per rank
11:37:24 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Preparing SFCProjection input fields ...
11:37:24 - INFO - pyhermes.pipeline.pipeline:SFCProjection - J=8, L=256, box_size=1000.0, phi_resolution=1024, wavelet_mode=db2, wavelet_level=10
11:37:24 - INFO - pyhermes.pipeline.pipeline:SFCProjection - No catalog_weight provided; using unit catalogue weights for 10000000 particles.
11:37:24 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Input particles ready | source=custom particle_pos array | particle_count=10000000 | catalog_weight_key=None | field_value_key=None | weight_normalization=catalog | catalog_weight_sum=1.000000e+07 | field_integral=1.000000e+00
11:37:24 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Saved particle positions, catalogue weights, and field values to ./data/random_1e7.npz
11:37:24 - INFO - pyhermes.pipeline.pipeline:SFCProjection - Single process mode
11:37:27 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for scaling function: 2.7688 sec
11:37:27 - WARNING - pyhermes.io.funcs:check_fout - Output file './output/random_sfc.pkl' already exists and will be overwritten (overwrite=True).
11:37:27 - INFO - pyhermes.io.base:SFCField - Writing SFCField data to ---> ./output/random_sfc.pkl <---
11:37:27 - INFO - pyhermes.pipeline.pipeline:SFCProjection - The time for task: 3.2446 sec
6. Inspect the companion data and generated products#
When save_particle_data: true, SFCProjection writes the exact positions and scalar factors used by that projection. This compact companion lets particle-centred statistics recover their catalogue consistently. It is a task product derived from the public FoF source, not another distributed source catalogue.
[14]:
particle_data = field_j8.get_particle_data()
print(f"companion arrays: {sorted(particle_data)}")
print(f"positions: {particle_data['pos'].shape}")
print(f"projection-weight sum: {particle_data['projection_weight'].sum():.6g}")
companion arrays: ['catalog_weight', 'field_value', 'normalized_catalog_weight', 'pos', 'projection_weight']
positions: (406728, 3)
projection-weight sum: 1
[15]:
generated_products = [
"./output/quijote8000_snap004_sfc.pkl",
"./output/quijote8000_snap004_particles.npz",
"./output/quijote8000_snap004_sfc_J9.pkl",
"./output/quijote8000_snap004_sfc_massweight.pkl",
"./output/quijote8000_snap004_rsd_sfc.pkl",
"./output/quijote8000_snap004_rsd_sfc_massweight.pkl",
"./output/quijote8000_snap004_rsd_diag_sfc.pkl",
"./output/random_sfc.pkl",
"./data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd.npz",
"./data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd_massweight.npz",
"./data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd_diag.npz",
"./data/random_1e7.npz",
]
for product in generated_products:
path = Path(product)
print(f"{path}: {path.stat().st_size / 1024**2:.1f} MiB")
output/quijote8000_snap004_sfc.pkl: 128.0 MiB
output/quijote8000_snap004_particles.npz: 10.9 MiB
output/quijote8000_snap004_sfc_J9.pkl: 1024.0 MiB
output/quijote8000_snap004_sfc_massweight.pkl: 128.0 MiB
output/quijote8000_snap004_rsd_sfc.pkl: 128.0 MiB
output/quijote8000_snap004_rsd_sfc_massweight.pkl: 128.0 MiB
output/quijote8000_snap004_rsd_diag_sfc.pkl: 128.0 MiB
output/random_sfc.pkl: 128.0 MiB
data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd.npz: 10.9 MiB
data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd_massweight.npz: 10.9 MiB
data/quijote_halos/8000/groups_004/group_tab_004.pos.rsd_diag.npz: 10.9 MiB
data/random_1e7.npz: 267.0 MiB
Next steps#
This notebook has generated the same reusable example fields as scripts/prepare_sfc_fields.py; running that script again is unnecessary. Continue with window.ipynb for the common field-operator language. From there, choose physical_fields.ipynb for differential and inverse-Laplacian operators, counting.ipynb for one-point distributions, or the correlation notebooks for two- and three-point statistics. The script remains useful as a non-interactive batch entry point.