look into spectra and HMF

halo_mass_functions.py

@@ -2,9 +2,8 @@ from math import log
 from pathlib import Path
 
 import numpy as np
-
-# from colossus.cosmology import cosmology
-# from colossus.lss import mass_function
+from colossus.cosmology import cosmology
+from colossus.lss import mass_function
 from matplotlib import pyplot as plt
 from matplotlib.axes import Axes
 from matplotlib.figure import Figure
@@ -33,14 +32,30 @@ def monofonic_tests():
     else:
         resolutions.append(512)
 
-    for i, waveform in enumerate(["DB2", "shannon"]):
+    plank_cosmo = cosmology.cosmologies["planck18"]
+    our_cosmo = {
+        # "sigma8": 0.807,
+        "H0": 0.67742 * 100,
+        "Om0": 0.3099,
+        "Ob0": 0.048891054,
+        "ns": 0.96822,
+    }
+
+    cosmology.addCosmology("ourCosmo", params={**plank_cosmo, **our_cosmo})
+    cosmology.setCosmology("ourCosmo")
+    print(cosmology.getCurrent())
+
+    for i, waveform in enumerate(["DB2", "shannon", "shannon_rehalo"]):
         for j, resolution in enumerate(resolutions):
+            if waveform == "shannon_rehalo" and resolution != 128:
+                continue
             print(waveform, resolution)
             dir = base_dir / f"{waveform}_{resolution}_100"
             halos = read_velo_halos(dir)
 
             # halos.to_csv("weird_halos.csv")
-            halo_masses: np.ndarray = halos["Mvir"].to_numpy()
+            halo_masses: np.ndarray = halos["Mvir"].to_numpy() * 1e10 * 0.67742
+            halo_masses = halo_masses[halo_masses > 0]
 
             (
                 Ns,
@@ -58,7 +73,7 @@ def monofonic_tests():
             name = f"{waveform} {resolution}"
             ax.step(
                 left_edges,
-                number_densities,
+                number_densities / (0.67742 ** 3),
                 where="post",
                 color=f"C{i}",
                 linestyle=linestyles[j],
@@ -76,6 +91,13 @@ def monofonic_tests():
 
             # break
         # break
+
+    mfunc = mass_function.massFunction(
+        left_edges, 1, mdef="vir", model="tinker08", q_out="dndlnM"
+    )
+
+    ax.plot(left_edges, mfunc, label="tinker08", color="C4")
+
     plt.legend()
     fig.savefig(Path(f"~/tmp/halo_mass_function.pdf").expanduser())
     plt.show()

hmf_test.py

@@ -0,0 +1,28 @@
+import matplotlib.pyplot as plt
+from velociraptor import load, tools
+
+from paths import base_dir
+from read_vr_files import read_velo_halos
+from unyt import Mpc
+import unyt
+
+for i, waveform in enumerate(["DB2", "DB4", "DB8", "shannon"]):
+    for j, resolution in enumerate([128, 256, 512]):
+        print(waveform, resolution)
+        dir = base_dir / f"{waveform}_{resolution}_100"
+        data = load(str(dir)+"/vroutput.properties.0")
+        print(data)
+        masses_200crit = data.masses.mass_200crit
+        masses_200crit.convert_to_units("msun")
+        lowest_halo_mass = 1e9 * unyt.msun
+        highest_halo_mass = 8e15 * unyt.msun
+        bin_centers, mass_function, error = tools.create_mass_function(
+            masses_200crit, lowest_halo_mass, highest_halo_mass, (100*Mpc)**3
+        )
+        plt.loglog(bin_centers,mass_function)
+        # halos = read_velo_halos(dir)
+        #
+        # print(halos["Rvir"].sort_values(ascending=False)[:4].to_numpy())
+
+
+plt.show()

spectra_test.py

@@ -0,0 +1,83 @@
+import numpy as np
+import pandas as pd
+import scipy.special as sf
+from matplotlib import pyplot as plt
+from matplotlib.axes import Axes
+from matplotlib.figure import Figure
+
+
+def read_monofonic_file(file: str):
+    return pd.read_csv(
+        file,
+        sep=" ",
+        skipinitialspace=True,
+        header=None,
+        names=[
+            "k", "P_dtot", "P_dcdm", "P_dbar", "P_tcdm", "P_tbar", "..", "...", "....", "....."
+        ],
+        skiprows=1,
+    )
+
+
+def read_spectra_file(file: str):
+    return pd.read_csv(
+        # f"/home/lukas/cosmos_data/monofonic_tests/spectra/DB8_100/DB8_100_ics_256_256_cross_spectrum.txt",
+        file,
+        sep=" ",
+        skipinitialspace=True,
+        header=None,
+        names=[
+            "k",
+            "Pcross",
+            "P1",
+            "err. P1",
+            "P2",
+            "err. P2",
+            "P2-1",
+            "err. P2-1",
+            "modes in bin",
+        ],
+        skiprows=1,
+    )
+
+
+monofonic_spectra = read_monofonic_file(
+    f"/home/lukas/cosmoca/cosmoca-reproducibility/monofonic_tests/input_powerspec.txt")
+print(monofonic_spectra)
+
+
+def Dplus(lambda0, a):
+    return a * np.sqrt(1.0 + lambda0 * a ** 3) * \
+        sf.hyp2f1(3.0 / 2.0, 5.0 / 6.0, 11.0 / 6.0, -lambda0 * a ** 3)
+
+
+OmegaLambda = 0.6901
+OmegaM = 0.3099
+lambda_0 = OmegaLambda / OmegaM
+a_ini = 0.02
+factor = Dplus(lambda_0, a_ini) ** 2 / Dplus(lambda_0, 1) ** 2
+
+sim_spectra = read_spectra_file(
+    "/home/lukas/cosmoca/cosmoca-reproducibility/monofonic_tests/out.txt_cross_spectrum.txt")
+sim_spectra_ics_new = read_spectra_file(
+    "/home/lukas/cosmoca/cosmoca-reproducibility/monofonic_tests/DB8_spectra.txt_cross_spectrum.txt")
+sim_spectra_ics_old = read_spectra_file(
+    "/home/lukas/cosmos_data/monofonic_tests/spectra/DB8_100/DB8_100_ics_256_256_cross_spectrum.txt", )
+sim_spectra = read_spectra_file(
+    "/home/lukas/cosmos_data/monofonic_tests/spectra/DB8_100/DB8_100_a4_256_256_cross_spectrum.txt", )
+sim_spectra_new = read_spectra_file(
+    "/home/lukas/cosmoca/cosmoca-reproducibility/monofonic_tests/DB8_spectra_end_cross_spectrum.txt")
+
+fig: Figure = plt.figure()
+ax: Axes = fig.gca(
+
+)
+ax.loglog(sim_spectra.k, sim_spectra.P1*factor, label="sim")
+ax.loglog(sim_spectra_new.k, sim_spectra_new.P1*factor, label="sim new")
+ax.loglog(sim_spectra_ics_old.k, sim_spectra_ics_old.P1, label="ics old")
+ax.loglog(sim_spectra_ics_new.k, sim_spectra_ics_new.P1, label="ics new")
+ax.loglog(monofonic_spectra.k, monofonic_spectra.P_dtot, linestyle="dotted", label="ref")
+# ax.loglog(monofonic_spectra.k, monofonic_spectra.P_dcdm, label="P_dcdm")
+# ax.loglog(monofonic_spectra.k, monofonic_spectra.P_dbar, label="P_dbar")
+plt.legend()
+plt.show()