2022-05-24 17:06:49 +02:00
|
|
|
from pathlib import Path
|
2022-05-09 15:20:10 +02:00
|
|
|
|
2022-06-10 11:06:32 +02:00
|
|
|
import numpy as np
|
2022-05-04 13:42:57 +02:00
|
|
|
import pandas as pd
|
|
|
|
from matplotlib import pyplot as plt
|
|
|
|
from matplotlib.axes import Axes
|
2022-07-18 19:27:56 +02:00
|
|
|
from matplotlib.collections import QuadMesh
|
2022-06-10 11:06:32 +02:00
|
|
|
from matplotlib.colors import LogNorm
|
2022-05-04 13:42:57 +02:00
|
|
|
from matplotlib.figure import Figure
|
|
|
|
|
2022-07-12 16:09:52 +02:00
|
|
|
# density like in Vr:
|
2022-07-18 19:27:56 +02:00
|
|
|
from halo_vis import get_comp_id
|
|
|
|
from paths import base_dir
|
2022-07-21 12:34:57 +02:00
|
|
|
from utils import figsize_from_page_fraction, rowcolumn_labels, waveforms
|
2022-07-18 19:27:56 +02:00
|
|
|
|
2022-07-12 16:09:52 +02:00
|
|
|
G = 43.022682 # in Mpc (km/s)^2 / (10^10 Msun)
|
|
|
|
|
|
|
|
|
2022-07-12 16:16:00 +02:00
|
|
|
def concentration(row, halo_type: str):
|
|
|
|
r_200crit = row[f'{halo_type}_R_200crit']
|
|
|
|
if r_200crit <= 0:
|
|
|
|
cnfw = -1
|
|
|
|
colour = 'orange'
|
|
|
|
return cnfw, colour
|
2022-07-12 16:09:52 +02:00
|
|
|
|
2022-07-12 16:16:00 +02:00
|
|
|
r_size = row[f'{halo_type}_R_size'] # largest difference from center of mass to any halo particle
|
|
|
|
m_200crit = row[f'{halo_type}_Mass_200crit']
|
|
|
|
vmax = row[f'{halo_type}_Vmax'] # largest velocity coming from enclosed mass profile calculation
|
|
|
|
rmax = row[f'{halo_type}_Rmax']
|
|
|
|
npart = row[f'{halo_type}_npart']
|
|
|
|
VmaxVvir2 = vmax ** 2 * r_200crit / (G * m_200crit)
|
|
|
|
if VmaxVvir2 <= 1.05:
|
|
|
|
if m_200crit == 0:
|
|
|
|
cnfw = r_size / rmax
|
2022-07-14 13:16:59 +02:00
|
|
|
colour = 'white'
|
2022-07-12 16:16:00 +02:00
|
|
|
else:
|
|
|
|
cnfw = r_200crit / rmax
|
2022-07-14 13:16:59 +02:00
|
|
|
colour = 'white'
|
2022-07-12 16:16:00 +02:00
|
|
|
else:
|
|
|
|
if npart >= 100: # only calculate cnfw for groups with more than 100 particles
|
|
|
|
cnfw = row[f'{halo_type}_cNFW']
|
|
|
|
colour = 'black'
|
2022-07-12 16:09:52 +02:00
|
|
|
else:
|
2022-07-12 16:16:00 +02:00
|
|
|
if m_200crit == 0:
|
|
|
|
cnfw = r_size / rmax
|
2022-07-14 13:16:59 +02:00
|
|
|
colour = 'white'
|
2022-07-12 16:09:52 +02:00
|
|
|
else:
|
2022-07-12 16:16:00 +02:00
|
|
|
cnfw = r_200crit / rmax
|
2022-07-14 13:16:59 +02:00
|
|
|
colour = 'white'
|
2022-07-12 16:16:00 +02:00
|
|
|
assert np.isclose(cnfw, row[f'{halo_type}_cNFW'])
|
|
|
|
|
2022-07-12 16:09:52 +02:00
|
|
|
return cnfw, colour
|
|
|
|
|
|
|
|
|
2022-07-21 12:34:57 +02:00
|
|
|
def plot_comparison_hist2d(ax: Axes, ax_scatter: Axes, file: Path, property: str, mode: str):
|
2022-07-08 13:43:01 +02:00
|
|
|
print("WARNING: Can only plot hist2d of properties with comp_ or ref_ right now!")
|
|
|
|
print(f" Selected property: {property}")
|
|
|
|
x_col = f"ref_{property}"
|
|
|
|
y_col = f"comp_{property}"
|
|
|
|
df = pd.read_csv(file)
|
|
|
|
if mode == 'concentration_analysis':
|
|
|
|
min_x = min([min(df[x_col]), min(df[y_col])])
|
|
|
|
max_x = max([max(df[x_col]), max(df[y_col])])
|
|
|
|
df = df.loc[2 * df.ref_cNFW < df.comp_cNFW]
|
|
|
|
else:
|
|
|
|
min_x = min([min(df[x_col]), min(df[y_col])])
|
|
|
|
max_x = max([max(df[x_col]), max(df[y_col])])
|
2022-07-21 12:34:57 +02:00
|
|
|
num_bins = 100
|
|
|
|
bins = np.geomspace(min_x, max_x, num_bins)
|
2022-07-12 16:12:34 +02:00
|
|
|
if mode == "concentration_bla" and property == 'cNFW':
|
2022-07-12 16:09:52 +02:00
|
|
|
colors = []
|
2022-07-12 16:12:34 +02:00
|
|
|
for i, row in df.iterrows():
|
2022-07-18 19:27:56 +02:00
|
|
|
comp_cnfw, comp_colour = concentration(row, halo_type="comp") # ref or comp
|
2022-07-14 13:16:59 +02:00
|
|
|
ref_cnfw, ref_colour = concentration(row, halo_type='ref')
|
|
|
|
if comp_colour == 'white' or ref_colour == 'white':
|
|
|
|
colors.append('white')
|
|
|
|
else:
|
|
|
|
colors.append('black')
|
2022-07-12 16:09:52 +02:00
|
|
|
ax.scatter(df[x_col], df[y_col], c=colors, s=1, alpha=.3)
|
|
|
|
else:
|
2022-07-21 12:34:57 +02:00
|
|
|
stds = []
|
|
|
|
for rep_row in range(num_bins):
|
|
|
|
rep_x_left = bins[rep_row]
|
|
|
|
rep_x_right = bins[rep_row] + 1
|
|
|
|
rep_bin = (rep_x_left < df[x_col]) & (df[x_col] < rep_x_right)
|
|
|
|
rep_values = df.loc[rep_bin][y_col]
|
|
|
|
if len(rep_values) > 30:
|
|
|
|
mean = rep_values.mean()
|
|
|
|
std = rep_values.std()
|
|
|
|
stds.append(std)
|
|
|
|
else:
|
|
|
|
stds.append(np.nan)
|
|
|
|
ax_scatter.step(bins, stds, label=f"{file.stem}")
|
|
|
|
|
2022-07-18 19:27:56 +02:00
|
|
|
image: QuadMesh
|
2022-07-21 12:34:57 +02:00
|
|
|
_, _, _, image = ax.hist2d(df[x_col], df[y_col], bins=(bins, bins), norm=LogNorm())
|
|
|
|
# ax.plot([rep_x_left, rep_x_left], [mean - std, mean + std], c="C1")
|
|
|
|
# ax.annotate(
|
|
|
|
# text=f"std={std:.2f}", xy=(rep_x_left, mean + std),
|
|
|
|
# textcoords="axes fraction", xytext=(0.1, 0.9),
|
|
|
|
# arrowprops={}
|
|
|
|
# )
|
2022-07-18 19:27:56 +02:00
|
|
|
print(mean - std, mean + std)
|
|
|
|
print("vmin/vmax", image.norm.vmin, image.norm.vmax)
|
|
|
|
# fig.colorbar(hist)
|
2022-07-08 13:43:01 +02:00
|
|
|
|
|
|
|
# ax.set_xscale("log")
|
|
|
|
# ax.set_yscale("log")
|
|
|
|
|
|
|
|
ax.loglog([min_x, max_x], [min_x, max_x], linewidth=1, color="C2")
|
2022-07-18 19:27:56 +02:00
|
|
|
return x_col, y_col
|
|
|
|
# ax.set_title(file.name)
|
2022-07-08 13:43:01 +02:00
|
|
|
# fig.savefig(Path(f"~/tmp/comparison_{file.stem}.pdf").expanduser())
|
2022-07-18 19:27:56 +02:00
|
|
|
# fig.suptitle
|
2022-07-08 13:43:01 +02:00
|
|
|
|
2022-07-12 16:09:52 +02:00
|
|
|
|
|
|
|
def plot_comparison_hist(file: Path, property: str, mode: str):
|
2022-07-08 13:43:01 +02:00
|
|
|
print("WARNING: Can only plot hist of properties w/o comp_ or ref_ right now!")
|
|
|
|
print(f" Selected property: {property}")
|
|
|
|
df = pd.read_csv(file)
|
|
|
|
if mode == 'concentration_analysis':
|
|
|
|
df = df.loc[2 * df.ref_cNFW < df.comp_cNFW]
|
|
|
|
fig2: Figure = plt.figure()
|
|
|
|
ax2: Axes = fig2.gca()
|
|
|
|
|
|
|
|
ax2.hist(df[property][df[property] < 50], bins=100)
|
|
|
|
ax2.set_xlabel(property)
|
|
|
|
ax2.set_title(file.name)
|
2022-07-18 19:27:56 +02:00
|
|
|
ax.set_aspect("scaled")
|
2022-07-08 13:43:01 +02:00
|
|
|
# fig2.savefig(Path(f"~/tmp/distances_{file.stem}.pdf").expanduser())
|
|
|
|
fig2.suptitle
|
|
|
|
plt.show()
|
2022-05-06 13:23:31 +02:00
|
|
|
|
2022-07-12 16:09:52 +02:00
|
|
|
|
2022-07-18 19:27:56 +02:00
|
|
|
comparisons_dir = base_dir / "comparisons"
|
2022-07-12 15:55:43 +02:00
|
|
|
|
2022-07-18 19:27:56 +02:00
|
|
|
# properties = ['group_size', 'Mass_200crit', 'Mass_tot', 'Mvir', 'R_200crit', 'Rvir', 'Vmax', 'cNFW', 'q',
|
|
|
|
# 's'] # Mass_FOF and cNFW_200crit don't work, rest looks normal except for cNFW
|
|
|
|
properties = ['Mvir']
|
2022-07-12 15:55:43 +02:00
|
|
|
# mode = 'concentration_analysis'
|
2022-07-18 19:27:56 +02:00
|
|
|
mode = 'normal'
|
2022-05-04 13:42:57 +02:00
|
|
|
|
2022-07-18 19:27:56 +02:00
|
|
|
comparisons = [(256, 512), (256, 1024)] # , (512, 1024)
|
2022-05-04 13:42:57 +02:00
|
|
|
|
2022-07-18 19:27:56 +02:00
|
|
|
for property in properties:
|
|
|
|
fig: Figure
|
|
|
|
fig, axes = plt.subplots(
|
|
|
|
len(waveforms), len(comparisons),
|
|
|
|
sharey="all", sharex="all",
|
|
|
|
figsize=figsize_from_page_fraction(columns=2)
|
|
|
|
)
|
2022-07-21 12:34:57 +02:00
|
|
|
fig_scatter: Figure = plt.figure(figsize=figsize_from_page_fraction())
|
|
|
|
ax_scatter: Axes = fig_scatter.gca()
|
|
|
|
ax_scatter.set_xscale("log")
|
2022-07-18 19:27:56 +02:00
|
|
|
for i, waveform in enumerate(waveforms):
|
|
|
|
for j, (ref_res, comp_res) in enumerate(comparisons):
|
|
|
|
file_id = get_comp_id(waveform, ref_res, waveform, comp_res)
|
|
|
|
file = comparisons_dir / file_id
|
|
|
|
print(file)
|
|
|
|
ax: Axes = axes[i, j]
|
2022-07-21 12:34:57 +02:00
|
|
|
x_col, y_col = plot_comparison_hist2d(ax, ax_scatter, file, property, mode)
|
2022-07-18 19:27:56 +02:00
|
|
|
if i == len(waveforms) - 1:
|
|
|
|
ax.set_xlabel(x_col)
|
|
|
|
if j == 0:
|
|
|
|
ax.set_ylabel(y_col)
|
|
|
|
pad = 5
|
|
|
|
rowcolumn_labels(axes, comparisons, isrow=False)
|
|
|
|
rowcolumn_labels(axes, waveforms, isrow=True)
|
|
|
|
|
|
|
|
fig.tight_layout()
|
|
|
|
fig.savefig(Path(f"~/tmp/comparison_{property}.pdf").expanduser())
|
2022-07-21 12:34:57 +02:00
|
|
|
ax_scatter.legend()
|
|
|
|
fig_scatter.tight_layout()
|
2022-07-18 19:27:56 +02:00
|
|
|
plt.show()
|
2022-07-08 13:43:01 +02:00
|
|
|
# axis_ratios = ['q', 's'] #they look normal
|
2022-05-24 17:06:49 +02:00
|
|
|
|
2022-07-08 13:43:01 +02:00
|
|
|
# for property in axis_ratios:
|
|
|
|
# plot_comparison_hist2d(file, property, 'no')
|
|
|
|
# plot_comparison_hist2d(file, property, mode)
|
2022-06-30 18:14:20 +02:00
|
|
|
|
2022-07-08 13:43:01 +02:00
|
|
|
# plot_comparison_hist2d(file, 'cNFW_200mean', mode)
|
2022-06-10 11:06:32 +02:00
|
|
|
|
2022-07-12 15:55:43 +02:00
|
|
|
# ref_property = 'ref_cNFW_200crit'
|
|
|
|
# comp_property = 'comp_cNFW_200crit'
|
|
|
|
|
|
|
|
# df = pd.read_csv(file)
|
|
|
|
# all_ref_structure_types: pd.DataFrame = df[ref_property]
|
|
|
|
# all_comp_structure_types: pd.DataFrame = df[comp_property]
|
|
|
|
|
|
|
|
# df_odd: pd.DataFrame = df.loc[2 * df.ref_cNFW < df.comp_cNFW]
|
|
|
|
# odd_ref_structure_types: pd.DataFrame = df_odd[ref_property]
|
|
|
|
# odd_comp_structure_types: pd.DataFrame = df_odd[comp_property]
|
|
|
|
|
|
|
|
# print(all_ref_structure_types.mean(), all_comp_structure_types.mean())
|
|
|
|
# print(odd_ref_structure_types.mean(), odd_comp_structure_types.mean())
|
|
|
|
|
2022-05-24 17:06:49 +02:00
|
|
|
|
2022-07-12 16:09:52 +02:00
|
|
|
# ref_colour = []
|
|
|
|
# comp_colour = []
|
|
|
|
# ref_cnfw = []
|
|
|
|
# comp_cnfw = []
|
|
|
|
# df = pd.read_csv(file)
|
|
|
|
#
|
|
|
|
# for index, row in df.iterrows():
|
|
|
|
# cnfw, colour = concentration(row)
|
|
|
|
# ref_cnfw.append(cnfw[0])
|
|
|
|
# ref_colour.append(colour[0])
|
|
|
|
# comp_cnfw.append(cnfw[1])
|
|
|
|
# comp_colour.append(colour[1])
|
|
|
|
#
|
|
|
|
# fig: Figure = plt.figure()
|
|
|
|
# ax: Axes = fig.gca()
|
|
|
|
#
|
|
|
|
# ax.scatter(ref_cnfw, comp_cnfw, s=1, c=comp_colour, alpha=.3)
|
|
|
|
# ax.set_xscale("log")
|
|
|
|
# ax.set_yscale("log")
|
|
|
|
# plt.show()
|
2022-05-04 13:42:57 +02:00
|
|
|
|
2022-07-08 13:43:01 +02:00
|
|
|
# #Maybe for later:
|
|
|
|
# if __name__ == '__main__':
|
|
|
|
# print('Run with sizes.py <Path to file> <property: str> <mode: str>')
|
|
|
|
# file = Path(argv[1])
|
|
|
|
# property = str(argv[2])
|
|
|
|
# mode = str(argv[3])
|
2022-05-09 15:20:10 +02:00
|
|
|
|
|
|
|
|
2022-07-12 16:09:52 +02:00
|
|
|
# #This is to find the median of the quality of our matches
|
|
|
|
# matches:pd.DataFrame=df["match"]
|
|
|
|
# print(matches)
|
|
|
|
# exit()
|
|
|
|
# print(matches.median())
|
|
|
|
# print(matches.std())
|
|
|
|
# exit()
|
2022-05-09 15:20:10 +02:00
|
|
|
|
2022-07-12 16:09:52 +02:00
|
|
|
# #This is to save weird concentration data to own csv
|
|
|
|
# df_odd: pd.DataFrame = df.loc[2 * df.ref_cNFW < df.comp_cNFW]
|
|
|
|
# df_odd.to_csv("weird_cnfw.csv")
|
|
|
|
# exit()
|