From ffdc64db74ad19d154821f894266be3f26a13654 Mon Sep 17 00:00:00 2001 From: Lukas Winkler Date: Tue, 13 Dec 2022 11:42:07 +0100 Subject: [PATCH] look into spectra and HMF --- halo_mass_functions.py | 34 ++++++++++++++--- hmf_test.py | 28 ++++++++++++++ spectra_test.py | 83 ++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 139 insertions(+), 6 deletions(-) create mode 100644 hmf_test.py create mode 100644 spectra_test.py diff --git a/halo_mass_functions.py b/halo_mass_functions.py index 10dcc1e..a4e1335 100644 --- a/halo_mass_functions.py +++ b/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() diff --git a/hmf_test.py b/hmf_test.py new file mode 100644 index 0000000..c69dbb9 --- /dev/null +++ b/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() diff --git a/spectra_test.py b/spectra_test.py new file mode 100644 index 0000000..a4729c5 --- /dev/null +++ b/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()