Angular momentum (#14)

* add ang mom

* Update TODO

* add ascii dump script

* add comment

* Updating TODO

* add angmom

* add garbage collector

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?

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`.

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))

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!")

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])