diff --git a/README.md b/README.md index 96c6f42..18590ca 100644 --- a/README.md +++ b/README.md @@ -1,18 +1,20 @@ # CSiBORG tools -## :scroll: Short-term TODO -- [x] Make sure I am not taking halos outside of the well-resolved region. +## Short-term TODO +- [ ] Get `powmes` working and compare its power spectrum to `Pylians`. - [ ] Check alignment of the angular momentum +- [x] Make sure I am not taking halos outside of the well-resolved region. - [x] Calculate the gravitational potential. - [x] Add the gravitational field calculation. -## :hourglass: Long-term TODO +## Long-term TODO - [ ] Calculate the cross-correlation in CSiBORG. Should see the scale of the constraints? -- [x] Calculate DM environmental properties. - [ ] Find the distribution of particles in the first snapshot - [ ] Implement a custom model for matchin galaxies to halos. +- [x] Calculate DM environmental properties. -## :bulb: Open questions +## Open questions - What scaling of the search region? No reason for it to be a multiple of $R_{200c}$. +- Begin extracting the DM environmental properties at galaxy positions? diff --git a/csiborgtools/read/make_cat.py b/csiborgtools/read/make_cat.py index c380175..73d2f4d 100644 --- a/csiborgtools/read/make_cat.py +++ b/csiborgtools/read/make_cat.py @@ -228,6 +228,18 @@ class HaloCatalogue: """ return numpy.vstack([self["v{}".format(p)] for p in ("x", "y", "z")]).T + @property + def angmomentum(self): + """ + Angular momentum of halos in the box coordinate system. + + Returns + ------- + angmom : 2-dimensional array + Array of shape `(n_halos, 3)`. + """ + return numpy.vstack([self["L{}".format(p)] for p in ("x", "y", "z")]).T + def radius_neigbours(self, X, radius): """ Return sorted nearest neigbours within `radius` or `X`. diff --git a/csiborgtools/read/outsim.py b/csiborgtools/read/outsim.py index f74522d..fc7b319 100644 --- a/csiborgtools/read/outsim.py +++ b/csiborgtools/read/outsim.py @@ -19,6 +19,7 @@ I/O functions for analysing the CSiBORG realisations. import numpy from os.path import (join, dirname, basename, isfile) +import gc from os import remove from tqdm import trange from astropy.io import ascii @@ -127,6 +128,7 @@ def combine_splits(n_splits, part_reader, cols_add, remove_splits=False, def make_ascii_powmes(particles, fout, verbose=True): """ Write an ASCII file with appropriate formatting for POWMES. + This is an extremely memory inefficient implementation. Parameters ---------- @@ -154,6 +156,9 @@ def make_ascii_powmes(particles, fout, verbose=True): print("Writing temporary file `{}`...".format(ftemp)) ascii.write(out, ftemp, overwrite=True, delimiter=",", fast_writer=True) + del out + gc.collect() + # Write to the first line the number of particles if verbose: print("Writing the full file `{}`...".format(fout)) diff --git a/scripts/run_asciipos.py b/scripts/run_asciipos.py new file mode 100644 index 0000000..03af85a --- /dev/null +++ b/scripts/run_asciipos.py @@ -0,0 +1,68 @@ +# Copyright (C) 2022 Richard Stiskalek +# This program is free software; you can redistribute it and/or modify it +# under the terms of the GNU General Public License as published by the +# Free Software Foundation; either version 3 of the License, or (at your +# option) any later version. +# +# This program is distributed in the hope that it will be useful, but +# WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +# Public License for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to the Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +"""A script to dump or remove files for POWMES.""" + +from argparse import ArgumentParser +import numpy +from datetime import datetime +from os.path import join, exists +from os import remove +from mpi4py import MPI +try: + import csiborgtools +except ModuleNotFoundError: + import sys + sys.path.append("../") + import csiborgtools +import utils + +parser = ArgumentParser() +parser.add_argument("--mode", type=str, choices=["dump", "remove"]) +args = parser.parse_args() + +F64 = numpy.float64 +I64 = numpy.int64 + +# Get MPI things +comm = MPI.COMM_WORLD +rank = comm.Get_rank() +nproc = comm.Get_size() + + +dumpdir = join(utils.dumpdir, "temp_powmes") +fout = join(dumpdir, "out_{}_{}.ascii") +paths = csiborgtools.read.CSiBORGPaths() + + +n_sims = paths.ic_ids[:1] +for i in csiborgtools.fits.split_jobs(len(n_sims), nproc)[rank]: + print("{}: calculating {}th simulation.".format(datetime.now(), i)) + n_sim = n_sims[i] + n_snap = paths.get_maximum_snapshot(n_sim) + paths.set_info(n_sim, n_snap) + + f = fout.format(n_sim, n_snap) + if args.mode == "dump": + # Read the particles + reader = csiborgtools.read.ParticleReader(paths) + particles = reader.read_particle(["x", "y", "z", "M"]) + csiborgtools.read.make_ascii_powmes(particles, f, verbose=True) + else: + if exists(f): + remove(f) + +comm.Barrier() +if rank == 0: + print("All finished! See you!") diff --git a/scripts/run_fit_halos.py b/scripts/run_fit_halos.py index 8aa208c..0f65e21 100644 --- a/scripts/run_fit_halos.py +++ b/scripts/run_fit_halos.py @@ -43,6 +43,7 @@ dumpdir = utils.dumpdir loaddir = join(utils.dumpdir, "temp") cols_collect = [("npart", I64), ("totpartmass", F64), ("Rs", F64), ("vx", F64), ("vy", F64), ("vz", F64), + ("Lx", F64), ("Ly", F64), ("Lz", F64), ("rho0", F64), ("conc", F64), ("rmin", F64), ("rmax", F64), ("r200", F64), ("r500", F64), ("m200", F64), ("m500", F64), ("lambda200c", F64)] @@ -66,6 +67,7 @@ for i, n_sim in enumerate(paths.ic_ids): cols = [("index", I64), ("npart", I64), ("totpartmass", F64), ("Rs", F64), ("rho0", F64), ("conc", F64), ("lambda200c", F64), ("vx", F64), ("vy", F64), ("vz", F64), + ("Lx", F64), ("Ly", F64), ("Lz", F64), ("rmin", F64), ("rmax", F64), ("r200", F64), ("r500", F64), ("m200", F64), ("m500", F64)] out = csiborgtools.utils.cols_to_structured(N, cols) @@ -84,6 +86,7 @@ for i, n_sim in enumerate(paths.ic_ids): out["vx"][n] = numpy.average(clump.vel[:, 0], weights=clump.m) out["vy"][n] = numpy.average(clump.vel[:, 1], weights=clump.m) out["vz"][n] = numpy.average(clump.vel[:, 2], weights=clump.m) + out["Lx"][n], out["Ly"][n], out["Lz"][n] = clump.angular_momentum # Spherical overdensity radii and masses rs, ms = clump.spherical_overdensity_mass([200, 500])