csiborgtools/scripts/concentration_fit.ipynb
Richard Stiskalek 8a56c22813
Within halo work and NFW fit (#4)
* add listing of snapshots

* change distance to comoving

* ignore cp files

* rename nb

* add str to list

* add NFW profile shapes

* add fits imports

* Rename to Nsnap

* in clumps_read only select props

* make clumpid int

* expand doc

* add import

* edit readme

* distribute halos

* add profile & posterior

* add import

* add import

* add documentation

* add rvs and init guess

* update todo

* update nb

* add file

* return end index too

* change clump_ids format to int32

* skeleton of dump particle

* update nb

* add func to drop 0 clump indxs parts

* add import

* add halo dump

* switch to float32

* Update TODO

* update TODO

* add func that loads a split

* add halo object

* Rename to clump

* make post work with a clump

* add optimiser

* add Nsplits

* ignore submission scripts

* ignore .out

* add dumppath

* add job splitting

* add split halos script

* rename file

* renaem files

* rm file

* rename imports

* edit desc

* add pick clump

* add number of particles

* update TODO

* update todo

* add script

* add dumping

* change dumpdir structure

* change dumpdir

* add import

* Remove tqdm

* Increase the number of splits

* rm shuffle option

* Change to remove split

* add emojis

* fix part counts in splits

* change num of splits

* rm with particle cut

* keep splits

* fit only if 10 part and more

* add min distance

* rm warning about not set vels

* update TODO

* calculate rho0 too

* add results collection

* add import

* add func to combine splits

* update TODO

* add extract cols

* update nb

* update TODO
2022-10-30 20:16:56 +00:00

420 KiB

In [1]:
import numpy as np
%matplotlib notebook
import matplotlib.pyplot as plt
try:
    import csiborgtools
except ModuleNotFoundError:
    import sys
    sys.path.append("../")
    import csiborgtools
import utils
%load_ext autoreload
%autoreload 2

import joblib
from os.path import join
from glob import glob
In [53]:
from astropy.cosmology import FlatLambdaCDM
In [ ]:

In [ ]:

In [54]:
cosmo = FlatLambdaCDM(H0=70.5, Om0=0.307, Tcmb0=2.728)
In [60]:
x = cosmo.Om0 - 1

18*np.pi**2 + 82 * x - 39 * x**2
Out[60]:
102.09716821960843
In [ ]:

In [ ]:

In [2]:
Nsim = 9844
simpath = csiborgtools.io.get_sim_path(Nsim)
Nsnap = 1016
In [40]:
fname = join(utils.dumpdir, "ramses_out_{}_{}.npy".format(str(Nsim).zfill(5), str(Nsnap).zfill(5)))
In [41]:
data = np.load(fname)
In [49]:
plt.figure()
plt.scatter(data["logRs"], data["rho0"], s=2.5, rasterized=True)
plt.yscale("log")

plt.ylabel(r"$\rho_0$")
plt.xlabel(r"$\log R_{\rm s}$")
# plt.savefig("../plots/rho0.png", dpi=450)
plt.show()
No description has been provided for this image
In [44]:
data["rho0"]
Out[44]:
array([287721.82008412, 249035.67807171,  69665.64384864, ...,
         2582.73496889,             nan,             nan])
In [ ]:

In [ ]:

In [24]:
out = csiborgtools.io.combine_splits(utils.Nsplits, Nsim, Nsnap, utils.dumpdir, )
100%|███████████████████████████████████████████████████████████████████████████████| 200/200 [00:24<00:00,  8.16it/s]
In [27]:
out.dtype.names
Out[27]:
('index',
 'level',
 'parent',
 'ncell',
 'peak_x',
 'peak_y',
 'peak_z',
 'rho-',
 'rho+',
 'rho_av',
 'mass_cl',
 'relevance',
 'npart',
 'totpartmass',
 'logRs')
In [36]:
plt.figure()

plt.scatter(out["mass_cl"], out["logRs"], s=3, rasterized=True)
# t = np.logspace(-8, -4, 1000)
# plt.plot(t, t, c="red", ls="--")


# plt.yscale("log")
plt.xscale("log")
plt.xlabel("mass_cl")
# plt.ylabel("summed mass of all particles")

plt.savefig("../plots/mass.png", dpi=400)
plt.show()
No description has been provided for this image
In [ ]:

In [22]:
mask = np.isin(out["index"], arr["index"])
In [23]:
mask
Out[23]:
array([ True,  True,  True, ..., False, False, False])
In [ ]:

In [ ]:
np.where()
In [ ]:

In [16]:

Out[16]:
array([       1,        2,        3, ..., 21795024, 21796985, 21797711])
In [ ]:

In [ ]:

In [ ]:
# clump_ids = csiborgtools.io.read_clumpid(Nsnap, simpath)
# clumps = csiborgtools.io.read_clumps(Nsnap, simpath, )

# particles = csiborgtools.io.read_particle(["x", "y", "z", "M", "level"], Nsnap, simpath)
# clump_ids, particles = csiborgtools.io.drop_zero_indx(clump_ids, particles)

# with_particles = csiborgtools.fits.clump_with_particles(clump_ids, clumps)
# clumps = clumps[with_particles]
In [ ]:
f = join(utils.dumpdir, "ramses_out_09844_01016_123.npy")
f = np.load(f)
In [ ]:
plt.figure()
plt.hist(f["logRs"], bins="auto")
plt.show()
In [ ]: