mirror of
https://github.com/Findus23/halo_comparison.git
synced 2024-09-19 16:03:50 +02:00
minor changes
This commit is contained in:
parent
e3a0e4dd3d
commit
e1c77488e1
GPG key ID:
54DE4D798D244853
7 changed files with 126 additions and 40 deletions
49
cic.py
49
cic.py
|
|
@ -1,3 +1,4 @@
|
|||
import sys
|
||||
from pathlib import Path
|
||||
from typing import Tuple
|
||||
|
||||
|
|
@ -8,18 +9,21 @@ from matplotlib.colors import LogNorm
|
|||
from matplotlib.figure import Figure
|
||||
from numba import njit
|
||||
|
||||
from paths import base_dir
|
||||
from readfiles import read_file
|
||||
|
||||
Extent = Tuple[float, float, float, float]
|
||||
|
||||
|
||||
@njit()
|
||||
def cic_deposit(X, Y, ngrid) -> np.ndarray:
|
||||
@njit(parallel=True)
|
||||
def cic_deposit(X, Y, ngrid, periodic=True) -> np.ndarray:
|
||||
rho = np.zeros((ngrid, ngrid))
|
||||
for x, y in zip(X, Y):
|
||||
x = np.fmod(1.0 + x, 1.0)
|
||||
y = np.fmod(1.0 + y, 1.0)
|
||||
if periodic:
|
||||
x = np.fmod(1.0 + x, 1.0)
|
||||
y = np.fmod(1.0 + y, 1.0)
|
||||
else:
|
||||
if not (0 < x < 1) or not (0 < y < 1):
|
||||
continue
|
||||
il = int(np.floor(x * ngrid))
|
||||
ir = (il + 1) % ngrid
|
||||
jl = int(np.floor(y * ngrid))
|
||||
|
|
@ -38,44 +42,49 @@ def cic_range(
|
|||
X: np.ndarray, Y: np.ndarray,
|
||||
ngrid: int,
|
||||
xmin: float, xmax: float,
|
||||
ymin: float, ymax: float
|
||||
ymin: float, ymax: float, *args, **kwargs
|
||||
) -> Tuple[np.ndarray, Extent]:
|
||||
xrange = xmax - xmin
|
||||
yrange = ymax - ymin
|
||||
Xs = (X - xmin) / xrange
|
||||
Ys = (Y - ymin) / yrange
|
||||
print(Xs.min(), Xs.max())
|
||||
print(Ys.min(), Ys.max())
|
||||
print((60 - ymin) / yrange)
|
||||
extent = (xmin, xmax, ymin, ymax)
|
||||
return cic_deposit(Xs, Ys, ngrid), extent
|
||||
rho = cic_deposit(Xs, Ys, ngrid, *args, **kwargs)
|
||||
print(rho)
|
||||
exit()
|
||||
return rho, extent
|
||||
|
||||
|
||||
def cic_from_radius(
|
||||
X: np.ndarray, Y: np.ndarray,
|
||||
ngrid: int,
|
||||
x_center: float, y_center: float,
|
||||
radius: float
|
||||
|
||||
radius: float, *args, **kwargs
|
||||
) -> Tuple[np.ndarray, Extent]:
|
||||
return cic_range(X, Y, ngrid, x_center - radius, x_center + radius, y_center - radius, y_center + radius)
|
||||
return cic_range(
|
||||
X, Y, ngrid,
|
||||
x_center - radius, x_center + radius,
|
||||
y_center - radius, y_center + radius,
|
||||
*args, **kwargs
|
||||
)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
reference_dir = Path(base_dir / f"shannon_512_100")
|
||||
reference_dir = Path(sys.argv[1])
|
||||
df_ref, _ = read_file(reference_dir)
|
||||
|
||||
# quick hack filter
|
||||
df_ref = df_ref[df_ref["X"] < 40]
|
||||
df_ref = df_ref[df_ref["X"] > 30]
|
||||
df_ref = df_ref[df_ref["Y"] < 40]
|
||||
df_ref = df_ref[df_ref["Y"] > 30]
|
||||
|
||||
fig: Figure = plt.figure()
|
||||
ax: Axes = fig.gca()
|
||||
|
||||
print("start cic")
|
||||
rho, extent = cic_range(df_ref.X.to_numpy(), df_ref.Y.to_numpy(), 1000, 30, 40, 30, 40)
|
||||
rho, extent = cic_from_radius(df_ref.X.to_numpy(), df_ref.Y.to_numpy(), 2000, 48.85, 56.985, 2, periodic=False)
|
||||
# rho, extent = cic_range(df_ref.X.to_numpy(), df_ref.Y.to_numpy(), 1000, 30, 70, 30, 70, periodic=False)
|
||||
print("finished cic")
|
||||
data = 1.001 + rho
|
||||
i = ax.imshow(data, norm=LogNorm(), extent=extent)
|
||||
data = 1.1 + rho
|
||||
i = ax.imshow(data.T, norm=LogNorm(), extent=extent, origin="lower")
|
||||
fig.colorbar(i)
|
||||
plt.show()
|
||||
|
||||
|
|
|
|||
|
|
@ -81,10 +81,13 @@ def compare_halo_resolutions(
|
|||
print(f"Memory ref: {memory_usage(df_ref):.2f} MB")
|
||||
print(f"Memory comp: {memory_usage(df_comp):.2f} MB")
|
||||
|
||||
comp_halo_masses=dict(df_comp_halo["Mvir"])
|
||||
|
||||
for index, original_halo in df_ref_halo.iterrows():
|
||||
print(f"{index} of {len(df_ref_halo)} original halos")
|
||||
halo_particle_ids = ref_halo_lookup[int(index)]
|
||||
ref_halo: pd.Series = df_ref_halo.loc[index]
|
||||
ref_halo_mass=ref_halo["Mvir"]
|
||||
if ref_halo["cNFW"] < 0:
|
||||
print("NEGATIVE")
|
||||
print(ref_halo["cNFW"])
|
||||
|
|
@ -196,12 +199,19 @@ def compare_halo_resolutions(
|
|||
# plt.show()
|
||||
best_halo = None
|
||||
best_halo_match = 0
|
||||
|
||||
num_skipped_for_mass = 0
|
||||
for i, halo_id in enumerate(nearby_halos):
|
||||
# print("----------", halo, "----------")
|
||||
# halo_data = df_comp_halo.loc[halo]
|
||||
# particles_in_comp_halo: DataFrame = df_comp.loc[df_comp["FOFGroupIDs"] == halo]
|
||||
particle_ids_in_comp_halo = comp_halo_lookup[halo_id]
|
||||
mass_factor_limit = 5
|
||||
|
||||
if not (1 / mass_factor_limit < (comp_halo_masses[halo_id] / ref_halo_mass) < mass_factor_limit):
|
||||
# print("mass not similar, skipping")
|
||||
num_skipped_for_mass += 1
|
||||
continue
|
||||
|
||||
halo_size = len(particle_ids_in_comp_halo)
|
||||
# df = particles_in_comp_halo.join(halo_particles, how="inner", rsuffix="ref")
|
||||
shared_particles = particle_ids_in_comp_halo.intersection(halo_particle_ids)
|
||||
|
|
@ -217,7 +227,6 @@ def compare_halo_resolutions(
|
|||
|
||||
ax.scatter(apply_offset_to_list(df["X"], offset_x), apply_offset_to_list(df["Y"], offset_y), s=1,
|
||||
alpha=.3, c=color)
|
||||
comp_halo = df_comp_halo.loc[halo_id]
|
||||
circle = Circle((apply_offset(comp_halo.X, offset_x), apply_offset(comp_halo.Y, offset_y)),
|
||||
comp_halo["Rvir"], zorder=10,
|
||||
linewidth=1, edgecolor=color, fill=None
|
||||
|
|
@ -228,7 +237,7 @@ def compare_halo_resolutions(
|
|||
if shared_size > best_halo_match:
|
||||
best_halo_match = shared_size
|
||||
best_halo = halo_id
|
||||
|
||||
print(f"skipped {num_skipped_for_mass} halos due to mass ratio")
|
||||
if not best_halo:
|
||||
skip_counter += 1
|
||||
continue
|
||||
|
|
@ -290,7 +299,7 @@ if __name__ == '__main__':
|
|||
comp_waveform="shannon",
|
||||
reference_resolution=128,
|
||||
comparison_resolution=256,
|
||||
plot=True,
|
||||
plot=False,
|
||||
plot3d=False,
|
||||
plot_cic=False,
|
||||
velo_halos=True,
|
||||
|
|
|
|||
|
|
@ -10,10 +10,10 @@ from pyvista import Axes
|
|||
|
||||
def main():
|
||||
rows = ["shannon", "DB8", "DB4", "DB2"]
|
||||
offset = 1.1
|
||||
offset = 2
|
||||
columns = [128, 256, 512]
|
||||
fig: Figure = plt.figure(figsize=(9, 9))
|
||||
axes: List[List[Axes]] = fig.subplots(4, 3, sharex=True, sharey=True)
|
||||
axes: List[List[Axes]] = fig.subplots(len(rows), len(columns), sharex=True, sharey=True)
|
||||
with h5py.File("vis.cache.hdf5") as vis_out:
|
||||
vmin, vmax = vis_out["vmin_vmax"]
|
||||
print(vmin, vmax)
|
||||
|
|
@ -27,7 +27,7 @@ def main():
|
|||
vmax_scaled = (vmax + offset) * mass
|
||||
rho = (rho + offset) * mass
|
||||
|
||||
img = ax.imshow(rho, norm=LogNorm(vmin=vmin_scaled, vmax=vmax_scaled), extent=extent)
|
||||
img = ax.imshow(rho.T, norm=LogNorm(vmin=vmin_scaled, vmax=vmax_scaled), extent=extent,origin="lower")
|
||||
print(img)
|
||||
pad = 5
|
||||
# based on https://stackoverflow.com/a/25814386/4398037
|
||||
|
|
|
|||
27
halo_vis.py
27
halo_vis.py
|
|
@ -1,3 +1,5 @@
|
|||
from typing import Tuple
|
||||
|
||||
import h5py
|
||||
import numpy as np
|
||||
import pandas as pd
|
||||
|
|
@ -11,8 +13,10 @@ from readfiles import read_file
|
|||
show_unbound = False
|
||||
all_in_area = True
|
||||
|
||||
Coords = Tuple[float, float, float, float] # radius, X, Y, Z
|
||||
|
||||
def load_halo_data(waveform: str, resolution: int, halo_id: int, radius=None, X=None, Y=None, Z=None):
|
||||
|
||||
def load_halo_data(waveform: str, resolution: int, halo_id: int, coords: Coords):
|
||||
dir = base_dir / f"{waveform}_{resolution}_100"
|
||||
df, meta = read_file(dir)
|
||||
df_halo, halo_lookup, unbound = read_velo_halos(dir, recursivly=False, skip_unbound=not show_unbound)
|
||||
|
|
@ -22,11 +26,14 @@ def load_halo_data(waveform: str, resolution: int, halo_id: int, radius=None, X=
|
|||
|
||||
halo = df_halo.loc[halo_id]
|
||||
if all_in_area:
|
||||
if not (radius and X and Y):
|
||||
if coords:
|
||||
radius, X, Y, Z = coords
|
||||
else:
|
||||
radius = halo["R_size"]
|
||||
X = halo["Xc"]
|
||||
Y = halo["Yc"]
|
||||
Z = halo["Zc"]
|
||||
coords: Coords = radius, X, Y, Z
|
||||
df = df[df["X"].between(X - radius, X + radius)]
|
||||
df = df[df["Y"].between(Y - radius, Y + radius)]
|
||||
halo_particles = df[df["Z"].between(Z - radius, Z + radius)]
|
||||
|
|
@ -35,7 +42,7 @@ def load_halo_data(waveform: str, resolution: int, halo_id: int, radius=None, X=
|
|||
del halo_lookup
|
||||
del unbound
|
||||
halo_particles = df.loc[list(halo_particle_ids)]
|
||||
return halo, halo_particles, meta #TODO: return and keep r,XYZ
|
||||
return halo, halo_particles, meta, coords
|
||||
|
||||
|
||||
def get_comp_id(ref_waveform: str, reference_resolution: int, comp_waveform: str, comp_resolution: int):
|
||||
|
|
@ -69,7 +76,7 @@ def main():
|
|||
rhos = {}
|
||||
ref_waveform = "shannon"
|
||||
ref_resolution = 128
|
||||
radius = None
|
||||
coords = None
|
||||
vmin = np.Inf
|
||||
vmax = -np.Inf
|
||||
with h5py.File("vis.cache.hdf5", "w") as vis_out:
|
||||
|
|
@ -83,17 +90,17 @@ def main():
|
|||
else:
|
||||
halo_id = map_halo_id(initial_halo_id, ref_waveform, ref_resolution, waveform, resolution)
|
||||
|
||||
halo, halo_particles, meta = load_halo_data(waveform, resolution, halo_id)
|
||||
if not radius:
|
||||
radius = halo["R_size"]
|
||||
X = halo["Xc"]
|
||||
Y = halo["Yc"]
|
||||
halo, halo_particles, meta, image_coords = load_halo_data(waveform, resolution, halo_id, coords)
|
||||
if not coords:
|
||||
coords = image_coords
|
||||
print(coords)
|
||||
print("mass", halo["Mvir"])
|
||||
# print("sleep")
|
||||
# sleep(100)
|
||||
radius, X, Y, Z = coords
|
||||
rho, extent = cic_from_radius(
|
||||
halo_particles.X.to_numpy(), halo_particles.Y.to_numpy(),
|
||||
1000, X, Y, radius)
|
||||
1000, X, Y, radius, periodic=False)
|
||||
rhos[(waveform, resolution)] = rho
|
||||
vmin = min(rho.min(), vmin)
|
||||
vmax = max(rho.max(), vmax)
|
||||
|
|
|
|||
62
poetry.lock
generated
62
poetry.lock
generated
|
|
@ -6,6 +6,17 @@ category = "main"
|
|||
optional = false
|
||||
python-versions = "*"
|
||||
|
||||
[[package]]
|
||||
name = "beniget"
|
||||
version = "0.4.1"
|
||||
description = "Extract semantic information about static Python code"
|
||||
category = "main"
|
||||
optional = false
|
||||
python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*"
|
||||
|
||||
[package.dependencies]
|
||||
gast = ">=0.5.0,<0.6.0"
|
||||
|
||||
[[package]]
|
||||
name = "cycler"
|
||||
version = "0.11.0"
|
||||
|
|
@ -44,6 +55,14 @@ ufo = ["fs (>=2.2.0,<3)"]
|
|||
unicode = ["unicodedata2 (>=14.0.0)"]
|
||||
woff = ["zopfli (>=0.1.4)", "brotlicffi (>=0.8.0)", "brotli (>=1.0.1)"]
|
||||
|
||||
[[package]]
|
||||
name = "gast"
|
||||
version = "0.5.3"
|
||||
description = "Python AST that abstracts the underlying Python version"
|
||||
category = "main"
|
||||
optional = false
|
||||
python-versions = ">=2.7, !=3.0.*, !=3.1.*, !=3.2.*, !=3.3.*"
|
||||
|
||||
[[package]]
|
||||
name = "h5py"
|
||||
version = "3.6.0"
|
||||
|
|
@ -196,6 +215,14 @@ python-versions = ">=3.7"
|
|||
docs = ["olefile", "sphinx (>=2.4)", "sphinx-copybutton", "sphinx-issues (>=3.0.1)", "sphinx-removed-in", "sphinx-rtd-theme (>=1.0)", "sphinxext-opengraph"]
|
||||
tests = ["check-manifest", "coverage", "defusedxml", "markdown2", "olefile", "packaging", "pyroma", "pytest", "pytest-cov", "pytest-timeout"]
|
||||
|
||||
[[package]]
|
||||
name = "ply"
|
||||
version = "3.11"
|
||||
description = "Python Lex & Yacc"
|
||||
category = "main"
|
||||
optional = false
|
||||
python-versions = "*"
|
||||
|
||||
[[package]]
|
||||
name = "posix-ipc"
|
||||
version = "1.0.5"
|
||||
|
|
@ -247,6 +274,23 @@ python-versions = "!=3.0.*,!=3.1.*,!=3.2.*,>=2.7"
|
|||
[package.dependencies]
|
||||
six = ">=1.5"
|
||||
|
||||
[[package]]
|
||||
name = "pythran"
|
||||
version = "0.11.0"
|
||||
description = "Ahead of Time compiler for numeric kernels"
|
||||
category = "main"
|
||||
optional = false
|
||||
python-versions = "*"
|
||||
|
||||
[package.dependencies]
|
||||
beniget = ">=0.4.0,<0.5.0"
|
||||
gast = ">=0.5.0,<0.6.0"
|
||||
numpy = "*"
|
||||
ply = ">=3.4"
|
||||
|
||||
[package.extras]
|
||||
doc = ["numpy", "nbsphinx", "scipy", "guzzle-sphinx-theme"]
|
||||
|
||||
[[package]]
|
||||
name = "pytz"
|
||||
version = "2022.1"
|
||||
|
|
@ -378,13 +422,17 @@ url = "https://github.com/pyvista/pyvista-wheels/raw/main/vtk-9.1.0.dev0-cp310-c
|
|||
[metadata]
|
||||
lock-version = "1.1"
|
||||
python-versions = "^3.9,<3.11"
|
||||
content-hash = "86d7667a37ddf1e9ba04e199b8a8ba3b2e14fe4b200c1e90edd320cbdbe1fb12"
|
||||
content-hash = "db1d65073afc49ea08202b1a47d7d0747b8db6c2c7235e2b771d6345ba3103fe"
|
||||
|
||||
[metadata.files]
|
||||
appdirs = [
|
||||
{file = "appdirs-1.4.4-py2.py3-none-any.whl", hash = "sha256:a841dacd6b99318a741b166adb07e19ee71a274450e68237b4650ca1055ab128"},
|
||||
{file = "appdirs-1.4.4.tar.gz", hash = "sha256:7d5d0167b2b1ba821647616af46a749d1c653740dd0d2415100fe26e27afdf41"},
|
||||
]
|
||||
beniget = [
|
||||
{file = "beniget-0.4.1-py3-none-any.whl", hash = "sha256:cb061256631313f9d06031b824f7f403baecaf609b2d3d14d43f23356cf143f2"},
|
||||
{file = "beniget-0.4.1.tar.gz", hash = "sha256:75554b3b8ad0553ce2f607627dad3d95c60c441189875b98e097528f8e23ac0c"},
|
||||
]
|
||||
cycler = [
|
||||
{file = "cycler-0.11.0-py3-none-any.whl", hash = "sha256:3a27e95f763a428a739d2add979fa7494c912a32c17c4c38c4d5f082cad165a3"},
|
||||
{file = "cycler-0.11.0.tar.gz", hash = "sha256:9c87405839a19696e837b3b818fed3f5f69f16f1eec1a1ad77e043dcea9c772f"},
|
||||
|
|
@ -431,6 +479,10 @@ fonttools = [
|
|||
{file = "fonttools-4.33.3-py3-none-any.whl", hash = "sha256:f829c579a8678fa939a1d9e9894d01941db869de44390adb49ce67055a06cc2a"},
|
||||
{file = "fonttools-4.33.3.zip", hash = "sha256:c0fdcfa8ceebd7c1b2021240bd46ef77aa8e7408cf10434be55df52384865f8e"},
|
||||
]
|
||||
gast = [
|
||||
{file = "gast-0.5.3-py3-none-any.whl", hash = "sha256:211aac1e58c167b25d3504998f2db694454a24bb1fb1225bce99420166f21d6a"},
|
||||
{file = "gast-0.5.3.tar.gz", hash = "sha256:cfbea25820e653af9c7d1807f659ce0a0a9c64f2439421a7bba4f0983f532dea"},
|
||||
]
|
||||
h5py = [
|
||||
{file = "h5py-3.6.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:a5320837c60870911645e9a935099bdb2be6a786fcf0dac5c860f3b679e2de55"},
|
||||
{file = "h5py-3.6.0-cp310-cp310-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:98646e659bf8591a2177e12a4461dced2cad72da0ba4247643fd118db88880d2"},
|
||||
|
|
@ -720,6 +772,10 @@ pillow = [
|
|||
{file = "Pillow-9.1.0-pp38-pypy38_pp73-win_amd64.whl", hash = "sha256:8d79c6f468215d1a8415aa53d9868a6b40c4682165b8cb62a221b1baa47db458"},
|
||||
{file = "Pillow-9.1.0.tar.gz", hash = "sha256:f401ed2bbb155e1ade150ccc63db1a4f6c1909d3d378f7d1235a44e90d75fb97"},
|
||||
]
|
||||
ply = [
|
||||
{file = "ply-3.11-py2.py3-none-any.whl", hash = "sha256:096f9b8350b65ebd2fd1346b12452efe5b9607f7482813ffca50c22722a807ce"},
|
||||
{file = "ply-3.11.tar.gz", hash = "sha256:00c7c1aaa88358b9c765b6d3000c6eec0ba42abca5351b095321aef446081da3"},
|
||||
]
|
||||
posix-ipc = [
|
||||
{file = "posix_ipc-1.0.5-cp27-cp27m-macosx_10_9_x86_64.whl", hash = "sha256:ccb36ba90efec56a1796f1566eee9561f355a4f45babbc4d18ac46fb2d0b246b"},
|
||||
{file = "posix_ipc-1.0.5-cp36-cp36m-macosx_10_6_intel.whl", hash = "sha256:613bf1afe90e84c06255ec1a6f52c9b24062492de66e5f0dbe068adf67fc3454"},
|
||||
|
|
@ -757,6 +813,10 @@ python-dateutil = [
|
|||
{file = "python-dateutil-2.8.2.tar.gz", hash = "sha256:0123cacc1627ae19ddf3c27a5de5bd67ee4586fbdd6440d9748f8abb483d3e86"},
|
||||
{file = "python_dateutil-2.8.2-py2.py3-none-any.whl", hash = "sha256:961d03dc3453ebbc59dbdea9e4e11c5651520a876d0f4db161e8674aae935da9"},
|
||||
]
|
||||
pythran = [
|
||||
{file = "pythran-0.11.0-py3-none-any.whl", hash = "sha256:8fea162e6096c6a3d45613a9a5389db739bf78805e713f853c9fcc7113ab9369"},
|
||||
{file = "pythran-0.11.0.tar.gz", hash = "sha256:0b2cba712e09f7630879dff69f268460bfe34a6d6000451b47d598558a92a875"},
|
||||
]
|
||||
pytz = [
|
||||
{file = "pytz-2022.1-py2.py3-none-any.whl", hash = "sha256:e68985985296d9a66a881eb3193b0906246245294a881e7c8afe623866ac6a5c"},
|
||||
{file = "pytz-2022.1.tar.gz", hash = "sha256:1e760e2fe6a8163bc0b3d9a19c4f84342afa0a2affebfaa84b01b978a02ecaa7"},
|
||||
|
|
|
|||
|
|
@ -16,6 +16,7 @@ scipy = "^1.8.0"
|
|||
pynbody = "^1.1.0"
|
||||
numba = "^0.55.1"
|
||||
tables = "^3.7.0"
|
||||
pythran = "^0.11.0"
|
||||
|
||||
[tool.poetry.dev-dependencies]
|
||||
|
||||
|
|
|
|||
4
sizes.py
4
sizes.py
|
|
@ -17,8 +17,8 @@ with pd.option_context('display.max_rows', None):
|
|||
fig: Figure = plt.figure()
|
||||
ax: Axes = fig.gca()
|
||||
|
||||
x_col = "ref_cNFW"
|
||||
y_col = "comp_cNFW"
|
||||
x_col = "ref_Mvir"
|
||||
y_col = "comp_Mvir"
|
||||
|
||||
# ax.scatter(df["ref_sizes"], df["comp_sizes"], s=1, alpha=.3)
|
||||
ax.scatter(df[x_col], df[y_col], s=1, alpha=.3)
|
||||
|
|
|
|||
Loading…
Reference in a new issue