diff --git a/analysis/__init__.py b/analysis/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/analysis/power_spectrum.py b/analysis/power_spectrum.py deleted file mode 100644 index 53fdc4f..0000000 --- a/analysis/power_spectrum.py +++ /dev/null @@ -1,335 +0,0 @@ -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 deleted file mode 100644 index 0456e7b..0000000 --- a/analysis/slices.py +++ /dev/null @@ -1,114 +0,0 @@ -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 diff --git a/main.py b/main.py index 54be5a7..ff18d4f 100644 --- a/main.py +++ b/main.py @@ -190,9 +190,7 @@ 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) diff --git a/parameters_card.py b/parameters_card.py index fe1768c..fa7a6c7 100644 --- a/parameters_card.py +++ b/parameters_card.py @@ -48,14 +48,12 @@ 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=False, help="Write gravitational potential.") + 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.") - 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): @@ -134,8 +132,6 @@ 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"], @@ -211,7 +207,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"]+main_dict["simname"]+"_tile" + card_dict["OutputTilesBase"] = main_dict["workdir"]+"sCOLA_tile" if card_dict["OutputLPTPotential1"] is None: card_dict["OutputLPTPotential1"] = main_dict["workdir"]+"initial_conditions_DM_phi.h5" if card_dict["OutputLPTPotential2"] is None: @@ -231,8 +227,6 @@ 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 @@ -358,8 +352,6 @@ 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, @@ -429,8 +421,6 @@ def create_parameter_card_dict( WriteDensity=int(WriteDensity), OutputDensityBase=OutputDensityBase, MeshDensity=MeshDensity, - LoadPhiBCs=int(LoadPhiBCs), - InputPhiBCsBase=InputPhiBCsBase, h=h, Omega_m=Omega_m, diff --git a/scola.py b/scola.py index 87e2f43..3d1af5c 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_"+main_dict["simname"]+".sh" + slurm_script = slurm_dict["scripts"]+"scola.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_"+main_dict["simname"]+"_"+str(missing_tiles[0])+"_"+str(missing_tiles[-1])+".sh" + slurm_script = slurm_dict["scripts"]+"scola_"+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)