From 72ce9a3b99f1f61dfca9b93336e239877aae57b2 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 11:37:30 +0100 Subject: [PATCH 01/13] 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 f52a6e39573eebf5167a0d01f4de0ea97f49a3e1 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 11:37:48 +0100 Subject: [PATCH 02/13] 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 de414576b8cd9470a000d5184a3233456ed18005 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 11:38:04 +0100 Subject: [PATCH 03/13] 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 32fa9a5afb0dff95e98bed0c2aa9b79ffb394527 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 11:39:46 +0100 Subject: [PATCH 04/13] scola output tile --- parameters_card.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parameters_card.py b/parameters_card.py index 223b5ac..df03911 100644 --- a/parameters_card.py +++ b/parameters_card.py @@ -192,7 +192,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: From 942f775425f813fc975368b6dc77719d950fe45c Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 13:54:54 +0100 Subject: [PATCH 05/13] 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 c80c020e4f27fae0cd0db0e4f3bb18e300b67092 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 16:36:09 +0100 Subject: [PATCH 06/13] 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 536d3df365a62c59581f618bb346be9232a45f39 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 17 Mar 2025 16:36:26 +0100 Subject: [PATCH 07/13] 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 06edf57e243b724f79b4387089fdc3d50dea33ee Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 24 Mar 2025 17:33:03 +0100 Subject: [PATCH 08/13] improving slice to have the ticks --- analysis/slices.py | 98 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 73 insertions(+), 25 deletions(-) diff --git a/analysis/slices.py b/analysis/slices.py index 0456e7b..581927d 100644 --- a/analysis/slices.py +++ b/analysis/slices.py @@ -1,13 +1,18 @@ import numpy as np +import sys +sys.path.append('/home/aubin/Simbelmyne/sbmy_control/') +from cosmo_params import register_arguments_cosmo, parse_arguments_cosmo + fs = 18 fs_titles = fs -4 def plot_imshow_with_reference( data_list, - reference, - titles, + reference=None, + titles=None, vmin=None, - vmax=None, + vmax=None, + L=None, cmap='viridis'): """ Plot the imshow of a list of 2D arrays with two rows: one for the data itself, @@ -23,38 +28,69 @@ def plot_imshow_with_reference( data_list, if titles is None: titles = [None for f in data_list] + + if L is None: + L = [len(data) for data in data_list] + elif isinstance(L, int) or isinstance(L, float): + L = [L for data in data_list] + + sep = 10 if L[0] < 100 else 20 if L[0] < 200 else 100 + ticks = [np.arange(0, l+1, sep)*len(dat)/l for l, dat in zip(L,data_list)] + tick_labels = [np.arange(0, l+1, sep) for l in L] 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)) + fig, axes = plt.subplots(1 if reference is None else 2, n, figsize=(5 * n, 5 if reference is None else 5*2), dpi=max(500, data_list[0].shape[0]//2)) 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) + + if reference is not None: 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') + if reference is not None: + # 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) + axes[0, i].set_xticks(ticks[i]) + axes[0, i].set_yticks(ticks[i]) + axes[0, i].set_xticklabels(tick_labels[i]) + axes[0, i].set_yticklabels(tick_labels[i]) + axes[0, i].set_xlabel('Mpc/h') + 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() + # 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) + axes[1, i].set_xticks(ticks[i]) + axes[1, i].set_yticks(ticks[i]) + axes[1, i].set_xticklabels(tick_labels[i]) + axes[1, i].set_yticklabels(tick_labels[i]) + axes[1, i].set_xlabel('Mpc/h') + 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() + 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[0, i].set_xticks(ticks[i]) + axes[0, i].set_yticks(ticks[i]) + axes[0, i].set_xticklabels(tick_labels[i]) + axes[0, i].set_yticklabels(tick_labels[i]) + fig.colorbar(im, ax=axes[:], orientation='vertical') return fig, axes @@ -76,36 +112,48 @@ if __name__ == "__main__": 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.') + register_arguments_cosmo(parser) + args = parser.parse_args() from pysbmy.field import read_field + from pysbmy.cosmology import d_plus - ref_field = read_field(args.directory+args.reference) + ref_field = read_field(args.directory+args.reference) if args.reference is not None else None 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 + + # args.labels=[f"a={f.time:.2f}" for f in fields] + L = [f.L0 for f in fields] match args.axis: case 0 | 'x': - reference = ref_field.data[index,:,:] + reference = ref_field.data[index,:,:] if ref_field is not None else None fields = [f.data[index,:,:] for f in fields] + # reference = ref_field.data[index,:,:]/d_plus(1e-3,ref_field.time,parse_arguments_cosmo(args)) + # fields = [f.data[index,:,:]/d_plus(1e-3,f.time,parse_arguments_cosmo(args)) for f in fields] + # reference = ref_field.data[index,:,:]/d_plus(1e-3,0.05,parse_arguments_cosmo(args)) + # fields = [f.data[index,:,:]/d_plus(1e-3,time,parse_arguments_cosmo(args)) for f,time in zip(fields,[0.05, 1.0])] case 1 | 'y': - reference = ref_field.data[:,index,:] + reference = ref_field.data[:,index,:] if ref_field is not None else None fields = [f.data[:,index,:] for f in fields] case 2 | 'z': - reference = ref_field.data[:,:,index] + reference = ref_field.data[:,:,index] if ref_field is not None else None 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) + reference = np.log10(2.+reference) if ref_field is not None else None 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, axes = plot_imshow_with_reference(fields,reference,args.labels, vmin=args.vmin, vmax=args.vmax,cmap=args.cmap, L=L) fig.suptitle(args.title) if args.output is not None: From 881bc7b2343f2ff118e0409f61304ba68f63fa56 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 24 Mar 2025 17:33:27 +0100 Subject: [PATCH 09/13] bugfix lightcone output --- parameters_card.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/parameters_card.py b/parameters_card.py index df03911..cb1a14b 100644 --- a/parameters_card.py +++ b/parameters_card.py @@ -188,9 +188,9 @@ def parse_arguments_card(parsed_args): if card_dict["OutputSnapshotsBase"] is None: card_dict["OutputSnapshotsBase"] = main_dict["resultdir"]+"particles_"+main_dict["simname"] if card_dict["OutputFinalSnapshot"] is None: - card_dict["OutputFinalSnapshot"] = main_dict["resultdir"]+ligthcone_prefix+"final_particles_"+main_dict["simname"]+".gadget3" + card_dict["OutputFinalSnapshot"] = main_dict["resultdir"]+ligthcone_prefix+"final_particles_"+main_dict["simname"]+("_lc" if card_dict["GenerateLightcone"] else "")+".gadget3" if card_dict["OutputFinalDensity"] is None: - card_dict["OutputFinalDensity"] = main_dict["resultdir"]+ligthcone_prefix+"final_density_"+main_dict["simname"]+".h5" + card_dict["OutputFinalDensity"] = main_dict["resultdir"]+ligthcone_prefix+"final_density_"+main_dict["simname"]+("_lc" if card_dict["GenerateLightcone"] else "")+".h5" if card_dict["OutputTilesBase"] is None: card_dict["OutputTilesBase"] = main_dict["workdir"]+main_dict["simname"]+"_tile" if card_dict["OutputLPTPotential1"] is None: From 0447279fcc8f428c6f7f3ae1ae5a6995f123b864 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 24 Mar 2025 17:33:42 +0100 Subject: [PATCH 10/13] prevent scola overwriting --- scola.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/scola.py b/scola.py index a6156ad..ff341ab 100644 --- a/scola.py +++ b/scola.py @@ -13,6 +13,11 @@ def main_scola(parsed_args): nboxes_tot = int(parsed_args.N_tiles**3) + if have_all_tiles(parsed_args) and not parsed_args.force: + print_message("All tiles already exist. Use -F to overwrite.", 1, "scola", verbose=parsed_args.verbose) + print_ending_module("scola", verbose=parsed_args.verbose) + return + if parsed_args.execution == "local": from parameters_card import parse_arguments_card from tqdm import tqdm From fa54f87866e7bea85f0811f7ef3a9df83085dda8 Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Mon, 24 Mar 2025 17:34:01 +0100 Subject: [PATCH 11/13] powerspectrum --- analysis/power_spectrum.py | 167 ++++++++++++++++++++++--------------- 1 file changed, 99 insertions(+), 68 deletions(-) diff --git a/analysis/power_spectrum.py b/analysis/power_spectrum.py index 53fdc4f..64dacb5 100644 --- a/analysis/power_spectrum.py +++ b/analysis/power_spectrum.py @@ -5,72 +5,73 @@ kmax = 2e0 Nk = 50 AliasingCorr=False -def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk): +def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None): 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, - ) + if G is None: + 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 + return G, k, Pk -def get_cross_correlations(field_A, field_B, kmin=kmin, kmax=kmax, Nk=Nk): +def get_cross_correlations(field_A, field_B, kmin=kmin, kmax=kmax, Nk=Nk, G=None): 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, - ) + if G is None: + 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 + return G, 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) +def add_power_spectrum_to_plot(ax, field, Pk_ref=None, G=None, plot_args={}, power_args={}): + G, k, Pk = get_power_spectrum(field, G=G, **power_args) if Pk_ref is not None: - ax.plot(k, Pk/Pk_ref-1, **plot_args) + ax.plot(k, Pk/Pk_ref, **plot_args) else: ax.plot(k, Pk, **plot_args) - return ax + return ax, G, k, Pk -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) +def add_cross_correlations_to_plot(ax, field_A, field_B, G=None, plot_args={}, power_args={}): + G, k, Rks = get_cross_correlations(field_A, field_B, G=G, **power_args) ax.plot(k, Rks, **plot_args) - return ax + return ax, G, k, Rks def plot_power_spectra(filenames, @@ -79,6 +80,7 @@ def plot_power_spectra(filenames, linestyles=None, markers=None, Pk_ref=None, + G=None, ylims=[0.9,1.1], yticks = np.linspace(0.9,1.1,11), bound1=0.01, @@ -108,9 +110,10 @@ def plot_power_spectra(filenames, for i, filename in enumerate(filenames): field = read_field(filename) - add_power_spectrum_to_plot(ax=ax, + _, G, k, _ = add_power_spectrum_to_plot(ax=ax, field=field, Pk_ref=Pk_ref, + G=G, plot_args=dict(label=labels[i], color=colors[i], linestyle=linestyles[i], @@ -120,15 +123,15 @@ def plot_power_spectra(filenames, Nk=Nk), ) ax.set_xscale('log') - ax.set_xlim(kmin, kmax) + ax.set_xlim(k.min(),k[-2]) 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_xlabel(r'$k$ [$h/\mathrm{Mpc}$]', labelpad=-10) if Pk_ref is not None: - ax.set_ylabel(r'$P(k)/P_\mathrm{ref}(k)-1$') + ax.set_ylabel(r'$P(k)/P_\mathrm{ref}(k)$') else: ax.set_ylabel('$P(k)$') @@ -145,6 +148,7 @@ def plot_power_spectra(filenames, def plot_cross_correlations(filenames_A, filename_B, + G=None, labels=None, colors=None, linestyles=None, @@ -180,9 +184,10 @@ def plot_cross_correlations(filenames_A, for i, filename_A in enumerate(filenames_A): field_A = read_field(filename_A) - add_cross_correlations_to_plot(ax=ax, + _, G, k, _ = add_cross_correlations_to_plot(ax=ax, field_A=field_A, field_B=field_B, + G=G, plot_args=dict(label=labels[i], color=colors[i], linestyle=linestyles[i], @@ -192,12 +197,12 @@ def plot_cross_correlations(filenames_A, Nk=Nk), ) ax.set_xscale('log') - ax.set_xlim(kmin, kmax) + ax.set_xlim(k.min(), k[-2]) 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_xlabel(r'$k$ [$h/\mathrm{Mpc}$]', labelpad=-10) ax.set_ylabel('$R(k)$') if bound1 is not None: @@ -211,6 +216,28 @@ def plot_cross_correlations(filenames_A, return fig, ax +def get_ylims_and_yticks(ylims): + + if ylims[0] == ylims[1]: + ylims = None + yticks = None + else: + diff_ylims = ylims[1] - ylims[0] + factor = 1 + while diff_ylims<5.: + diff_ylims *= 10 + factor *= 10 + if diff_ylims<12.: + yticks = np.linspace(int(ylims[0]*factor)/factor,int(factor*ylims[1])/factor, int(diff_ylims)+1) + else: + while diff_ylims>12.: + diff_ylims /= 2. + factor /= 2. + yticks = np.linspace(int(ylims[0]*factor)/factor,int(factor*ylims[1])/factor, int(diff_ylims)+1) + + return ylims, yticks + + if __name__ == "__main__": from argparse import ArgumentParser parser = ArgumentParser(description='Plot power spectra of fields') @@ -226,6 +253,8 @@ if __name__ == "__main__": 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.') + 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.') args = parser.parse_args() @@ -235,12 +264,18 @@ if __name__ == "__main__": 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] + G, _, Pk_ref = get_power_spectrum(read_field(args.directory+args.reference), kmin=kmin, kmax=kmax, Nk=Nk) else: Pk_ref = None + G = None + filenames = [args.directory+f for f in args.filenames] + ylims_power, yticks_power = get_ylims_and_yticks(args.ylim_power) + ylims_corr, yticks_corr = get_ylims_and_yticks(args.ylim_corr) + + if args.power_spectrum and args.cross_correlation: import matplotlib.pyplot as plt fig, axes = plt.subplots(2, 1, figsize=(8,8)) @@ -250,14 +285,11 @@ if __name__ == "__main__": 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, + G=G, + ylims=ylims_power, + yticks = yticks_power, + bound1=0.01, + bound2=0.02, kmin=kmin, kmax=kmax, Nk=Nk, @@ -266,18 +298,15 @@ if __name__ == "__main__": plot_cross_correlations(filenames_A=filenames, filename_B=args.directory+args.reference, + G=G, 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, + ylims=ylims_corr, + yticks = yticks_corr, + bound1=0.001, + bound2=0.002, kmin=kmin, kmax=kmax, Nk=Nk, @@ -298,8 +327,9 @@ if __name__ == "__main__": linestyles=args.linestyles, markers=args.markers, Pk_ref=Pk_ref, - ylims=[0.9,1.1], - yticks = np.linspace(0.9,1.1,11), + G=G, + ylims=ylims_power, + yticks = yticks_power, bound1=0.01, bound2=0.02, kmin=kmin, @@ -312,12 +342,13 @@ if __name__ == "__main__": elif args.cross_correlation: fig, ax = plot_cross_correlations(filenames_A=filenames, filename_B=args.reference, + G=G, 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), + ylims=ylims_corr, + yticks = yticks_corr, bound1=0.001, bound2=0.002, kmin=kmin, From 4698184ca2e68f1d91eecb41889a39321ba605bc Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Wed, 26 Mar 2025 14:15:27 +0100 Subject: [PATCH 12/13] default cosmo --- cosmo_params.py | 30 +++++++++++++++++++++++++++++- 1 file changed, 29 insertions(+), 1 deletion(-) diff --git a/cosmo_params.py b/cosmo_params.py index 8ca1636..14e48dc 100644 --- a/cosmo_params.py +++ b/cosmo_params.py @@ -63,4 +63,32 @@ def parse_arguments_cosmo(parsed_args): return cosmo_dict def z2a(z): - return 1.0/(1.0+z) \ No newline at end of file + return 1.0/(1.0+z) + + +cosmo_defaults = { + "h": 0.6732, + "Omega_m": 0.302, + "Omega_b": 0.049, + "Omega_q": 0.6842, + "Omega_k":0.0, + "Omega_r": 0.0, + "n_s": 0.968, + "sigma8": 0.815, + "A_s": 2.148752e-09, + "Tcmb": 2.7255, + "k_p": 0.05, + "N_ur": 2.046, + "m_nu1": 0.06, + "m_nu2": 0.0, + "m_nu3": 0.0, + "w_0": -1.0, + "w_a": 0.0, + "fnl": 100.0, + "gnl": 0.0, + "k_max":10.0, + "tau_reio":0.06, + "WhichSpectrum":"EH", + "w0_fld":-1.0, + "wa_fld":0.0, + } \ No newline at end of file From 64d4ad0ce0f773f7ae31caf697074724fa506d5c Mon Sep 17 00:00:00 2001 From: Mayeul Aubin Date: Wed, 26 Mar 2025 14:15:37 +0100 Subject: [PATCH 13/13] better slice ticks --- analysis/slices.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/analysis/slices.py b/analysis/slices.py index 581927d..8b260cc 100644 --- a/analysis/slices.py +++ b/analysis/slices.py @@ -34,7 +34,7 @@ def plot_imshow_with_reference( data_list, elif isinstance(L, int) or isinstance(L, float): L = [L for data in data_list] - sep = 10 if L[0] < 100 else 20 if L[0] < 200 else 100 + sep = 10 if L[0] < 50 else 20 if L[0] < 200 else 100 ticks = [np.arange(0, l+1, sep)*len(dat)/l for l, dat in zip(L,data_list)] tick_labels = [np.arange(0, l+1, sep) for l in L] @@ -86,10 +86,11 @@ def plot_imshow_with_reference( data_list, 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[0, i].set_xticks(ticks[i]) - axes[0, i].set_yticks(ticks[i]) - axes[0, i].set_xticklabels(tick_labels[i]) - axes[0, i].set_yticklabels(tick_labels[i]) + 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 @@ -123,7 +124,7 @@ if __name__ == "__main__": fields = [read_field(args.directory+f) for f in args.filenames] if args.index is None: - index = ref_field.N0//2 + index = fields[0].N0//2 else: index=args.index