mirror of
https://github.com/Findus23/halo_comparison.git
synced 2024-09-19 16:03:50 +02:00
last week's work
This commit is contained in:
parent
a317408b8b
commit
db4fec38cd
GPG key ID:
54DE4D798D244853
7 changed files with 136 additions and 24 deletions
|
|
@ -12,16 +12,31 @@ from cic import cic_from_radius
|
|||
from cumulative_mass_profiles import cumulative_mass_profile
|
||||
from paths import auriga_dir
|
||||
from readfiles import read_file, read_halo_file
|
||||
from utils import read_swift_config
|
||||
|
||||
|
||||
# softening_length = 0.026041666666666668
|
||||
|
||||
def dir_name_to_parameter(dir_name: str):
|
||||
return map(int, dir_name.lstrip("auriga6_halo").split("_"))
|
||||
|
||||
|
||||
def levelmax_to_softening_length(levelmax: int) -> float:
|
||||
# TODO: replace with real calculation
|
||||
if levelmax == 9:
|
||||
return 0.00651041667
|
||||
if levelmax == 10:
|
||||
return 0.00325520833
|
||||
raise ValueError(f"unknown levelmax {levelmax}")
|
||||
|
||||
softening_length = 0.026041666666666668
|
||||
|
||||
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 / R$_\mathrm{group}$')
|
||||
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}$]")
|
||||
|
||||
|
|
@ -35,10 +50,16 @@ class Result:
|
|||
images = []
|
||||
vmin = np.Inf
|
||||
vmax = -np.Inf
|
||||
dirs = [d for d in auriga_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)
|
||||
|
||||
for i, dir in enumerate(auriga_dir.glob("*")):
|
||||
if "auriga" not in dir.name:
|
||||
continue
|
||||
if not is_by_adrian:
|
||||
levelmin, levelmin_TF, levelmax = dir_name_to_parameter(dir.name)
|
||||
print(levelmin, levelmin_TF, levelmax)
|
||||
# if levelmin_TF != 8:
|
||||
# continue
|
||||
Xc_adrian = 56.50153741810241
|
||||
Yc_adrian = 49.40761085700951
|
||||
Zc_adrian = 49.634393647291695
|
||||
|
|
@ -46,16 +67,20 @@ for i, dir in enumerate(auriga_dir.glob("*")):
|
|||
Yc = 51.34632916228137
|
||||
Zc = 51.68749302578122
|
||||
|
||||
is_by_adrian = "arj" in dir.name
|
||||
if not dir.is_dir():
|
||||
continue
|
||||
input_file = dir / "output_0007.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"]
|
||||
print(levelmax_to_softening_length(levelmax))
|
||||
assert softening_length == levelmax_to_softening_length(levelmax)
|
||||
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
|
||||
|
|
@ -66,6 +91,7 @@ for i, dir in enumerate(auriga_dir.glob("*")):
|
|||
halo_id += 1
|
||||
|
||||
halo = df_halos.loc[halo_id]
|
||||
# halo = halos.loc[1]
|
||||
log_radial_bins, bin_masses, bin_densities, group_radius = cumulative_mass_profile(
|
||||
df, halo, particles_meta, plot=False
|
||||
)
|
||||
|
|
@ -73,8 +99,9 @@ for i, dir in enumerate(auriga_dir.glob("*")):
|
|||
|
||||
ax2.loglog(log_radial_bins[:-1], bin_densities, label=str(dir.name), c=f"C{i}")
|
||||
|
||||
for ax in [ax1,ax2]:
|
||||
ax.axvline(4 * softening_length, color=f"C{i}", linestyle="dotted")
|
||||
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()
|
||||
|
|
@ -115,7 +142,7 @@ for result, ax in zip(images, axes):
|
|||
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")
|
||||
origin="lower")
|
||||
ax.set_title(result.title)
|
||||
|
||||
fig3.tight_layout()
|
||||
|
|
@ -123,8 +150,8 @@ 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("~/tmp/auriga1.pdf").expanduser())
|
||||
fig2.savefig(Path("~/tmp/auriga2.pdf").expanduser())
|
||||
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()
|
||||
|
|
|
|||
|
|
@ -7,6 +7,7 @@ import pandas as pd
|
|||
from matplotlib.axes import Axes
|
||||
from matplotlib.figure import Figure
|
||||
|
||||
from find_center import find_center
|
||||
from readfiles import ParticlesMeta, read_file, read_halo_file
|
||||
|
||||
|
||||
|
|
@ -14,14 +15,15 @@ def V(r):
|
|||
return 4 * np.pi * r ** 3 / 3
|
||||
|
||||
|
||||
def cumulative_mass_profile(particles_in_halos: pd.DataFrame, halo: pd.Series,
|
||||
def cumulative_mass_profile(particles: pd.DataFrame, halo: pd.Series,
|
||||
particles_meta: ParticlesMeta, plot=False, num_bins=30):
|
||||
print(type(particles_in_halos))
|
||||
centre = np.array([halo.X, halo.Y, halo.Z])
|
||||
positions = particles_in_halos[["X", "Y", "Z"]].to_numpy()
|
||||
print(type(particles))
|
||||
center = np.array([halo.X, halo.Y, halo.Z])
|
||||
center = find_center(particles, center)
|
||||
positions = particles[["X", "Y", "Z"]].to_numpy()
|
||||
print(positions)
|
||||
print(positions.shape)
|
||||
distances = np.linalg.norm(positions - centre, axis=1)
|
||||
distances = np.linalg.norm(positions - center, axis=1)
|
||||
group_radius = distances.max()
|
||||
# normalized_distances = distances / group_radius
|
||||
|
||||
|
|
|
|||
31
find_center.py
Normal file
31
find_center.py
Normal file
|
|
@ -0,0 +1,31 @@
|
|||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
||||
|
||||
def find_center(df: pd.DataFrame, center: np.ndarray, initial_radius=1):
|
||||
# plt.figure()
|
||||
all_particles = df[["X", "Y", "Z"]].to_numpy()
|
||||
radius = initial_radius
|
||||
center_history = []
|
||||
while True:
|
||||
center_history.append(center)
|
||||
distances = np.linalg.norm(all_particles - center, axis=1)
|
||||
in_radius_particles = all_particles[distances < radius]
|
||||
num_particles = in_radius_particles.shape[0]
|
||||
print("num_particles", num_particles)
|
||||
if num_particles < 10:
|
||||
break
|
||||
center_of_mass = in_radius_particles.mean(axis=0)
|
||||
new_center = (center_of_mass + center) / 2
|
||||
print("new center", new_center)
|
||||
shift = np.linalg.norm(center - new_center)
|
||||
radius = max(2 * shift, radius * 0.9)
|
||||
print("radius", radius)
|
||||
center = new_center
|
||||
center_history = np.array(center_history)
|
||||
# print(center_history)
|
||||
# plt.scatter(center_history[::, 0], center_history[::, 1], c=range(len(center_history[::, 1])))
|
||||
# plt.colorbar(label="step")
|
||||
# plt.show()
|
||||
return center
|
||||
46
poetry.lock
generated
46
poetry.lock
generated
|
|
@ -320,6 +320,14 @@ all = ["matplotlib", "colorcet", "cmocean", "meshio"]
|
|||
colormaps = ["matplotlib", "colorcet", "cmocean"]
|
||||
io = ["meshio (>=5.2)"]
|
||||
|
||||
[[package]]
|
||||
name = "pyyaml"
|
||||
version = "6.0"
|
||||
description = "YAML parser and emitter for Python"
|
||||
category = "main"
|
||||
optional = false
|
||||
python-versions = ">=3.6"
|
||||
|
||||
[[package]]
|
||||
name = "scipy"
|
||||
version = "1.8.1"
|
||||
|
|
@ -418,10 +426,11 @@ numpy = ["numpy (>=1.9)"]
|
|||
[package.source]
|
||||
type = "url"
|
||||
url = "https://github.com/pyvista/pyvista-wheels/raw/main/vtk-9.1.0.dev0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl"
|
||||
|
||||
[metadata]
|
||||
lock-version = "1.1"
|
||||
python-versions = "^3.9,<3.11"
|
||||
content-hash = "db1d65073afc49ea08202b1a47d7d0747b8db6c2c7235e2b771d6345ba3103fe"
|
||||
content-hash = "33a9a1ef640e87125af6637b5bb5b154f8e1ef79710d2e5e70300fcb3a41d37e"
|
||||
|
||||
[metadata.files]
|
||||
appdirs = [
|
||||
|
|
@ -827,6 +836,41 @@ pyvista = [
|
|||
{file = "pyvista-0.34.1-py3-none-any.whl", hash = "sha256:93d725b84eb037b44cccc45ea65f59bb25f9515238781217614a297f9cdff5fd"},
|
||||
{file = "pyvista-0.34.1.tar.gz", hash = "sha256:3fe322d76c5e74797242026d338f948c4cebeee99efcc817fbf97e39a09ccfb0"},
|
||||
]
|
||||
pyyaml = [
|
||||
{file = "PyYAML-6.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:d4db7c7aef085872ef65a8fd7d6d09a14ae91f691dec3e87ee5ee0539d516f53"},
|
||||
{file = "PyYAML-6.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:9df7ed3b3d2e0ecfe09e14741b857df43adb5a3ddadc919a2d94fbdf78fea53c"},
|
||||
{file = "PyYAML-6.0-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:77f396e6ef4c73fdc33a9157446466f1cff553d979bd00ecb64385760c6babdc"},
|
||||
{file = "PyYAML-6.0-cp310-cp310-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:a80a78046a72361de73f8f395f1f1e49f956c6be882eed58505a15f3e430962b"},
|
||||
{file = "PyYAML-6.0-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:f84fbc98b019fef2ee9a1cb3ce93e3187a6df0b2538a651bfb890254ba9f90b5"},
|
||||
{file = "PyYAML-6.0-cp310-cp310-win32.whl", hash = "sha256:2cd5df3de48857ed0544b34e2d40e9fac445930039f3cfe4bcc592a1f836d513"},
|
||||
{file = "PyYAML-6.0-cp310-cp310-win_amd64.whl", hash = "sha256:daf496c58a8c52083df09b80c860005194014c3698698d1a57cbcfa182142a3a"},
|
||||
{file = "PyYAML-6.0-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:897b80890765f037df3403d22bab41627ca8811ae55e9a722fd0392850ec4d86"},
|
||||
{file = "PyYAML-6.0-cp36-cp36m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:50602afada6d6cbfad699b0c7bb50d5ccffa7e46a3d738092afddc1f9758427f"},
|
||||
{file = "PyYAML-6.0-cp36-cp36m-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:48c346915c114f5fdb3ead70312bd042a953a8ce5c7106d5bfb1a5254e47da92"},
|
||||
{file = "PyYAML-6.0-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:98c4d36e99714e55cfbaaee6dd5badbc9a1ec339ebfc3b1f52e293aee6bb71a4"},
|
||||
{file = "PyYAML-6.0-cp36-cp36m-win32.whl", hash = "sha256:0283c35a6a9fbf047493e3a0ce8d79ef5030852c51e9d911a27badfde0605293"},
|
||||
{file = "PyYAML-6.0-cp36-cp36m-win_amd64.whl", hash = "sha256:07751360502caac1c067a8132d150cf3d61339af5691fe9e87803040dbc5db57"},
|
||||
{file = "PyYAML-6.0-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:819b3830a1543db06c4d4b865e70ded25be52a2e0631ccd2f6a47a2822f2fd7c"},
|
||||
{file = "PyYAML-6.0-cp37-cp37m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:473f9edb243cb1935ab5a084eb238d842fb8f404ed2193a915d1784b5a6b5fc0"},
|
||||
{file = "PyYAML-6.0-cp37-cp37m-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:0ce82d761c532fe4ec3f87fc45688bdd3a4c1dc5e0b4a19814b9009a29baefd4"},
|
||||
{file = "PyYAML-6.0-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:231710d57adfd809ef5d34183b8ed1eeae3f76459c18fb4a0b373ad56bedcdd9"},
|
||||
{file = "PyYAML-6.0-cp37-cp37m-win32.whl", hash = "sha256:c5687b8d43cf58545ade1fe3e055f70eac7a5a1a0bf42824308d868289a95737"},
|
||||
{file = "PyYAML-6.0-cp37-cp37m-win_amd64.whl", hash = "sha256:d15a181d1ecd0d4270dc32edb46f7cb7733c7c508857278d3d378d14d606db2d"},
|
||||
{file = "PyYAML-6.0-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:0b4624f379dab24d3725ffde76559cff63d9ec94e1736b556dacdfebe5ab6d4b"},
|
||||
{file = "PyYAML-6.0-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:213c60cd50106436cc818accf5baa1aba61c0189ff610f64f4a3e8c6726218ba"},
|
||||
{file = "PyYAML-6.0-cp38-cp38-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:9fa600030013c4de8165339db93d182b9431076eb98eb40ee068700c9c813e34"},
|
||||
{file = "PyYAML-6.0-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:277a0ef2981ca40581a47093e9e2d13b3f1fbbeffae064c1d21bfceba2030287"},
|
||||
{file = "PyYAML-6.0-cp38-cp38-win32.whl", hash = "sha256:d4eccecf9adf6fbcc6861a38015c2a64f38b9d94838ac1810a9023a0609e1b78"},
|
||||
{file = "PyYAML-6.0-cp38-cp38-win_amd64.whl", hash = "sha256:1e4747bc279b4f613a09eb64bba2ba602d8a6664c6ce6396a4d0cd413a50ce07"},
|
||||
{file = "PyYAML-6.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:055d937d65826939cb044fc8c9b08889e8c743fdc6a32b33e2390f66013e449b"},
|
||||
{file = "PyYAML-6.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:e61ceaab6f49fb8bdfaa0f92c4b57bcfbea54c09277b1b4f7ac376bfb7a7c174"},
|
||||
{file = "PyYAML-6.0-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:d67d839ede4ed1b28a4e8909735fc992a923cdb84e618544973d7dfc71540803"},
|
||||
{file = "PyYAML-6.0-cp39-cp39-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:cba8c411ef271aa037d7357a2bc8f9ee8b58b9965831d9e51baf703280dc73d3"},
|
||||
{file = "PyYAML-6.0-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:40527857252b61eacd1d9af500c3337ba8deb8fc298940291486c465c8b46ec0"},
|
||||
{file = "PyYAML-6.0-cp39-cp39-win32.whl", hash = "sha256:b5b9eccad747aabaaffbc6064800670f0c297e52c12754eb1d976c57e4f74dcb"},
|
||||
{file = "PyYAML-6.0-cp39-cp39-win_amd64.whl", hash = "sha256:b3d267842bf12586ba6c734f89d1f5b871df0273157918b0ccefa29deb05c21c"},
|
||||
{file = "PyYAML-6.0.tar.gz", hash = "sha256:68fb519c14306fec9720a2a5b45bc9f0c8d1b9c72adf45c37baedfcd949c35a2"},
|
||||
]
|
||||
scipy = [
|
||||
{file = "scipy-1.8.1-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:65b77f20202599c51eb2771d11a6b899b97989159b7975e9b5259594f1d35ef4"},
|
||||
{file = "scipy-1.8.1-cp310-cp310-macosx_12_0_arm64.whl", hash = "sha256:e013aed00ed776d790be4cb32826adb72799c61e318676172495383ba4570aa4"},
|
||||
|
|
|
|||
|
|
@ -17,6 +17,7 @@ pynbody = "^1.1.0"
|
|||
numba = "^0.55.1"
|
||||
tables = "^3.7.0"
|
||||
pythran = "^0.11.0"
|
||||
PyYAML = "^6.0"
|
||||
|
||||
[tool.poetry.dev-dependencies]
|
||||
|
||||
|
|
|
|||
|
|
@ -83,13 +83,13 @@ def cached_particles_in_halo(file: Path, *args, **kwargs) -> HaloParticleMapping
|
|||
return halo_particle_ids
|
||||
|
||||
|
||||
def read_velo_halos(directory: Path):
|
||||
def read_velo_halos(directory: Path,veloname="vroutput"):
|
||||
"""
|
||||
Returns a dataframe containing all scalar properties of the halos
|
||||
(https://velociraptor-stf.readthedocs.io/en/latest/output.html),
|
||||
"""
|
||||
group_catalog = h5py.File(directory / "vroutput.catalog_groups")
|
||||
group_properties = h5py.File(directory / "vroutput.properties")
|
||||
group_catalog = h5py.File(directory / f"{veloname}.catalog_groups")
|
||||
group_properties = h5py.File(directory / f"{veloname}.properties")
|
||||
scalar_properties = {}
|
||||
for k, v in group_properties.items():
|
||||
if not isinstance(v, Dataset):
|
||||
|
|
|
|||
7
utils.py
7
utils.py
|
|
@ -1,6 +1,8 @@
|
|||
from pathlib import Path
|
||||
from typing import Tuple
|
||||
|
||||
import pandas as pd
|
||||
import yaml
|
||||
from matplotlib import pyplot as plt
|
||||
from matplotlib.axes import Axes
|
||||
from matplotlib.figure import Figure
|
||||
|
|
@ -14,6 +16,7 @@ def memory_usage(df: pd.DataFrame):
|
|||
bytes_used = df.memory_usage(index=True).sum()
|
||||
return bytes_used / 1024 / 1024
|
||||
|
||||
|
||||
def create_figure() -> Tuple[Figure, Axes]:
|
||||
"""
|
||||
helper function for matplotlib OOP interface with proper typing
|
||||
|
|
@ -22,3 +25,7 @@ def create_figure() -> Tuple[Figure, Axes]:
|
|||
ax: Axes = fig.gca()
|
||||
return fig, ax
|
||||
|
||||
|
||||
def read_swift_config(dir: Path):
|
||||
with (dir / "used_parameters.yml").open() as f:
|
||||
return yaml.safe_load(f)
|
||||
|
|
|
|||
Loading…
Reference in a new issue