From a03f9a288eec8f97a84a71742514c1ae524d5ec1 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Wed, 4 Jun 2025 13:44:21 +0200 Subject: [PATCH 1/6] monofonic support for A_s or sigma8 --- sbmy_control/parameters_monofonic.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/sbmy_control/parameters_monofonic.py b/sbmy_control/parameters_monofonic.py index e4a2169..d9cf896 100644 --- a/sbmy_control/parameters_monofonic.py +++ b/sbmy_control/parameters_monofonic.py @@ -95,8 +95,6 @@ def get_config_from_dict(monofonic_dict): "Omega_L": monofonic_dict["Omega_q"], "H0": monofonic_dict["h"]*100., "n_s": monofonic_dict["n_s"], - "sigma_8": monofonic_dict["sigma8"], - "A_s": monofonic_dict["A_s"], "Tcmb": monofonic_dict["Tcmb"], "k_p": monofonic_dict["k_p"], "N_ur": monofonic_dict["N_ur"], @@ -110,7 +108,14 @@ def get_config_from_dict(monofonic_dict): "transfer": "CLASS", "ztarget": 2.5, "ZeroRadiation": False - } + } + + if "A_s" in monofonic_dict.keys() and monofonic_dict["A_s"] is not None: + config["cosmology"]["A_s"] = monofonic_dict["A_s"] + elif "sigma8" in monofonic_dict.keys() and monofonic_dict["sigma8"] is not None: + config["cosmology"]["sigma8"] = monofonic_dict["sigma8"] + else: + raise KeyError("Either A_s or sigma8 must be provided in the cosmology parameters.") config["random"] = { "generator": "NGENIC", From 64e90c84a2f9213ee2f7b8d52bf718569da3c33f Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Tue, 24 Jun 2025 11:42:54 +0200 Subject: [PATCH 2/6] convert snapshots --- .../scripts/convert_snapshot_to_density.py | 73 +++++++++++++++++-- 1 file changed, 66 insertions(+), 7 deletions(-) diff --git a/sbmy_control/scripts/convert_snapshot_to_density.py b/sbmy_control/scripts/convert_snapshot_to_density.py index dadefed..869ee90 100644 --- a/sbmy_control/scripts/convert_snapshot_to_density.py +++ b/sbmy_control/scripts/convert_snapshot_to_density.py @@ -34,6 +34,55 @@ def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, F.write(output_path) print("Density written to", output_path) print("Done.") + + +def convert_snapshots_to_density(snapshot_base_path, output_path, N=None, corner=(0.0, 0.0, 0.0)): + """ + Convert multiple snapshots to density fields. + + Parameters + ---------- + snapshot_base_path : str + Base path for the snapshot files. + output_path : str + Path to the output density file. + N : int + Size of the density field grid (N x N x N). + corner : tuple of float + Corner of the box (x, y, z). + """ + # Get all snapshot files with path "snapshot_base_path*" + import glob + snapshot_files = glob.glob(snapshot_base_path + "*") + if not snapshot_files: + raise FileNotFoundError(f"No snapshot files found at {snapshot_base_path}") + + A = None + F = None + + for snapshot_path in snapshot_files: + # Read the snapshot + print(f"Reading snapshot {snapshot_path}...") + snap = read_snapshot(snapshot_path) + + if N is None: + N = snap.Np0 + + # Calculate density + print("Calculating density...") + F = get_density_pm_snapshot(snap, N, N, N, corner[0], corner[1], corner[2]) + + if A is None: + A = F.data + else: + A += F.data + + F.data = A # Combine densities from all snapshots + # Write density to file + print("Writing density...") + F.write(output_path) + print("Density written to", output_path) + print("Done.") def console_main(): @@ -70,13 +119,23 @@ def console_main(): args = parser.parse_args() - convert_snapshot_to_density( - snapshot_path=args.snapshot, - output_path=args.output, - N=args.N, - corner=args.corner, - ) + if args.snapshot.endswith("h5") or args.snapshot.endswith("hdf5") or args.snapshot.endswith("gadget3"): + # If the snapshot is a single file, convert it to density + convert_snapshot_to_density( + snapshot_path=args.snapshot, + output_path=args.output, + N=args.N, + corner=args.corner, + ) + else: + # If the snapshot is a base path, convert all snapshots to density + convert_snapshots_to_density( + snapshot_base_path=args.snapshot, + output_path=args.output, + N=args.N, + corner=args.corner, + ) if __name__ == "__main__": - console_main() \ No newline at end of file + console_main() From 1db239707dce63fdab8ad622c9fb01d63be1c100 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 30 Jun 2025 14:27:20 +0200 Subject: [PATCH 3/6] tpoy --- sbmy_control/simbelmyne.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbmy_control/simbelmyne.py b/sbmy_control/simbelmyne.py index a106ae9..ac2a6ba 100644 --- a/sbmy_control/simbelmyne.py +++ b/sbmy_control/simbelmyne.py @@ -16,7 +16,7 @@ def main_simbelmyne(parsed_args): if isfile(log_file): # Remove the preexisting log file to allow for the progress_bar to be run normally from os import remove - oremove(log_file) + remove(log_file) command_args = ["simbelmyne", paramfile, log_file] From 6b22752ca0d22a809c8cde006cff3a3c926268b9 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 30 Jun 2025 14:27:46 +0200 Subject: [PATCH 4/6] monofonic slurm script with sim name --- sbmy_control/monofonic.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sbmy_control/monofonic.py b/sbmy_control/monofonic.py index 134b55b..861f14b 100644 --- a/sbmy_control/monofonic.py +++ b/sbmy_control/monofonic.py @@ -32,7 +32,7 @@ def main_monofonic(parsed_args): print_message("Running monofonic in slurm mode.", 1, "monofonic", verbose=parsed_args.verbose) slurm_dict=parse_arguments_slurm(parsed_args) main_dict=parse_arguments_main(parsed_args) - slurm_script = slurm_dict["scripts"]+"monofonic.sh" + slurm_script = slurm_dict["scripts"]+f"monofonic_{main_dict['simname']}.sh" if not isfile(slurm_script) or parsed_args.force: print_message(f"SLURM script {slurm_script} does not exist (or forced). Creating it.", 2, "monofonic", verbose=parsed_args.verbose) From 18caa0d30aa854ce811ec811280e8d50cdf27e7f Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Tue, 1 Jul 2025 15:15:35 +0200 Subject: [PATCH 5/6] better k_modes for power spectra --- sbmy_control/analysis/power_spectrum.py | 49 +++++++++++++++++-------- 1 file changed, 33 insertions(+), 16 deletions(-) diff --git a/sbmy_control/analysis/power_spectrum.py b/sbmy_control/analysis/power_spectrum.py index 7a92370..d97e24f 100644 --- a/sbmy_control/analysis/power_spectrum.py +++ b/sbmy_control/analysis/power_spectrum.py @@ -1,7 +1,7 @@ import numpy as np import os -kmin = 1e-1 +nfirst_kmodes = 10 kmax = 2e0 Nk = 50 AliasingCorr=False @@ -24,12 +24,17 @@ def crop_field(field, Ncrop): field.L2 = field.N2*d2 -def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None): +def get_power_spectrum(field, nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, G=None): from pysbmy.power import PowerSpectrum from pysbmy.fft import FourierGrid from pysbmy.correlations import get_autocorrelation if G is None: + default_k_modes = PowerSpectrum(field.L0,field.L1,field.L2,int(field.N0),int(field.N1),int(field.N2),).FourierGrid.k_modes + if kmax is None: + kmax = default_k_modes[-1] + if kmax > default_k_modes[-1]: + raise ValueError(f"kmax ({kmax}) is larger than the maximum k mode available in the field ({default_k_modes[-1]}).") G = FourierGrid( field.L0, field.L1, @@ -37,8 +42,8 @@ def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None): int(field.N0), int(field.N1), int(field.N2), - k_modes=np.concat([PowerSpectrum(field.L0,field.L1,field.L2,int(field.N0),int(field.N1),int(field.N2),).FourierGrid.k_modes[:10],np.logspace( - np.log10(kmin), + k_modes=np.concat([default_k_modes[:nfirst_kmodes],np.logspace( + np.log10(default_k_modes[nfirst_kmodes]), np.log10(kmax), Nk, )]), @@ -51,13 +56,18 @@ def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None): return G, k, Pk -def get_cross_correlations(field_A, field_B, kmin=kmin, kmax=kmax, Nk=Nk, G=None): +def get_cross_correlations(field_A, field_B, nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, G=None): from pysbmy.power import PowerSpectrum from pysbmy.fft import FourierGrid from pysbmy.correlations import get_crosscorrelation if G is None: + default_k_modes = PowerSpectrum(field_A.L0,field_A.L1,field_A.L2,int(field_A.N0),int(field_A.N1),int(field_A.N2),).FourierGrid.k_modes + if kmax is None: + kmax = default_k_modes[-1] + if kmax > default_k_modes[-1]: + raise ValueError(f"kmax ({kmax}) is larger than the maximum k mode available in the field ({default_k_modes[-1]}).") G = FourierGrid( field_A.L0, field_A.L1, @@ -65,8 +75,8 @@ def get_cross_correlations(field_A, field_B, kmin=kmin, kmax=kmax, Nk=Nk, G=None int(field_A.N0), int(field_A.N1), int(field_A.N2), - k_modes=np.concat([PowerSpectrum(field_A.L0,field_A.L1,field_A.L2,int(field_A.N0),int(field_A.N1),int(field_A.N2),).FourierGrid.k_modes[:10],np.logspace( - np.log10(kmin), + k_modes=np.concat([default_k_modes[:nfirst_kmodes],np.logspace( + np.log10(default_k_modes[nfirst_kmodes]), np.log10(kmax), Nk, )]), @@ -104,7 +114,7 @@ def plot_power_spectra(filenames, yticks = np.linspace(0.9,1.1,11), bound1=0.01, bound2=0.02, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, figsize=(8,4), @@ -139,7 +149,7 @@ def plot_power_spectra(filenames, color=colors[i], linestyle=linestyles[i], marker=markers[i],), - power_args=dict(kmin=kmin, + power_args=dict(nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk), ) @@ -178,7 +188,7 @@ def plot_cross_correlations(filenames_A, yticks = np.linspace(0.99,1.001,12), bound1=0.001, bound2=0.002, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, figsize=(8,4), @@ -217,7 +227,7 @@ def plot_cross_correlations(filenames_A, color=colors[i], linestyle=linestyles[i], marker=markers[i],), - power_args=dict(kmin=kmin, + power_args=dict(nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk), ) @@ -283,8 +293,15 @@ def console_main(): parser.add_argument('-yrp', '--ylim_power', type=float, nargs=2, default=[0.9,1.1], help='Y-axis limits.') parser.add_argument('-yrc', '--ylim_corr', type=float, nargs=2, default=[0.99,1.001], help='Y-axis limits.') parser.add_argument('--crop', type=int, default=None, help='Remove the outter N pixels of the fields.') + parser.add_argument('--nfirst_kmodes', type=float, default=10, help='First k modes from all available k modes to be used in the power spectrum before the geomspace part.') + parser.add_argument('--kmax', type=float, default=None, help='Maximum k value for the power spectrum geomspace part. If None, it will be set to the maximum k mode available in the field.') + parser.add_argument('--Nk', type=int, default=50, help='Number of k values for the power spectrum geomspace part.') args = parser.parse_args() + + nfirst_kmodes = args.nfirst_kmodes + kmax = args.kmax + Nk = args.Nk if not args.power_spectrum and not args.cross_correlation: print('You must choose between power_spectrum and cross_correlation.') @@ -298,7 +315,7 @@ def console_main(): from pysbmy.field import read_field F_ref = read_field(args.directory+args.reference) crop_field(F_ref, args.crop) - G, _, Pk_ref = get_power_spectrum(F_ref, kmin=kmin, kmax=kmax, Nk=Nk) + G, _, Pk_ref = get_power_spectrum(F_ref, nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk) else: Pk_ref = None G = None @@ -325,7 +342,7 @@ def console_main(): yticks = yticks_power, bound1=0.01, bound2=0.02, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, ax=axes[0], @@ -344,7 +361,7 @@ def console_main(): yticks = yticks_corr, bound1=0.001, bound2=0.002, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, ax=axes[1], @@ -371,7 +388,7 @@ def console_main(): yticks = yticks_power, bound1=0.01, bound2=0.02, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, Ncrop=args.crop, @@ -392,7 +409,7 @@ def console_main(): yticks = yticks_corr, bound1=0.001, bound2=0.002, - kmin=kmin, + nfirst_kmodes=nfirst_kmodes, kmax=kmax, Nk=Nk, Ncrop=args.crop, From b4913dd9ffec5294a7aea5b19ff781759655adde Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Tue, 1 Jul 2025 15:15:58 +0200 Subject: [PATCH 6/6] bugfix convert snapshots to density --- .../scripts/convert_snapshot_to_density.py | 54 ++++++++++++++----- 1 file changed, 40 insertions(+), 14 deletions(-) diff --git a/sbmy_control/scripts/convert_snapshot_to_density.py b/sbmy_control/scripts/convert_snapshot_to_density.py index 869ee90..f3935d9 100644 --- a/sbmy_control/scripts/convert_snapshot_to_density.py +++ b/sbmy_control/scripts/convert_snapshot_to_density.py @@ -1,9 +1,7 @@ -from pysbmy.density import get_density_pm_snapshot -from pysbmy.snapshot import read_snapshot import argparse -def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, 0.0, 0.0)): +def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, 0.0, 0.0), time=1.0): """ Convert a snapshot to a density field. @@ -18,6 +16,9 @@ def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, corner : tuple of float Corner of the box (x, y, z). """ + from pysbmy.density import get_density_pm_snapshot + from pysbmy.snapshot import read_snapshot + # Read the snapshot print("Reading snapshot...") snap = read_snapshot(snapshot_path) @@ -28,6 +29,7 @@ def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, # Calculate density print("Calculating density...") F = get_density_pm_snapshot(snap, N, N, N, corner[0], corner[1], corner[2]) + F.time = time # Set the time for the field # Write density to file print("Writing density...") @@ -36,7 +38,7 @@ def convert_snapshot_to_density(snapshot_path, output_path, N=None, corner=(0.0, print("Done.") -def convert_snapshots_to_density(snapshot_base_path, output_path, N=None, corner=(0.0, 0.0, 0.0)): +def convert_snapshots_to_density(snapshot_base_path, output_path, N=None, corner=(0.0, 0.0, 0.0), time=1.0): """ Convert multiple snapshots to density fields. @@ -51,33 +53,57 @@ def convert_snapshots_to_density(snapshot_base_path, output_path, N=None, corner corner : tuple of float Corner of the box (x, y, z). """ + from pysbmy.density import get_density_pm, density_to_delta + from pysbmy.snapshot import read_snapshot + from pysbmy.field import Field + import numpy as np + # Get all snapshot files with path "snapshot_base_path*" import glob + snapshot_files = glob.glob(snapshot_base_path + "*") if not snapshot_files: raise FileNotFoundError(f"No snapshot files found at {snapshot_base_path}") A = None - F = None + L0 = L1 = L2 = None + d0 = d1 = d2 = None for snapshot_path in snapshot_files: # Read the snapshot print(f"Reading snapshot {snapshot_path}...") snap = read_snapshot(snapshot_path) - + if N is None: N = snap.Np0 - - # Calculate density - print("Calculating density...") - F = get_density_pm_snapshot(snap, N, N, N, corner[0], corner[1], corner[2]) + + if L0 is None: + # Get the grid parameters from the snapshot + L0 = snap.L0 + L1 = snap.L1 + L2 = snap.L2 + d0 = L0 / N + d1 = L1 / N + d2 = L2 / N + else: + # Ensure the grid parameters match the first snapshot + if not (snap.L0 == L0 and snap.L1 == L1 and snap.L2 == L2): + raise ValueError(f"All snapshots must have the same grid parameters. Got {snap.L0}, {snap.L1}, {snap.L2} but expected {L0}, {L1}, {L2}.") if A is None: - A = F.data - else: - A += F.data + # Initialize the density field + A = np.zeros((N, N, N), dtype=np.float32) + + # Calculate density + print("Calculating density for this snapshot...") + A = get_density_pm(snap.pos, A, d0=d0, d1=d1, d2=d2) - F.data = A # Combine densities from all snapshots + # Convert density to delta field + print("Converting density to delta field...") + A = density_to_delta(A) + + # Create a Field object + F = Field(L0=L0,L1=L1,L2=L2,corner0=corner[0],corner1=corner[1],corner2=corner[2],rank=1,N0=N,N1=N,N2=N,time=time,data=A) # Write density to file print("Writing density...") F.write(output_path)