mirror of
https://github.com/Richard-Sti/csiborgtools.git
synced 2024-12-22 22:38:03 +00:00
9e4b34f579
* Update README * Update density field reader * Update name of SDSSxALFAFA * Fix quick bug * Add little fixes * Update README * Put back fit_init * Add paths to initial snapshots * Add export * Remove some choices * Edit README * Add Jens' comments * Organize imports * Rename snapshot * Add additional print statement * Add paths to initial snapshots * Add masses to the initial files * Add normalization * Edit README * Update README * Fix bug in CSiBORG1 so that does not read fof_00001 * Edit README * Edit README * Overwrite comments * Add paths to init lag * Fix Quijote path * Add lagpatch * Edit submits * Update README * Fix numpy int problem * Update README * Add a flag to keep the snapshots open when fitting * Add a flag to keep snapshots open * Comment out some path issue * Keep snapshots open * Access directly snasphot * Add lagpatch for CSiBORG2 * Add treatment of x-z coordinates flipping * Add radial velocity field loader * Update README * Add lagpatch to Quijote * Fix typo * Add setter * Fix typo * Update README * Add output halo cat as ASCII * Add import * Add halo plot * Update README * Add evaluating field at radial distanfe * Add field shell evaluation * Add enclosed mass computation * Add BORG2 import * Add BORG boxsize * Add BORG paths * Edit run * Add BORG2 overdensity field * Add bulk flow clauclation * Update README * Add new plots * Add nbs * Edit paper * Update plotting * Fix overlap paths to contain simname * Add normalization of positions * Add default paths to CSiBORG1 * Add overlap path simname * Fix little things * Add CSiBORG2 catalogue * Update README * Add import * Add TNG density field constructor * Add TNG density * Add draft of calculating BORG ACL * Fix bug * Add ACL of enclosed density * Add nmean acl * Add galaxy bias calculation * Add BORG acl notebook * Add enclosed mass calculation * Add TNG300-1 dir * Add TNG300 and BORG1 dir * Update nb
829 KiB
829 KiB
In [1]:
import sys
import numpy as np
import h5py
import matplotlib.pyplot as plt
import scienceplots
from glob import glob
sys.path.append("../")
import csiborgtools
%matplotlib widget
%load_ext autoreload
%autoreload 2
In [5]:
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
files = glob("/mnt/extraspace/rstiskalek/CSiBORG/overlap/overlap_0*")
nsim0 = 7444
files = [f for f in files if str(nsim0) in f]
nsimxs = np.sort([int(f.split("_")[-1].split(".")[0][1:]) for f in files])
In [7]:
csiborgtools.read.get_cross_sims(7444, paths, True)
Out[7]:
In [ ]:
In [3]:
cat0 = csiborgtools.read.HaloCatalogue(nsim0, paths, minmass=("totpartmass", 1e12))
catxs = [csiborgtools.read.HaloCatalogue(nsimx, paths, minmass=("totpartmass", 1e12)) for nsimx in nsimxs]
In [4]:
reader = csiborgtools.read.NPairsOverlap(cat0, catxs, paths)
In [5]:
x = reader.cat0("totpartmass")
y1 = reader.summed_overlap(False)
In [6]:
y2 = reader.prob_nomatch(False)
In [8]:
plt.figure()
plt.scatter(1 - np.mean(y1, axis=1), np.mean(y2, axis=1), c=np.log10(x), s=2)
plt.colorbar(label=r"$\log_{10} M_{\rm halo} / M_\odot$")
t = np.linspace(0.3, 1, 100)
plt.plot(t, t, color="red", linestyle="--")
plt.xlabel(r"$1 - \langle \mathcal{O}_a^{\mathcal{A} \mathcal{B}} \rangle_{\mathcal{B}}$")
plt.ylabel(r"$\langle \eta_a^{\mathcal{A} \mathcal{B}} \rangle_{\mathcal{B}}$")
plt.tight_layout()
plt.savefig("../plots/prob_nomatch_vs_overlap_nosmooth.png", dpi=450)
plt.show()
In [ ]:
In [ ]:
In [9]:
plt.close("all")
ks = [1, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100]
for k in ks:
plt.figure()
plt.title(r"Crossing with ${}$ simulations.".format(k))
plt.hexbin(x, np.mean(y1[:, :k], axis=1), mincnt=1, xscale="log", bins="log", gridsize=50)
plt.colorbar()
plt.xlabel(r"$M_{\rm tot} / M_\odot$")
plt.ylabel("Mean overlap")
plt.tight_layout()
plt.ylim(0, 1)
# plt.savefig("../plots/overlap_std_{}.png".format(k), dpi=300, bbox_inches="tight")
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [136]:
x = reader.cat0("totpartmass")
mu, std = reader.counterpart_mass(True, 0.1, return_full=False)
In [137]:
plt.figure()
plt.scatter(x, mu, s=1)
# plt.hexbin(x, np.mean(y, axis=1), mincnt=1, xscale="log", bins="log", gridsize=50)
t = np.linspace(1e12, 1e15, 100)
plt.plot(t, t, c="red", ls="--")
plt.xscale("log")
plt.yscale("log")
plt.xlabel(r"Reference $M_{\rm tot} / M_\odot$")
plt.ylabel(r"Cross $M_{\rm tot} / M_\odot$")
plt.tight_layout()
# plt.savefig("../plots/prob_nomatch_vs_mass.png", dpi=300, bbox_inches="tight")
plt.show()
In [126]:
plt.figure()
plt.hexbin(x, np.mean(y, axis=1), mincnt=1, xscale="log", bins="log", gridsize=50)
plt.xlabel(r"$M_{\rm tot} / M_\odot$")
plt.ylabel("Mean probability of no match")
plt.tight_layout()
plt.savefig("../plots/prob_nomatch_vs_mass.png", dpi=300, bbox_inches="tight")
plt.show()
In [ ]:
In [84]:
In [118]:
In [117]:
import imageio
images = []
for k in ks:
filename = "../plots/overlap_{}.png".format(k)
images.append(imageio.imread(filename))
imageio.mimsave("../plots/overlap_movie.gif", images, duration=0.5)
In [ ]:
In [ ]:
In [ ]:
In [21]:
In [22]:
clumpmap0 = csiborgtools.read.read_h5(paths.particles_path(nsim0))["clumpmap"]
parts0 = csiborgtools.read.read_h5(paths.initmatch_path(nsim0, "particles"))["particles"]
clid2map0 = {clid: i for i, clid in enumerate(clumpmap0[:, 0])}
# # clumpmapx = csiborgtools.read.read_h5(paths.particles_path(nsimx))["clumpmap"]
# # partsx = csiborgtools.read.read_h5(paths.initmatch_path(nsimx, "particles"))["particles"]
# # clid2mapx = {clid: i for i, clid in enumerate(clumpmapx[:, 0])}
cat0 = csiborgtools.read.HaloCatalogue(nsim0, paths, load_initial=True, minmass=None, with_lagpatch=True, load_clumps_cat=True)
catx = csiborgtools.read.HaloCatalogue(nsimx, paths, load_initial=True, minmass=None, with_lagpatch=True, load_clumps_cat=True)
In [24]:
f = np.load(paths.velocity_field_path("PCS", 7444))
In [28]:
plt.figure()
plt.imshow(np.sum(f[0, ...], axis=0))
plt.show()
In [ ]:
In [76]:
x = pair.dist(in_initial=True, norm_kind="sum_patch")
y = pair.overlap(False)
x = np.concatenate(x)
y = np.concatenate(y)
In [80]:
from scipy.special import comb
from itertools import combinations
In [81]:
ics = paths.get_ics()
In [88]:
import random
combs = list(combinations(ics, 2))
random.shuffle(combs)
In [91]:
random.Random(32).shuffle(combs)
In [92]:
combs
Out[92]:
In [ ]:
In [77]:
plt.figure()
plt.scatter(x, y, s=0.1)
plt.show()
In [ ]:
In [53]:
for i in list(clid2map0.keys())[50000:]:
print(i)
k = clid2map0[i]
X = csiborgtools.read.load_parent_particles(k, parts0, clumpmap0, clid2map0, cat0.clumps_cat)
if X is not None:
break
X
Out[53]:
In [54]:
plt.figure()
plt.scatter(X[:, 0], X[:, 1], s=1)
plt.show()
In [ ]:
In [5]:
parts = csiborgtools.read.read_h5(paths.particles_path(nsim0))["particles"]
In [6]:
box = csiborgtools.read.CSiBORGBox(nsnap0, nsim0, paths)
field = csiborgtools.field.DensityField(box, "CIC")
In [7]:
rho = np.load(paths.density_field_path("PCS", 7444, False))
rho_rsp = np.load(paths.density_field_path("PCS", 7444, True))
In [14]:
plt.figure()
plt.imshow(np.sum(rho, axis=1))
plt.show()
In [13]:
plt.figure()
plt.imshow(np.sum(rho_rsp, axis=1))
plt.show()
In [ ]:
In [25]:
np.abs(x - 1).max()
Out[25]:
In [ ]:
In [94]:
k0 = clumpmap0[clid2map0[0], 2] + 1
# parts.shape[0]
In [96]:
parts.shape[0] - k0
Out[96]:
In [44]:
clumps_cat = csiborgtools.read.ClumpsCatalogue(nsim0, paths, load_fitted=True, minmass=None, maxdist=155 / 0.705)
In [119]:
cat = csiborgtools.read.HaloCatalogue(nsim0, paths, minmass=None)
In [124]:
cat['ID']
Out[124]:
In [114]:
len(cat)
Out[114]:
In [118]:
numpy.isnan(cat["lagpatch"]).sum()
Out[118]:
In [ ]:
In [88]:
cat['x0']
Out[88]:
In [ ]:
In [46]:
len(clumps_cat)
Out[46]:
In [50]:
m = numpy.isfinite(clumps_cat["lambda200c"])
plt.figure()
plt.hist(np.log10(clumps_cat["lambda200c"][m]), bins="auto")
plt.show()
In [ ]:
In [ ]:
In [5]:
clumps_cat["index"][:10]
Out[5]:
In [6]:
X = csiborgtools.read.load_parent_particles(3, parts, clumpmap0, clid2map0, clumps_cat)
In [29]:
clump = csiborgtools.fits.Clump(X, clumps_cat[2], box)
clumps_cat[2]["index"]
Out[29]:
In [37]:
clump.spherical_overdensity_mass(200)
Out[37]:
In [26]:
def f(mass):
return numpy.cumsum(mass)
@jit(nopython=True)
def g(mass):
return numpy.cumsum(mass)
x = clump['M']
In [25]:
%timeit f(x)
In [28]:
%timeit g(x)
In [ ]:
In [19]:
%timeit clump.spherical_overdensity_mass(200)
In [14]:
plt.figure()
plt.plot(rs, y)
plt.yscale("log")
plt.xscale("log")
plt.show()
In [16]:
%timeit clump.r()
In [87]:
x, y, z = clump['x'], clump['y'], clump['z']
x0, y0, z0 = clumps_cat[3]['x'], clumps_cat[3]['y'], clumps_cat[3]['z']
In [93]:
%timeit ((x - x0)**2 + (y - y0)**2 + (z - z0)**2)**0.5
In [92]:
%timeit clump.r()
In [267]:
delta_bckg *= 2
In [260]:
%timeit delta1, delta2, cellmins, nonzero = overlapper.make_deltas(pos, pos2, mass, mass)
# __ = overlapper.make_deltas(pos, pos2, mass, mass)
In [271]:
%timeit overlapper(pos, pos2, mass, mass, delta_bckg)
In [259]:
fig, axs = plt.subplots(ncols=2, figsize=(10, 5))
axs[0].imshow(np.sum(delta1, axis=0))
axs[1].imshow(np.sum(delta2, axis=0))
fig.tight_layout()
In [93]:
delta = overlapper.make_bckg_delta(pos, mass / mass)
In [282]:
from functools import cache
from functools import lru_cache
In [283]:
from time import sleep
@lru_cache(maxsize=1024)
def f(x):
sleep(0.5)
return x * x
In [286]:
%timeit f(2)
In [77]:
clumpmap0[:]
Out[77]:
In [72]:
dens = overlapper.make_delta(pos, mass, subbox=True)
In [75]:
In [ ]:
In [ ]:
In [55]:
x.shape, x0.shape
Out[55]:
In [56]:
plt.figure()
# plt.scatter(x[:, 0], x[:, 1], s=0.1)
plt.scatter(x0[:, 0], x0[:, 1], s=0.1)
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [11]:
unique_clump_ids = numpy.unique(clump_ids)
In [22]:
unique_clump_ids[:100]
Out[22]:
In [26]:
%timeit numpy.nan == numpy.nan
In [27]:
%timeit numpy.isnan(numpy.nan)
In [29]:
test = {"asd": 1}
with h5py.File("test.hp", 'w') as f:
f.create_dataset("test", data=test)
In [ ]:
f
In [23]:
k0, kf = minmax_clump(259, clump_ids)
kf - k0 + 1
Out[23]:
In [24]:
clump_ids[k0:kf+1]
Out[24]:
In [ ]:
In [3]:
import numpy
nsim = nsim0
nsnap = nsnap0
part0 = h5py.File(paths.initmatch_path(nsim, "particles"), 'r')['particles']
cmap = h5py.File(paths.particles_path(nsim, "clumpmap"), 'r')
partf = h5py.File(paths.particles_path(nsim), 'r')["particles"]
In [8]:
part0[cmap[str(558655)], :]
In [5]:
cat = csiborgtools.read.ClumpsCatalogue(nsim, paths, rawdata=True)
In [44]:
from tqdm import tqdm
y = np.zeros((10000, 2), dtype=np.int32)
for i, hid in tqdm(enumerate(cat["index"][cat.ismain][:10000])):
n = np.sum(cat["parent"] == hid)
y[i, :] = hid, n
In [8]:
cat["index"][cat["parent"] == 558655]
Out[8]:
In [45]:
y[np.argmax(y[:, 1]), :]
Out[45]:
In [8]:
keys = list(cmap.keys())
In [11]:
csiborgtools.fits.load_parent_particles(558655, part0, cmap, cat)
# Xf = csiborgtools.fits.load_parent_particles(378, partf, cmap, cat)
In [6]:
plt.figure()
plt.scatter(X[:, 0], X[:, 1], s=0.1)
# plt.scatter(Xf[:, 0], Xf[:, 1], s=0.1)
plt.show()
In [ ]:
In [3]:
pid_final = np.load("pid_final.npy")
pid_init = np.load("pid_initial.npy")
In [4]:
pid_order = np.argsort(pid_final).astype(np.int32)
pid_order_2 = np.argsort(pid_order).astype(np.int32)
In [6]:
# Sort them
pid_init = pid_init[np.argsort(pid_init)]
In [11]:
np.alltrue(pid_init[pid_order_2] == pid_final)
Out[11]:
In [10]:
Out[10]:
In [ ]:
In [ ]:
In [ ]:
In [25]:
pid_final[pid_order][pid_order_2]
Out[25]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [18]:
part0['ID']
Out[18]:
In [ ]:
In [ ]:
In [8]:
indx1 = numpy.argsort(part0['ID'])
indx2 = numpy.argsort(pid)
indx1_inv = numpy.argsort(indx1)
part0 = part0[indx2][indx1_inv]
# pid = pid[indx2]
In [9]:
part0
Out[9]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [10]:
pid
Out[10]:
In [6]:
pid
Out[6]:
In [9]:
cat = csiborgtools.read.ClumpsCatalogue(nsim, paths, load_fitted=False, rawdata=True)
parent_ids = cat["index"][cat.ismain]
parent_ids[:50]
Out[9]:
In [54]:
i = cat["index"][np.argsort(cat["mass_cl"])[::-3][213]]
i
Out[54]:
In [55]:
# clid = 10170336
# clid = 10170336 + 12
clid = i
mmain_indxs = cat["index"][cat["parent"] == clid]
# mask = clump_ids == clid
# print(mask.sum())
print(mmain_indxs)
print([numpy.sum(clump_ids == j) for j in mmain_indxs])
In [56]:
mask = numpy.isin(clump_ids, mmain_indxs)
mask.sum()
Out[56]:
In [57]:
# coords = part0[mask, :]
plt.figure()
plt.scatter(part0['x'][mask], part0['y'][mask], s=1)
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [193]:
particles0 = h5py.File(paths.initmatch_path(nsim0, "particles"), 'r')['particles']
particlesx = h5py.File(paths.initmatch_path(nsimx, "particles") , 'r')['particles']
In [4]:
# delta_bckg = overlapper.make_bckg_delta(particles0, verbose=True)
# delta_bckg = overlapper.make_bckg_delta(particlesx, delta=delta_bckg, verbose=True)
# np.save("./bckg_{}_{}.npy".format(nsim0, nsimx), delta_bckg)
# delta_bckg = np.load("./bckg_{}_{}.npy".format(nsim0, nsimx))
In [194]:
map0 = h5py.File(paths.initmatch_path(nsim0, "halomap"), 'r')
mapx = h5py.File(paths.initmatch_path(nsimx, "halomap"), 'r')
In [174]:
keys = np.array(sorted([int(x) for x in mapx.keys()]))
In [175]:
paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
halocat = csiborgtools.read.HaloCatalogue(nsimx, paths, rawdata=True)
clumpcat = csiborgtools.read.ClumpsCatalogue(nsimx, paths, rawdata=True)
In [176]:
clid = 2
np.sum(clumpcat["parent"] == clid)
# mmain_indxs = clumpcat["index"][clumpcat["parent"] == clid]
# mmain_indxs
# mmain_mask = numpy.isin(clump_ids, mmain_indxs, assume_unique=True)
Out[176]:
In [195]:
fs = list(map0.keys())
In [196]:
from tqdm import tqdm
tot = 0
for f in tqdm(fs):
tot += map0[f].shape[0]
print(tot)
In [197]:
particles0
Out[197]:
In [198]:
tot / 182510395
Out[198]:
In [172]:
particlesx
Out[172]:
In [177]:
i = 10170336
k = np.where(cat["ID"] == i)[0][0]
k
Out[177]:
In [178]:
# i = np.random.choice(keys[keys<100000])
# # i = 10040155 #+ 11
print(i)
# i = 32
mask = map0[str(i)]
mask
Out[178]:
In [183]:
halo = particlesx[mask, :]
In [184]:
np.sum(halo[:, 3]) / 1.1641532e-10
Out[184]:
In [186]:
plt.figure()
plt.scatter(halo[:, 0], halo[:, 1], s=2)
# plt.scatter(cat['x'][k], cat['y'][k], s=10, c='r')
plt.xlabel(r"x")
plt.ylabel(r"y")
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [21]:
plt.figure()
plt.imshow(np.log10(np.sum(delta_bckg, axis=0)))
plt.show()
In [2]:
f = h5py.File("../data/particles_7444.h5", "r")
In [8]:
f["particles"][0, :]
Out[8]:
In [ ]: