halo_comparison/auriga_comparison.py

import pickle
from dataclasses import dataclass
from enum import Enum
from pathlib import Path
from typing import List

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.axes import Axes
from matplotlib.colors import LogNorm
from matplotlib.figure import Figure

from cic import cic_from_radius
from halo_mass_profile import halo_mass_profile
from paths import auriga_dir, richings_dir
from readfiles import read_file, read_halo_file
from utils import read_swift_config


class Mode(Enum):
    richings = 1
    auriga6 = 2


mode = Mode.richings


def dir_name_to_parameter(dir_name: str):
    return map(int, dir_name.lstrip("auriga6_halo").lstrip("richings21_").split("_"))


def levelmax_to_softening_length(levelmax: int) -> float:
    box_size = 100
    return box_size / 30 / 2 ** levelmax


fig1: Figure = plt.figure(figsize=(9, 6))
ax1: Axes = fig1.gca()
fig2: Figure = plt.figure(figsize=(9, 6))
ax2: Axes = fig2.gca()

for ax in [ax1, ax2]:
    ax.set_xlabel(r'R [Mpc]')
ax1.set_ylabel(r'M [$10^{10} \mathrm{M}_\odot$]')
ax2.set_ylabel("density [$\\frac{10^{10} \\mathrm{M}_\\odot}{Mpc^3}$]")

part_numbers = []

reference_file = Path("auriga_reference.pickle")


@dataclass
class Result:
    title: str
    rho: np.ndarray


images = []
vmin = np.Inf
vmax = -np.Inf
root_dir = auriga_dir if mode == Mode.auriga6 else richings_dir
dirs = [d for d in root_dir.glob("*") if d.is_dir() and "bak" not in d.name]
for i, dir in enumerate(sorted(dirs)):
    is_by_adrian = "arj" in dir.name
    print(dir.name)

    if not is_by_adrian:
        levelmin, levelmin_TF, levelmax = dir_name_to_parameter(dir.name)
        print(levelmin, levelmin_TF, levelmax)
        # if levelmax != 12:
        #     continue
    Xc_adrian = 56.50153741810241
    Yc_adrian = 49.40761085700951
    Zc_adrian = 49.634393647291695
    Xc = 58.25576087992683
    Yc = 51.34632916228137
    Zc = 51.68749302578122

    input_file = dir / "output_0007.hdf5"
    if mode == Mode.richings:
        input_file = dir / "output_0004.hdf5"
    if is_by_adrian:
        input_file = dir / "output_0000.hdf5"
        softening_length = None
    else:
        swift_conf = read_swift_config(dir)
        softening_length = swift_conf["Gravity"]["comoving_DM_softening"]
        assert softening_length == swift_conf["Gravity"]["max_physical_DM_softening"]
        ideal_softening_length = levelmax_to_softening_length(levelmax)
        # if not np.isclose(softening_length, levelmax_to_softening_length(levelmax)):
        #     raise ValueError(f"softening length for levelmax {levelmax} should be {ideal_softening_length} "
        #                      f"but is {softening_length}")
    print(input_file)
    df, particles_meta = read_file(input_file)
    df_halos = read_halo_file(input_file.with_name("fof_" + input_file.name))
    # halos = read_velo_halos(dir, veloname="velo_out")
    # particles_in_halo = df.loc[df["FOFGroupIDs"] == 3]

    halo_id = 1
    while True:
        particles_in_halo = df.loc[df["FOFGroupIDs"] == halo_id]
        if len(particles_in_halo) > 1:
            break
        halo_id += 1

    halo = df_halos.loc[halo_id]
    part_numbers.append(len(df) * particles_meta.particle_mass)
    # halo = halos.loc[1]
    log_radial_bins, bin_masses, bin_densities, center = halo_mass_profile(
        df, halo, particles_meta, plot=False, num_bins=100,
        vmin=0.002, vmax=6.5
    )
    if is_by_adrian:
        with reference_file.open("wb") as f:
            pickle.dump([log_radial_bins, bin_masses, bin_densities], f)
    ax1.loglog(log_radial_bins[:-1], bin_masses, label=str(dir.name), c=f"C{i}")

    ax2.loglog(log_radial_bins[:-1], bin_densities, label=str(dir.name), c=f"C{i}")

    if reference_file.exists() and not is_by_adrian and mode == Mode.auriga6:
        with reference_file.open("rb") as f:
            data: List[np.ndarray] = pickle.load(f)
            ref_log_radial_bins, ref_bin_masses, ref_bin_densities = data
        mass_deviation: np.ndarray = np.abs(bin_masses - ref_bin_masses)
        density_deviation: np.ndarray = np.abs(bin_densities - ref_bin_densities)
        ax1.loglog(log_radial_bins[:-1], mass_deviation, c=f"C{i}",
                   linestyle="dotted")

        ax2.loglog(log_radial_bins[:-1], density_deviation, c=f"C{i}",
                   linestyle="dotted")
        accuracy = mass_deviation / ref_bin_masses
        print(accuracy)
        print("mean accuracy", accuracy.mean())

    if softening_length:
        for ax in [ax1, ax2]:
            ax.axvline(4 * softening_length, color=f"C{i}", linestyle="dotted")

    X, Y, Z = df.X.to_numpy(), df.Y.to_numpy(), df.Z.to_numpy()
    print()
    print(Yc - Yc_adrian)
    # shift: (-6, 0, -12)
    # if not is_by_adrian:
    #     xshift = Xc - Xc_adrian
    #     yshift = Yc - Yc_adrian
    #     zshift = Zc - Zc_adrian
    #     print("shift", xshift, yshift, zshift)

    X -= center[0]
    Y -= center[1]
    Z -= center[2]

    rho, extent = cic_from_radius(X, Z, 500, 0, 0, 5, periodic=False)

    vmin = min(vmin, rho.min())
    vmax = max(vmax, rho.max())

    images.append(Result(
        rho=rho,
        title=str(dir.name)
    ))

    # plot_cic(
    #     rho, extent,
    #     title=str(dir.name)
    # )
ax1.legend()
ax2.legend()

# fig3: Figure = plt.figure(figsize=(9, 9))
# axes: List[Axes] = fig3.subplots(3, 3, sharex=True, sharey=True).flatten()
fig3: Figure = plt.figure(figsize=(6, 3))
axes: List[Axes] = fig3.subplots(1, 2, sharex=True, sharey=True).flatten()

for result, ax in zip(images, axes):
    data = 1.1 + result.rho
    vmin_scaled = 1.1 + vmin
    vmax_scaled = 1.1 + vmax
    img = ax.imshow(data.T, norm=LogNorm(vmin=vmin_scaled, vmax=vmax_scaled), extent=extent,
                    origin="lower")
    ax.set_title(result.title)

fig3.tight_layout()
fig3.subplots_adjust(right=0.825)
cbar_ax = fig3.add_axes([0.85, 0.15, 0.05, 0.7])
fig3.colorbar(img, cax=cbar_ax)

fig1.savefig(Path(f"~/tmp/auriga1_{8}.pdf").expanduser())
fig2.savefig(Path(f"~/tmp/auriga2_{8}.pdf").expanduser())
fig3.savefig(Path("~/tmp/auriga3.pdf").expanduser())

plt.show()
print(part_numbers)