From 928ea5eee86b35fe5c5911277150ebedef1601a5 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Fri, 14 Mar 2025 11:56:37 +0100 Subject: [PATCH 01/12] include parameters for phiBCs_from_coarse_tCOLA --- parameters_card.py | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/parameters_card.py b/parameters_card.py index 223b5ac..fa7a6c7 100644 --- a/parameters_card.py +++ b/parameters_card.py @@ -47,6 +47,13 @@ def register_arguments_card(parser:ArgumentParser): parser.add_argument("--OutputFCsDensity", type=str, default=None, help="Output FCs density file.") parser.add_argument("--OutputFCsSnapshot", type=str, default=None, help="Output FCs snapshot file.") parser.add_argument("--OutputRngStateLPT", type=str, default=None, help="Output RNG state file.") + ## Tests with phiBCs and density + parser.add_argument("--WriteGravPot", type=bool, default=True, help="Write gravitational potential.") + parser.add_argument("--OutputGravitationalPotentialBase", type=str, default=None, help="Output gravitational potential base.") + parser.add_argument("--MeshGravPot", type=int, default=None, help="Mesh for gravitational potential.") + parser.add_argument("--WriteDensity", type=bool, default=False, help="Write density.") + parser.add_argument("--OutputDensityBase", type=str, default=None, help="Output density base.") + parser.add_argument("--MeshDensity", type=int, default=None, help="Mesh for density.") def register_arguments_card_for_ICs(parser:ArgumentParser): @@ -118,6 +125,14 @@ def parse_arguments_card(parsed_args): OutputFCsDensity=parsed_args.OutputFCsDensity, OutputFCsSnapshot=parsed_args.OutputFCsSnapshot, OutputRngStateLPT=parsed_args.OutputRngStateLPT, + ## Tests with phiBCs and density + WriteGravPot=parsed_args.WriteGravPot, + OutputGravitationalPotentialBase=parsed_args.OutputGravitationalPotentialBase, + MeshGravPot=parsed_args.MeshGravPot, + WriteDensity=parsed_args.WriteDensity, + OutputDensityBase=parsed_args.OutputDensityBase, + MeshDensity=parsed_args.MeshDensity, + ## Cosmological parameters h=cosmo_dict["h"], Omega_m=cosmo_dict["Omega_m"], Omega_b=cosmo_dict["Omega_b"], @@ -203,6 +218,15 @@ def parse_arguments_card(parsed_args): card_dict["OutputFCsSnapshot"] = main_dict["resultdir"]+"final_particles_"+main_dict["simname"]+".gadget3" if card_dict["OutputRngStateLPT"] is None: card_dict["OutputRngStateLPT"] = main_dict["workdir"]+"rng_state.h5" + ## Tests with phiBCs and density + if card_dict["OutputGravitationalPotentialBase"] is None: + card_dict["OutputGravitationalPotentialBase"] = main_dict["workdir"]+"gravpot_"+main_dict["simname"] + if card_dict["MeshGravPot"] is None: + card_dict["MeshGravPot"] = card_dict["N_PM_mesh"] + if card_dict["OutputDensityBase"] is None: + card_dict["OutputDensityBase"] = main_dict["workdir"]+"density_"+main_dict["simname"] + if card_dict["MeshDensity"] is None: + card_dict["MeshDensity"] = card_dict["N_PM_mesh"] return card_dict @@ -321,6 +345,14 @@ def create_parameter_card_dict( OutputFCsDensity:str = 'fcs_density.h5', OutputFCsSnapshot:str = 'fcs_particles.gadget3', + ## Tests with phiBCs and density + WriteGravPot:bool = True, + OutputGravitationalPotentialBase:str = 'gravitational_potential.h5', + MeshGravPot:int = 128, + WriteDensity:bool = False, + OutputDensityBase:str = 'density.h5', + MeshDensity:int = 128, + ## Cosmological parameters h:float = 0.6732, Omega_m:float = 0.302, @@ -381,6 +413,15 @@ def create_parameter_card_dict( OutputLPTPotential1=OutputLPTPotential1, OutputLPTPotential2=OutputLPTPotential2, OutputTilesBase=OutputTilesBase, + + # Tests with phiBCs and density + WriteGravPot=int(WriteGravPot), + OutputGravitationalPotentialBase=OutputGravitationalPotentialBase, + MeshGravPot=MeshGravPot, + WriteDensity=int(WriteDensity), + OutputDensityBase=OutputDensityBase, + MeshDensity=MeshDensity, + h=h, Omega_m=Omega_m, Omega_b=Omega_b, From 9d6f8f9de20295ac079ee239dd4b8702222d5573 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 11:37:30 +0100 Subject: [PATCH 02/12] added analysis --- analysis/__init__.py | 0 analysis/power_spectrum.py | 335 +++++++++++++++++++++++++++++++++++++ analysis/slices.py | 114 +++++++++++++ 3 files changed, 449 insertions(+) create mode 100644 analysis/__init__.py create mode 100644 analysis/power_spectrum.py create mode 100644 analysis/slices.py diff --git a/analysis/__init__.py b/analysis/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/analysis/power_spectrum.py b/analysis/power_spectrum.py new file mode 100644 index 0000000..53fdc4f --- /dev/null +++ b/analysis/power_spectrum.py @@ -0,0 +1,335 @@ +import numpy as np + +kmin = 1e-1 +kmax = 2e0 +Nk = 50 +AliasingCorr=False + +def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk): + from pysbmy.power import PowerSpectrum + from pysbmy.fft import FourierGrid + from pysbmy.correlations import get_autocorrelation + + + G = FourierGrid( + field.L0, + field.L1, + field.L2, + field.N0, + field.N1, + field.N2, + k_modes=np.concat([PowerSpectrum(field.L0,field.L1,field.L2,field.N0,field.N1,field.N2,).FourierGrid.k_modes[:10],np.logspace( + np.log10(kmin), + np.log10(kmax), + Nk, + )]), + kmax=kmax, + ) + k = G.k_modes[1:] + Pk, _ = get_autocorrelation(field, G, AliasingCorr) + Pk = Pk[1:] + + return k, Pk + + +def get_cross_correlations(field_A, field_B, kmin=kmin, kmax=kmax, Nk=Nk): + from pysbmy.power import PowerSpectrum + from pysbmy.fft import FourierGrid + from pysbmy.correlations import get_crosscorrelation + + + G = FourierGrid( + field_A.L0, + field_A.L1, + field_A.L2, + field_A.N0, + field_A.N1, + field_A.N2, + k_modes=np.concat([PowerSpectrum(field_A.L0,field_A.L1,field_A.L2,field_A.N0,field_A.N1,field_A.N2,).FourierGrid.k_modes[:10],np.logspace( + np.log10(kmin), + np.log10(kmax), + Nk, + )]), + kmax=kmax, + ) + k = G.k_modes[1:] + _, _, Rks, _ = get_crosscorrelation(field_A, field_B, G, AliasingCorr) + Rks = Rks[1:] + + return k, Rks + + +def add_power_spectrum_to_plot(ax, field, Pk_ref=None, plot_args={}, power_args={}): + k, Pk = get_power_spectrum(field, **power_args) + if Pk_ref is not None: + ax.plot(k, Pk/Pk_ref-1, **plot_args) + else: + ax.plot(k, Pk, **plot_args) + return ax + +def add_cross_correlations_to_plot(ax, field_A, field_B, plot_args={}, power_args={}): + k, Rks = get_cross_correlations(field_A, field_B, **power_args) + ax.plot(k, Rks, **plot_args) + return ax + + +def plot_power_spectra(filenames, + labels=None, + colors=None, + linestyles=None, + markers=None, + Pk_ref=None, + ylims=[0.9,1.1], + yticks = np.linspace(0.9,1.1,11), + bound1=0.01, + bound2=0.02, + kmin=kmin, + kmax=kmax, + Nk=Nk, + figsize=(8,4), + dpi=300, + ax=None, + fig=None,): + + import matplotlib.pyplot as plt + from pysbmy.field import read_field + + if fig is None or ax is None: + fig, ax = plt.subplots(figsize=figsize, dpi=dpi) + + if labels is None: + labels = [None for f in filenames] + if colors is None: + colors = [None for f in filenames] + if linestyles is None: + linestyles = [None for f in filenames] + if markers is None: + markers = [None for f in filenames] + + for i, filename in enumerate(filenames): + field = read_field(filename) + add_power_spectrum_to_plot(ax=ax, + field=field, + Pk_ref=Pk_ref, + plot_args=dict(label=labels[i], + color=colors[i], + linestyle=linestyles[i], + marker=markers[i],), + power_args=dict(kmin=kmin, + kmax=kmax, + Nk=Nk), + ) + ax.set_xscale('log') + ax.set_xlim(kmin, kmax) + if ylims is not None: + ax.set_ylim(ylims) + if yticks is not None: + ax.set_yticks(yticks) + ax.set_xlabel(r'$k$ [$h/\mathrm{Mpc}$]') + + if Pk_ref is not None: + ax.set_ylabel(r'$P(k)/P_\mathrm{ref}(k)-1$') + else: + ax.set_ylabel('$P(k)$') + + if bound1 is not None: + ax.axhspan(1 - bound1, 1 + bound1, color="grey", alpha=0.2) + if bound2 is not None: + ax.axhspan(1 - bound2, 1 + bound2, color="grey", alpha=0.1) + + ax.grid(which="major",alpha=0.5) + ax.grid(which="minor",alpha=0.2) + + return fig, ax + + +def plot_cross_correlations(filenames_A, + filename_B, + labels=None, + colors=None, + linestyles=None, + markers=None, + ylims=[0.99, 1.001], + yticks = np.linspace(0.99,1.001,12), + bound1=0.001, + bound2=0.002, + kmin=kmin, + kmax=kmax, + Nk=Nk, + figsize=(8,4), + dpi=300, + ax=None, + fig=None,): + + import matplotlib.pyplot as plt + from pysbmy.field import read_field + + if fig is None or ax is None: + fig, ax = plt.subplots(figsize=figsize, dpi=dpi) + + if labels is None: + labels = [None for f in filenames_A] + if colors is None: + colors = [None for f in filenames_A] + if linestyles is None: + linestyles = [None for f in filenames_A] + if markers is None: + markers = [None for f in filenames_A] + + field_B = read_field(filename_B) + + for i, filename_A in enumerate(filenames_A): + field_A = read_field(filename_A) + add_cross_correlations_to_plot(ax=ax, + field_A=field_A, + field_B=field_B, + plot_args=dict(label=labels[i], + color=colors[i], + linestyle=linestyles[i], + marker=markers[i],), + power_args=dict(kmin=kmin, + kmax=kmax, + Nk=Nk), + ) + ax.set_xscale('log') + ax.set_xlim(kmin, kmax) + if ylims is not None: + ax.set_ylim(ylims) + if yticks is not None: + ax.set_yticks(yticks) + ax.set_xlabel(r'$k$ [$h/\mathrm{Mpc}$]') + ax.set_ylabel('$R(k)$') + + if bound1 is not None: + ax.axhspan(1 - bound1, 1 + bound1, color="grey", alpha=0.2) + if bound2 is not None: + ax.axhspan(1 - bound2, 1 + bound2, color="grey", alpha=0.1) + + ax.grid(which="major",alpha=0.5) + ax.grid(which="minor",alpha=0.2) + + return fig, ax + + +if __name__ == "__main__": + from argparse import ArgumentParser + parser = ArgumentParser(description='Plot power spectra of fields') + + parser.add_argument('-ps', '--power_spectrum', action='store_true', help='Plot power spectra.') + parser.add_argument('-cc', '--cross_correlation', action='store_true', help='Plot cross correlations.') + parser.add_argument('-d', '--directory', type=str, required=True, help='Directory containing the fields files.') + parser.add_argument('-ref', '--reference', type=str, default=None, help='Reference field file.') + parser.add_argument('-f', '--filenames', type=str, nargs='+', required=True, help='Field files to be plotted.') + parser.add_argument('-o', '--output', type=str, default=None, help='Output plot file name.') + parser.add_argument('-l', '--labels', type=str, nargs='+', default=None, help='Labels for each field.') + parser.add_argument('-c', '--colors', type=str, nargs='+', default=None, help='Colors for each field.') + parser.add_argument('-ls', '--linestyles', type=str, nargs='+', default=None, help='Linestyles for each field.') + parser.add_argument('-m', '--markers', type=str, nargs='+', default=None, help='Markers for each field.') + parser.add_argument('-t','--title', type=str, default=None, help='Title of the plot.') + + args = parser.parse_args() + + if not args.power_spectrum and not args.cross_correlation: + print('You must choose between power_spectrum and cross_correlation.') + exit() + + if args.reference is not None: + from pysbmy.field import read_field + Pk_ref = get_power_spectrum(read_field(args.directory+args.reference), kmin=kmin, kmax=kmax, Nk=Nk)[1] + else: + Pk_ref = None + + filenames = [args.directory+f for f in args.filenames] + + if args.power_spectrum and args.cross_correlation: + import matplotlib.pyplot as plt + fig, axes = plt.subplots(2, 1, figsize=(8,8)) + plot_power_spectra(filenames=filenames, + labels=args.labels, + colors=args.colors, + linestyles=args.linestyles, + markers=args.markers, + Pk_ref=Pk_ref, + # ylims=[0.9,1.1], + # yticks = np.linspace(0.9,1.1,11), + # bound1=0.01, + # bound2=0.02, + ylims=None, + yticks = None, + bound1=None, + bound2=None, + kmin=kmin, + kmax=kmax, + Nk=Nk, + ax=axes[0], + fig=fig) + + plot_cross_correlations(filenames_A=filenames, + filename_B=args.directory+args.reference, + labels=args.labels, + colors=args.colors, + linestyles=args.linestyles, + markers=args.markers, + # ylims=[0.99, 1.001], + # yticks = np.linspace(0.99,1.001,12), + # bound1=0.001, + # bound2=0.002, + ylims=None, + yticks = None, + bound1=None, + bound2=None, + kmin=kmin, + kmax=kmax, + Nk=Nk, + ax=axes[1], + fig=fig) + + axes[1].legend(loc='lower left') + axes[0].set_title("Power Spectrum") + axes[1].set_title("Cross Correlations") + + if args.title is not None: + fig.suptitle(args.title) + + elif args.power_spectrum: + fig, ax = plot_power_spectra(filenames=filenames, + labels=args.labels, + colors=args.colors, + linestyles=args.linestyles, + markers=args.markers, + Pk_ref=Pk_ref, + ylims=[0.9,1.1], + yticks = np.linspace(0.9,1.1,11), + bound1=0.01, + bound2=0.02, + kmin=kmin, + kmax=kmax, + Nk=Nk) + ax.legend() + if args.title is not None: + ax.set_title(args.title) + + elif args.cross_correlation: + fig, ax = plot_cross_correlations(filenames_A=filenames, + filename_B=args.reference, + labels=args.labels, + colors=args.colors, + linestyles=args.linestyles, + markers=args.markers, + ylims=[0.99, 1.001], + yticks = np.linspace(0.99,1.001,12), + bound1=0.001, + bound2=0.002, + kmin=kmin, + kmax=kmax, + Nk=Nk) + ax.legend(loc='lower left') + if args.title is not None: + ax.set_title(args.title) + + if args.output is not None: + fig.savefig(args.output) + else: + fig.savefig(args.directory+'power_spectrum.png') + + diff --git a/analysis/slices.py b/analysis/slices.py new file mode 100644 index 0000000..0456e7b --- /dev/null +++ b/analysis/slices.py @@ -0,0 +1,114 @@ +import numpy as np + +fs = 18 +fs_titles = fs -4 + +def plot_imshow_with_reference( data_list, + reference, + titles, + vmin=None, + vmax=None, + cmap='viridis'): + """ + Plot the imshow of a list of 2D arrays with two rows: one for the data itself, + one for the data compared to a reference. Each row will have a common colorbar. + + Parameters: + - data_list: list of 2D arrays to be plotted + - reference: 2D array to be used as reference for comparison + - titles: list of titles for each subplot + - cmap: colormap to be used for plotting + """ + import matplotlib.pyplot as plt + + if titles is None: + titles = [None for f in data_list] + + def score(data, reference): + return np.linalg.norm(data-reference)/np.linalg.norm(reference) + + n = len(data_list) + fig, axes = plt.subplots(2, n, figsize=(5 * n, 10)) + + if vmin is None or vmax is None: + vmin = min(np.quantile(data,0.01) for data in data_list) + vmax = max(np.quantile(data,0.99) for data in data_list) + vmin_diff = min(np.quantile((data-reference),0.01) for data in data_list) + vmax_diff = max(np.quantile((data-reference),0.99) for data in data_list) + else: + vmin_diff = vmin + vmax_diff = vmax + + # Plot the data itself + for i, data in enumerate(data_list): + im = axes[0, i].imshow(data, cmap=cmap, origin='lower', vmin=vmin, vmax=vmax) + axes[0, i].set_title(titles[i], fontsize=fs_titles) + fig.colorbar(im, ax=axes[0, :], orientation='vertical') + + # Plot the data compared to the reference + for i, data in enumerate(data_list): + im = axes[1, i].imshow(data - reference, cmap=cmap, origin='lower', vmin=vmin_diff, vmax=vmax_diff) + axes[1, i].set_title(f'{titles[i]} - Reference', fontsize=fs_titles) + fig.colorbar(im, ax=axes[1, :], orientation='vertical') + + # Add the score on the plots + for i, data in enumerate(data_list): + axes[1, i].text(0.5, 0.9, f"Score: {score(data, reference):.2e}", fontsize=10, transform=axes[1, i].transAxes, color='white') + # plt.tight_layout() + + return fig, axes + + +if __name__ == "__main__": + from argparse import ArgumentParser + parser = ArgumentParser(description='Comparisons of fields slices.') + + parser.add_argument('-a','--axis', type=int, default=0, help='Axis along which the slices will be taken.') + parser.add_argument('-i','--index', type=int, default=None, help='Index of the slice along the axis.') + parser.add_argument('-d', '--directory', type=str, required=True, help='Directory containing the fields files.') + parser.add_argument('-ref', '--reference', type=str, default=None, help='Reference field file.') + parser.add_argument('-f', '--filenames', type=str, nargs='+', required=True, help='Field files to be plotted.') + parser.add_argument('-o', '--output', type=str, default=None, help='Output plot file name.') + parser.add_argument('-l', '--labels', type=str, nargs='+', default=None, help='Labels for each field.') + parser.add_argument('-c', '--cmap', type=str, default='viridis', help='Colormap to be used for plotting.') + parser.add_argument('-vmin', type=float, default=None, help='Minimum value for the colorbar.') + parser.add_argument('-vmax', type=float, default=None, help='Maximum value for the colorbar.') + parser.add_argument('-t', '--title', type=str, default=None, help='Title for the plot.') + parser.add_argument('-log','--log_scale', action='store_true', help='Use log scale for the data.') + + args = parser.parse_args() + + from pysbmy.field import read_field + + ref_field = read_field(args.directory+args.reference) + fields = [read_field(args.directory+f) for f in args.filenames] + + if args.index is None: + index = ref_field.N0//2 + else: + index=args.index + + match args.axis: + case 0 | 'x': + reference = ref_field.data[index,:,:] + fields = [f.data[index,:,:] for f in fields] + case 1 | 'y': + reference = ref_field.data[:,index,:] + fields = [f.data[:,index,:] for f in fields] + case 2 | 'z': + reference = ref_field.data[:,:,index] + fields = [f.data[:,:,index] for f in fields] + case _: + raise ValueError(f"Wrong axis provided : {args.axis}") + + if args.log_scale: + reference = np.log10(2.+reference) + fields = [np.log10(2.+f) for f in fields] + + fig, axes = plot_imshow_with_reference(fields,reference,args.labels, vmin=args.vmin, vmax=args.vmax,cmap=args.cmap) + fig.suptitle(args.title) + + if args.output is not None: + fig.savefig(args.output) + else: + fig.savefig(args.directory+'slices.png') \ No newline at end of file From b0c106bd8bd7ee369fbf70720dd4dd1a94dac915 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 11:37:48 +0100 Subject: [PATCH 03/12] value error unknown mode --- main.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/main.py b/main.py index ff18d4f..54be5a7 100644 --- a/main.py +++ b/main.py @@ -190,7 +190,9 @@ def main(parsed_args): print_message("Running allsCOLA finished.", 1, "control", verbose=parsed_args.verbose) - + case _: + raise ValueError(f"Unknown mode: {main_dict['mode']}") + print_ending_module("control", verbose=parsed_args.verbose) From 78e1b4cacf500cef0488b5dbcfcbd2625f042898 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 11:38:04 +0100 Subject: [PATCH 04/12] scola output --- scola.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scola.py b/scola.py index 3d1af5c..87e2f43 100644 --- a/scola.py +++ b/scola.py @@ -52,7 +52,7 @@ def main_scola(parsed_args): if have_no_tiles(parsed_args) or parsed_args.force: ## Submit all boxes print_message("Submitting all boxes.", 2, "scola", verbose=parsed_args.verbose) - slurm_script = slurm_dict["scripts"]+"scola.sh" + slurm_script = slurm_dict["scripts"]+"scola_"+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, "scola", verbose=parsed_args.verbose) @@ -85,7 +85,7 @@ def main_scola(parsed_args): print_message(f"Submitting missing boxes: {missing_tiles_arrays}.", 2, "scola", verbose=parsed_args.verbose) for missing_tiles in missing_tiles_arrays: - slurm_script = slurm_dict["scripts"]+"scola_"+str(missing_tiles[0])+"_"+str(missing_tiles[-1])+".sh" + slurm_script = slurm_dict["scripts"]+"scola_"+main_dict["simname"]+"_"+str(missing_tiles[0])+"_"+str(missing_tiles[-1])+".sh" if not isfile(slurm_script): print_message(f"SLURM script {slurm_script} does not exist. Creating it.", 2, "scola", verbose=parsed_args.verbose) From 1ca9b4fb09530b6450c1a1be619c983be8828083 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 11:38:30 +0100 Subject: [PATCH 05/12] grav pot output --- parameters_card.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/parameters_card.py b/parameters_card.py index fa7a6c7..fe1768c 100644 --- a/parameters_card.py +++ b/parameters_card.py @@ -48,12 +48,14 @@ def register_arguments_card(parser:ArgumentParser): parser.add_argument("--OutputFCsSnapshot", type=str, default=None, help="Output FCs snapshot file.") parser.add_argument("--OutputRngStateLPT", type=str, default=None, help="Output RNG state file.") ## Tests with phiBCs and density - parser.add_argument("--WriteGravPot", type=bool, default=True, help="Write gravitational potential.") + parser.add_argument("--WriteGravPot", type=bool, default=False, help="Write gravitational potential.") parser.add_argument("--OutputGravitationalPotentialBase", type=str, default=None, help="Output gravitational potential base.") parser.add_argument("--MeshGravPot", type=int, default=None, help="Mesh for gravitational potential.") parser.add_argument("--WriteDensity", type=bool, default=False, help="Write density.") parser.add_argument("--OutputDensityBase", type=str, default=None, help="Output density base.") parser.add_argument("--MeshDensity", type=int, default=None, help="Mesh for density.") + parser.add_argument("--LoadPhiBCs", type=bool, default=False, help="Load phiBCs.") + parser.add_argument("--InputPhiBCsBase", type=str, default=None, help="Input phiBCs file base.") def register_arguments_card_for_ICs(parser:ArgumentParser): @@ -132,6 +134,8 @@ def parse_arguments_card(parsed_args): WriteDensity=parsed_args.WriteDensity, OutputDensityBase=parsed_args.OutputDensityBase, MeshDensity=parsed_args.MeshDensity, + LoadPhiBCs=parsed_args.LoadPhiBCs, + InputPhiBCsBase=parsed_args.InputPhiBCsBase, ## Cosmological parameters h=cosmo_dict["h"], Omega_m=cosmo_dict["Omega_m"], @@ -207,7 +211,7 @@ def parse_arguments_card(parsed_args): if card_dict["OutputFinalDensity"] is None: card_dict["OutputFinalDensity"] = main_dict["resultdir"]+ligthcone_prefix+"final_density_"+main_dict["simname"]+".h5" if card_dict["OutputTilesBase"] is None: - card_dict["OutputTilesBase"] = main_dict["workdir"]+"sCOLA_tile" + card_dict["OutputTilesBase"] = main_dict["workdir"]+main_dict["simname"]+"_tile" if card_dict["OutputLPTPotential1"] is None: card_dict["OutputLPTPotential1"] = main_dict["workdir"]+"initial_conditions_DM_phi.h5" if card_dict["OutputLPTPotential2"] is None: @@ -227,6 +231,8 @@ def parse_arguments_card(parsed_args): card_dict["OutputDensityBase"] = main_dict["workdir"]+"density_"+main_dict["simname"] if card_dict["MeshDensity"] is None: card_dict["MeshDensity"] = card_dict["N_PM_mesh"] + if card_dict["InputPhiBCsBase"] is None: + card_dict["InputPhiBCsBase"] = main_dict["workdir"]+"gravpot_tCOLA" return card_dict @@ -352,6 +358,8 @@ def create_parameter_card_dict( WriteDensity:bool = False, OutputDensityBase:str = 'density.h5', MeshDensity:int = 128, + LoadPhiBCs:bool = False, + InputPhiBCsBase:str = 'gravitational_potential.h5', ## Cosmological parameters h:float = 0.6732, @@ -421,6 +429,8 @@ def create_parameter_card_dict( WriteDensity=int(WriteDensity), OutputDensityBase=OutputDensityBase, MeshDensity=MeshDensity, + LoadPhiBCs=int(LoadPhiBCs), + InputPhiBCsBase=InputPhiBCsBase, h=h, Omega_m=Omega_m, From fdaf5e90868ef13173a0b661f554b9deafea2de0 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 13:54:54 +0100 Subject: [PATCH 06/12] fixing progress bars with log files --- low_level.py | 4 ++++ scola.py | 11 +++++++++++ simbelmyne.py | 10 ++++++++-- 3 files changed, 23 insertions(+), 2 deletions(-) diff --git a/low_level.py b/low_level.py index e721508..e987cf0 100644 --- a/low_level.py +++ b/low_level.py @@ -240,6 +240,10 @@ def progress_bar_from_logfile(filename:str, desc:str="", verbose:int=1, **kwargs current_operation, total_operations = get_progress_from_logfile(filename) previous_operation = 0 + while total_operations == 0: # Wait for the simulation to be initialised, and to run the first operation + sleep(2) + current_operation, total_operations = get_progress_from_logfile(filename) + if current_operation == total_operations: # print_message("Progress bar not needed, the process is already finished.", 3, "low level", verbose=verbose) return diff --git a/scola.py b/scola.py index 87e2f43..a6156ad 100644 --- a/scola.py +++ b/scola.py @@ -22,6 +22,9 @@ def main_scola(parsed_args): tile_bar = tqdm(total=nboxes_tot, desc="sCOLA", unit="box", disable=(parsed_args.verbose!=1)) # The progress bar cannot work if there are print statements in the loop for b in range(1,nboxes_tot+1): if not isfile(card_dict["OutputTilesBase"]+str(b)+".h5") or parsed_args.force: + if isfile(log_file+f"_{b}"): # Remove the preexisting log file to allow for the progress_bar to be run normally + from os import remove + remove(log_file+f"_{b}") if parsed_args.verbose > 1: print_message(f"Running box {b}/{nboxes_tot}.", 2, "scola", verbose=parsed_args.verbose) command_args = ["scola", paramfile, log_file, "-b", str(b)] @@ -72,6 +75,10 @@ def main_scola(parsed_args): print_message(f"Submitting scola job to SLURM.", 1, "scola", verbose=parsed_args.verbose) command_args = ["sbatch", slurm_script] + for b in range(1,nboxes_tot+1): + if isfile(log_file+f"_{b}"): # Remove the preexisting log file to allow for the progress_bar to be run normally + os.remove(log_file+f"_{b}") + if parsed_args.verbose < 2: subprocess.run(command_args, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) else: @@ -105,6 +112,10 @@ def main_scola(parsed_args): print_message(f"Submitting scola job to SLURM.", 1, "scola", verbose=parsed_args.verbose) command_args = ["sbatch", slurm_script] + for b in range(missing_tiles[0],missing_tiles[1]+1): + if isfile(log_file+f"_{b}"): # Remove the preexisting log file to allow for the progress_bar to be run normally + os.remove(log_file+f"_{b}") + if parsed_args.verbose < 2: subprocess.run(command_args, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) else: diff --git a/simbelmyne.py b/simbelmyne.py index 417845d..3ef9a3f 100644 --- a/simbelmyne.py +++ b/simbelmyne.py @@ -14,7 +14,10 @@ def main_simbelmyne(parsed_args): if parsed_args.execution == "local": print_message("Running Simbelyne in local mode.", 1, "simbelmyne", verbose=parsed_args.verbose) - + 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) + command_args = ["simbelmyne", paramfile, log_file] if parsed_args.verbose < 2: @@ -51,6 +54,10 @@ def main_simbelmyne(parsed_args): command_args = ["sbatch", slurm_script] + 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) + if parsed_args.verbose < 2: subprocess.run(command_args, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) else: @@ -75,7 +82,6 @@ if __name__ == "__main__": from slurm_submission import register_arguments_slurm parser = ArgumentParser(description="Run Simbelmyne.") - # TODO: reduce the volume of arguments register_arguments_main(parser) register_arguments_timestepping(parser) register_arguments_simbelmyne(parser) From bd0e5f9ffc2076dd2cef50419e0b51b557c0c539 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 16:36:09 +0100 Subject: [PATCH 07/12] stop progress bar if error is found --- low_level.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/low_level.py b/low_level.py index e987cf0..8d56a98 100644 --- a/low_level.py +++ b/low_level.py @@ -223,6 +223,8 @@ def get_progress_from_logfile(filename): total_operations = int(splitted_line.split("/")[1]) except: pass + elif "Fatal" in line or "Error" in line: + return -1, -1 return current_operation, total_operations @@ -246,6 +248,8 @@ def progress_bar_from_logfile(filename:str, desc:str="", verbose:int=1, **kwargs if current_operation == total_operations: # print_message("Progress bar not needed, the process is already finished.", 3, "low level", verbose=verbose) + if current_operation == -1: + print_message("Error appeared in log file, skipping it.",level=3,module="low level",verbose=verbose) return with tqdm(desc=desc, total=total_operations, disable=(verbose==0), **kwargs) as pbar: @@ -255,6 +259,11 @@ def progress_bar_from_logfile(filename:str, desc:str="", verbose:int=1, **kwargs if current_operation > previous_operation: pbar.update(current_operation-previous_operation) previous_operation = current_operation + if current_operation == -1: + print_message("Error appeared in log file, skipping it.",level=3,module="low level",verbose=verbose) + return k+=1 if k*update_interval >= limit: - print_message(f"Progress bar timed out after {limit} seconds.", 3, "low level", verbose=verbose) \ No newline at end of file + print_message(f"Progress bar timed out after {limit} seconds.", 3, "low level", verbose=verbose) + + return \ No newline at end of file From a12e1a5ab48ca50691b8f1cec8d54ab0a8bd6b20 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 16:36:26 +0100 Subject: [PATCH 08/12] no progress bar for sbmy when scola --- simbelmyne.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/simbelmyne.py b/simbelmyne.py index 3ef9a3f..effb375 100644 --- a/simbelmyne.py +++ b/simbelmyne.py @@ -22,7 +22,8 @@ def main_simbelmyne(parsed_args): if parsed_args.verbose < 2: subprocess.Popen(command_args, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) # Use Popen instead of run to be able to display the progress bar - progress_bar_from_logfile(log_file, desc=main_dict["simname"], verbose=parsed_args.verbose, leave=False) + if not main_dict["mode"] in ["pre_sCOLA", "post_sCOLA", "allsCOLA"]: # No progress bar for pre/post sCOLA + progress_bar_from_logfile(log_file, desc=main_dict["simname"], verbose=parsed_args.verbose, leave=False) else: subprocess.run(command_args) @@ -56,7 +57,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) if parsed_args.verbose < 2: subprocess.run(command_args, stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) @@ -64,7 +65,8 @@ def main_simbelmyne(parsed_args): subprocess.run(command_args) print_message("Simbelmyne job submitted.", 2, "simbelmyne", verbose=parsed_args.verbose) - progress_bar_from_logfile(log_file, desc=main_dict["simname"], verbose=parsed_args.verbose) + if not main_dict["mode"] in ["pre_sCOLA", "post_sCOLA", "allsCOLA"]: # No progress bar for pre/post sCOLA + progress_bar_from_logfile(log_file, desc=main_dict["simname"], verbose=parsed_args.verbose) else: raise ValueError(f"Execution mode {parsed_args.execution} not recognized.") From 21723980c2c8c063b72f0635569b0ff457bf6669 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Wed, 9 Apr 2025 14:45:49 +0200 Subject: [PATCH 09/12] fixed fnl=0 --- cosmo_params.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cosmo_params.py b/cosmo_params.py index 14e48dc..d973866 100644 --- a/cosmo_params.py +++ b/cosmo_params.py @@ -19,7 +19,7 @@ def register_arguments_cosmo(parser:ArgumentParser): parser.add_argument("--m_nu3", type=float, default=0.0, help="Mass of the third neutrino species.") parser.add_argument("--w_0", type=float, default=-1.0, help="Dark energy equation of state parameter.") parser.add_argument("--w_a", type=float, default=0.0, help="Dark energy equation of state parameter.") - parser.add_argument("--fnl", type=float, default=100.0, help="Local non-Gaussianity parameter.") + parser.add_argument("--fnl", type=float, default=0.0, help="Local non-Gaussianity parameter.") parser.add_argument("--gnl", type=float, default=0.0, help="Equilateral non-Gaussianity parameter.") @@ -84,7 +84,7 @@ cosmo_defaults = { "m_nu3": 0.0, "w_0": -1.0, "w_a": 0.0, - "fnl": 100.0, + "fnl": 0.0, "gnl": 0.0, "k_max":10.0, "tau_reio":0.06, From 78d747d5aa265806a62b5625e4d9375785d3628e Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Wed, 9 Apr 2025 14:46:15 +0200 Subject: [PATCH 10/12] minor changes --- analysis/slices.py | 31 ++++++++++++++++++++++--------- main.py | 4 ++-- monofonic.py | 1 + parameters_monofonic.py | 1 - 4 files changed, 25 insertions(+), 12 deletions(-) diff --git a/analysis/slices.py b/analysis/slices.py index 8b260cc..b3fa8e1 100644 --- a/analysis/slices.py +++ b/analysis/slices.py @@ -83,15 +83,28 @@ def plot_imshow_with_reference( data_list, axes[1, i].text(0.5, 0.9, f"Score: {score(data, reference):.2e}", fontsize=10, transform=axes[1, i].transAxes, color='white') # plt.tight_layout() else: - for i, data in enumerate(data_list): - im = axes[i].imshow(data, cmap=cmap, origin='lower', vmin=vmin, vmax=vmax) - axes[i].set_title(titles[i], fontsize=fs_titles) - axes[i].set_xticks(ticks[i]) - axes[i].set_yticks(ticks[i]) - axes[i].set_xticklabels(tick_labels[i]) - axes[i].set_yticklabels(tick_labels[i]) - axes[i].set_xlabel('Mpc/h') - fig.colorbar(im, ax=axes[:], orientation='vertical') + + if len(data_list) == 1: + data_list = data_list[0] + im = axes.imshow(data_list, cmap=cmap, origin='lower', vmin=vmin, vmax=vmax) + axes.set_title(titles[0], fontsize=fs_titles) + axes.set_xticks(ticks[0]) + axes.set_yticks(ticks[0]) + axes.set_xticklabels(tick_labels[0]) + axes.set_yticklabels(tick_labels[0]) + axes.set_xlabel('Mpc/h') + fig.colorbar(im, ax=axes, orientation='vertical') + + else: + for i, data in enumerate(data_list): + im = axes[i].imshow(data, cmap=cmap, origin='lower', vmin=vmin, vmax=vmax) + axes[i].set_title(titles[i], fontsize=fs_titles) + axes[i].set_xticks(ticks[i]) + axes[i].set_yticks(ticks[i]) + axes[i].set_xticklabels(tick_labels[i]) + axes[i].set_yticklabels(tick_labels[i]) + axes[i].set_xlabel('Mpc/h') + fig.colorbar(im, ax=axes[:], orientation='vertical') return fig, axes diff --git a/main.py b/main.py index 54be5a7..357f93a 100644 --- a/main.py +++ b/main.py @@ -136,7 +136,7 @@ def main(parsed_args): from parameters_card import parse_arguments_card, main_parameter_card from timestepping import main_timestepping from ICs import main_ICs - from scola import main_scola, main_pre_scola, main_post_scola + from scola import main_scola, main_pre_scola, main_post_scola, have_all_tiles card_dict = main_parameter_card(parsed_args) @@ -177,7 +177,7 @@ def main(parsed_args): print_message("Running sCOLA.", 1, "control", verbose=parsed_args.verbose) main_scola(parsed_args) - if parsed_args.execution == "slurm": + if parsed_args.execution == "slurm" and not have_all_tiles(parsed_args): from tqdm import tqdm from low_level import progress_bar_from_logfile for b in tqdm(range(1,parsed_args.N_tiles**3+1), desc="sCOLA", unit="box", disable=(parsed_args.verbose==0)): diff --git a/monofonic.py b/monofonic.py index 6c0da8e..4fe11ee 100644 --- a/monofonic.py +++ b/monofonic.py @@ -9,6 +9,7 @@ def main_monofonic(parsed_args): print_starting_module("monofonic", verbose=parsed_args.verbose) monofonic_dict=main_parameters_monofonic(parsed_args) + # raise SystemExit # Exit the program after writing the config file to avoid running monofonic if parsed_args.execution == "local": print_message("Running monofonic in local mode.", 1, "monofonic", verbose=parsed_args.verbose) diff --git a/parameters_monofonic.py b/parameters_monofonic.py index 5554517..ead1055 100644 --- a/parameters_monofonic.py +++ b/parameters_monofonic.py @@ -140,7 +140,6 @@ def main_parameters_monofonic(parsed_args): return monofonic_dict create_monofonic_config(monofonic_dict["config"], get_config_from_dict(monofonic_dict)) print_message(f"Configuration file written to {monofonic_dict["config"]}", 2, "monofonic", verbose=parsed_args.verbose) - return monofonic_dict From 7d142b35b56ecf9280452aa4a7908026af7fdd4f Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 2 Jun 2025 18:18:53 +0200 Subject: [PATCH 11/12] merge fixes --- sbmy_control/main.py | 7 ------- sbmy_control/parameters_card.py | 3 --- 2 files changed, 10 deletions(-) diff --git a/sbmy_control/main.py b/sbmy_control/main.py index 350c97d..8bb7ff0 100644 --- a/sbmy_control/main.py +++ b/sbmy_control/main.py @@ -133,17 +133,10 @@ def main(parsed_args): case "allsCOLA": print_message(f"Running ICs, pre_sCOLA, sCOLA and post_sCOLA.", 1, "control", verbose=parsed_args.verbose) -<<<<<<< HEAD:main.py - from parameters_card import parse_arguments_card, main_parameter_card - from timestepping import main_timestepping - from ICs import main_ICs - from scola import main_scola, main_pre_scola, main_post_scola, have_all_tiles -======= from .parameters_card import parse_arguments_card, main_parameter_card from .timestepping import main_timestepping from .ICs import main_ICs from .scola import main_scola, main_pre_scola, main_post_scola, have_all_tiles ->>>>>>> main:sbmy_control/main.py card_dict = main_parameter_card(parsed_args) diff --git a/sbmy_control/parameters_card.py b/sbmy_control/parameters_card.py index 7381759..27a776f 100644 --- a/sbmy_control/parameters_card.py +++ b/sbmy_control/parameters_card.py @@ -47,7 +47,6 @@ def register_arguments_card(parser:ArgumentParser): parser.add_argument("--OutputFCsDensity", type=str, default=None, help="Output FCs density file.") parser.add_argument("--OutputFCsSnapshot", type=str, default=None, help="Output FCs snapshot file.") parser.add_argument("--OutputRngStateLPT", type=str, default=None, help="Output RNG state file.") -<<<<<<< HEAD:parameters_card.py ## Tests with phiBCs and density parser.add_argument("--WriteGravPot", type=bool, default=False, help="Write gravitational potential.") parser.add_argument("--OutputGravitationalPotentialBase", type=str, default=None, help="Output gravitational potential base.") @@ -57,12 +56,10 @@ def register_arguments_card(parser:ArgumentParser): parser.add_argument("--MeshDensity", type=int, default=None, help="Mesh for density.") parser.add_argument("--LoadPhiBCs", type=bool, default=False, help="Load phiBCs.") parser.add_argument("--InputPhiBCsBase", type=str, default=None, help="Input phiBCs file base.") -======= parser.add_argument("--WriteReferenceFrame", type=bool, default=False, help="Write reference frame (COCA).") parser.add_argument("--OutputMomentaBase", type=str, default=None, help="Output momenta base (COCA).") parser.add_argument("--ReadReferenceFrame", type=bool, default=False, help="Read reference frame (COCA).") parser.add_argument("--InputMomentaBase", type=str, default=None, help="Read momenta base (COCA).") ->>>>>>> main:sbmy_control/parameters_card.py def register_arguments_card_for_ICs(parser:ArgumentParser): From 8c92adcbcf4fb964e568c50f132634e6c0b62c75 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 2 Jun 2025 18:19:32 +0200 Subject: [PATCH 12/12] merge fixes --- sbmy_control/parameters_card.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/sbmy_control/parameters_card.py b/sbmy_control/parameters_card.py index 27a776f..5384adb 100644 --- a/sbmy_control/parameters_card.py +++ b/sbmy_control/parameters_card.py @@ -230,7 +230,6 @@ def parse_arguments_card(parsed_args): card_dict["OutputFCsSnapshot"] = main_dict["resultdir"]+"final_particles_"+main_dict["simname"]+".gadget3" if card_dict["OutputRngStateLPT"] is None: card_dict["OutputRngStateLPT"] = main_dict["workdir"]+"rng_state.h5" -<<<<<<< HEAD:parameters_card.py ## Tests with phiBCs and density if card_dict["OutputGravitationalPotentialBase"] is None: card_dict["OutputGravitationalPotentialBase"] = main_dict["workdir"]+"gravpot_"+main_dict["simname"] @@ -242,12 +241,10 @@ def parse_arguments_card(parsed_args): card_dict["MeshDensity"] = card_dict["N_PM_mesh"] if card_dict["InputPhiBCsBase"] is None: card_dict["InputPhiBCsBase"] = main_dict["workdir"]+"gravpot_tCOLA" -======= if card_dict["OutputMomentaBase"] is None: card_dict["OutputMomentaBase"] = main_dict["workdir"]+"momenta_"+main_dict["simname"]+"_" if card_dict["InputMomentaBase"] is None: card_dict["InputMomentaBase"] = main_dict["workdir"]+"momenta_"+main_dict["simname"]+"_" ->>>>>>> main:sbmy_control/parameters_card.py return card_dict @@ -366,7 +363,6 @@ def create_parameter_card_dict( OutputFCsDensity:str = 'fcs_density.h5', OutputFCsSnapshot:str = 'fcs_particles.gadget3', -<<<<<<< HEAD:parameters_card.py ## Tests with phiBCs and density WriteGravPot:bool = True, OutputGravitationalPotentialBase:str = 'gravitational_potential.h5', @@ -376,13 +372,11 @@ def create_parameter_card_dict( MeshDensity:int = 128, LoadPhiBCs:bool = False, InputPhiBCsBase:str = 'gravitational_potential.h5', -======= ## COCA parameters WriteReferenceFrame:bool = False, OutputMomentaBase:str = 'momenta_', ReadReferenceFrame:bool = False, InputMomentaBase:str = 'momenta_', ->>>>>>> main:sbmy_control/parameters_card.py ## Cosmological parameters h:float = 0.6732,