analysis/power_spectrum.py

@@ -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')
-    
-

analysis/slices.py

@@ -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')

main.py

@@ -190,8 +190,6 @@ 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)
 

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,

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)