convert snapshot to density script

analysis/power_spectrum.py

@@ -15,10 +15,10 @@ def get_power_spectrum(field, kmin=kmin, kmax=kmax, Nk=Nk, G=None):
             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(
+            int(field.N0),
+            int(field.N1),
+            int(field.N2),
+            k_modes=np.concat([PowerSpectrum(field.L0,field.L1,field.L2,int(field.N0),int(field.N1),int(field.N2),).FourierGrid.k_modes[:10],np.logspace(
                 np.log10(kmin),
                 np.log10(kmax),
                 Nk,
@@ -43,10 +43,10 @@ def get_cross_correlations(field_A, field_B, kmin=kmin, kmax=kmax, Nk=Nk, G=None
             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(
+            int(field_A.N0),
+            int(field_A.N1),
+            int(field_A.N2),
+            k_modes=np.concat([PowerSpectrum(field_A.L0,field_A.L1,field_A.L2,int(field_A.N0),int(field_A.N1),int(field_A.N2),).FourierGrid.k_modes[:10],np.logspace(
                 np.log10(kmin),
                 np.log10(kmax),
                 Nk,

analysis/slices.py

@@ -25,6 +25,7 @@ def plot_imshow_with_reference( data_list,
     - cmap: colormap to be used for plotting
     """
     import matplotlib.pyplot as plt
+    from matplotlib import ticker
 
     if titles is None:
         titles = [None for f in data_list]
@@ -65,6 +66,8 @@ def plot_imshow_with_reference( data_list,
             axes[0, i].set_xticklabels(tick_labels[i])
             axes[0, i].set_yticklabels(tick_labels[i])
             axes[0, i].set_xlabel('Mpc/h')
+            axes[0, i].xaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
+            axes[0, i].yaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
         fig.colorbar(im, ax=axes[0, :], orientation='vertical')
 
         # Plot the data compared to the reference
@@ -76,11 +79,13 @@ def plot_imshow_with_reference( data_list,
             axes[1, i].set_xticklabels(tick_labels[i])
             axes[1, i].set_yticklabels(tick_labels[i])
             axes[1, i].set_xlabel('Mpc/h')
+            axes[1, i].xaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
+            axes[1, i].yaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
         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')
+            axes[1, i].text(0.5, 0.9, f"RMS: {score(data, reference):.2e}", fontsize=10, transform=axes[1, i].transAxes, color='white')
         # plt.tight_layout()
     else:
 
@@ -93,6 +98,8 @@ def plot_imshow_with_reference( data_list,
             axes.set_xticklabels(tick_labels[0])
             axes.set_yticklabels(tick_labels[0])
             axes.set_xlabel('Mpc/h')
+            axes.xaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
+            axes.yaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
             fig.colorbar(im, ax=axes, orientation='vertical')
 
         else:
@@ -104,6 +111,8 @@ def plot_imshow_with_reference( data_list,
                 axes[i].set_xticklabels(tick_labels[i])
                 axes[i].set_yticklabels(tick_labels[i])
                 axes[i].set_xlabel('Mpc/h')
+                axes[i].xaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
+                axes[i].yaxis.set_major_formatter(ticker.FormatStrFormatter('%d'))
             fig.colorbar(im, ax=axes[:], orientation='vertical')
 
     return fig, axes
@@ -171,6 +180,6 @@ if __name__ == "__main__":
     fig.suptitle(args.title)
     
     if args.output is not None:
-        fig.savefig(args.output)
+        fig.savefig(args.output,bbox_inches='tight')
     else:
-        fig.savefig(args.directory+'slices.png')
+        fig.savefig(args.directory+'slices.jpg',bbox_inches='tight')

convert_snapshot_to_density.py

@@ -0,0 +1,54 @@
+from pysbmy.density import get_density_pm_snapshot
+from pysbmy.snapshot import read_snapshot
+import argparse
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description="Convert snapshot to density.")
+    parser.add_argument(
+        "-S",
+        "--snapshot",
+        type=str,
+        required=True,
+        help="Path to the snapshot file.",
+    )
+    parser.add_argument(
+        "-o",
+        "--output",
+        type=str,
+        required=True,
+        help="Path to the output density file.",
+    )
+    parser.add_argument(
+        "-N",
+        "--N",
+        type=int,
+        default=None,
+        help="Size of the density field grid (N x N x N).",
+    )
+    parser.add_argument(
+        "-c",
+        "--corner",
+        type=float,
+        nargs=3,
+        default=[0.0, 0.0, 0.0],
+        help="Corner of the box (x, y, z).",
+    )
+
+    args = parser.parse_args()
+
+    # Read the snapshot
+    print("Reading snapshot...")
+    snap = read_snapshot(args.snapshot)
+    
+    if args.N is None:
+        N = snap.Np0
+    else:
+        N = args.N
+
+    print("Calculating density...")
+    F=get_density_pm_snapshot(snap, N,N,N, args.corner[0],args.corner[1],args.corner[2])
+
+    print("Writing density...")
+    F.write(args.output)
+    print("Density written to", args.output)
+    print("Done.")

low_level.py

@@ -225,6 +225,8 @@ def get_progress_from_logfile(filename):
                 pass
         elif "Fatal" in line or "Error" in line:
             return -1, -1
+        elif "Everything done successfully, exiting." in line:
+            current_operation = total_operations
     return current_operation, total_operations
 
 

parameters_monofonic.py

@@ -71,7 +71,7 @@ def get_config_from_dict(monofonic_dict):
     config["setup"] = {
         "GridRes": monofonic_dict["gridres"],
         "BoxLength": monofonic_dict["boxlength"],
-        "zstart": 999.0,
+        "zstart": 99.0,
         "LPTorder": 2,
         "DoBaryons": False,
         "DoBaryonVrel": False,