Aquila-Consortium/csiborgtools
0
0
Fork 0
mirror of https://github.com/Richard-Sti/csiborgtools.git synced 2024-12-22 18:28:02 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Improve paths and documentation (#38)

Browse source 
* pep8

* style issue

* Documentation edits

* Remove old files

* Remove more old files

* Doc & stop setting snap and nsim in paths

* add nsnap, nsim support

* Remove blank space

* docs

* Docs edits

* Reduce redudant code

* docs

* Documentation

* Docs

* Simplify

* add nsnap nsim arguments

* Remove redundant docs

* Fix typos

* Shorten catalogue

* Remove import

* Remove blank line

* Docs only edits

* Rearrange imports

* Remove blank space

* Rearrange imports

* Rearrange imports

* Remove blank space

* Update submission scritp

* Edit docs

* Remove blank line

* new paths

* Update submission scripts

* Update file handling

* Remove blank line

* Move things around

* Ne paths

* Edit submission script

* edit paths

* Fix typo

* Fix bug

* Remove import

* Update nb
This commit is contained in:
Richard Stiskalek 2023-04-04 15:22:00 +01:00 committed by GitHub
parent 035d7e0071
commit 8d34b832af
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
32 changed files with 954 additions and 5358 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
csiborgtools/__init__.py
View file

@ -12,5 +12,4 @@
# 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.
from csiborgtools import (read, match, utils, units, fits, field) # noqa

1
csiborgtools/field/__init__.py
View file

@ -12,5 +12,4 @@
# 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.
from .density import DensityField # noqa

125
csiborgtools/field/density.py
View file

@ -15,12 +15,12 @@
"""
Density field and cross-correlation calculations.
"""
from warnings import warn
from tqdm import trange
import numpy
import MAS_library as MASL
import Pk_library as PKL
import smoothing_library as SL
from warnings import warn
from tqdm import trange
from ..units import (BoxUnits, radec_to_cartesian)
@ -55,9 +55,12 @@ class DensityField:
def __init__(self, particles, boxsize, box, MAS="CIC"):
self.particles = particles
assert boxsize > 0
self._boxsize = boxsize
assert isinstance(box, BoxUnits)
self.box = box
self.boxsize = boxsize
self.MAS = MAS
assert MAS in ["NGP", "CIC", "TSC", "PCS"]
self._MAS = MAS
@property
def particles(self):
@ -89,15 +92,10 @@ class DensityField:
"""
return self._boxsize
@boxsize.setter
def boxsize(self, boxsize):
"""Set `self.boxsize`."""
self._boxsize = boxsize
@property
def box(self):
"""
The simulation box information and transformations.
Simulation box information and transformations.
Returns
-------
@ -105,17 +103,10 @@ class DensityField:
"""
return self._box
@box.setter
def box(self, box):
"""Set the simulation box."""
if not isinstance(box, BoxUnits):
raise TypeError("`box` must be `BoxUnits` instance.")
self._box = box
@property
def MAS(self):
"""
The mass-assignment scheme.
Mass-assignment scheme.
Returns
-------
@ -123,15 +114,6 @@ class DensityField:
"""
return self._MAS
@MAS.setter
def MAS(self, MAS):
"""Sets `self.MAS`."""
opts = ["NGP", "CIC", "TSC", "PCS"]
if MAS not in opts:
raise ValueError("Invalid MAS `{}`. Options are: `{}`."
.format(MAS, opts))
self._MAS = MAS
@staticmethod
def _force_f32(x, name):
if x.dtype != numpy.float32:
@ -147,17 +129,16 @@ class DensityField:
Parameters
----------
grid : int
The grid size.
Grid size.
smooth_scale : float, optional
Scale to smoothen the density field, in units matching
`self.boxsize`. By default `None`, i.e. no smoothing is applied.
verbose : float, optional
A verbosity flag. By default `True`.
`self.boxsize`. By default no smoothing is applied.
verbose : bool
Verbosity flag.
Returns
-------
rho : 3-dimensional array of shape `(grid, grid, grid)`.
Density field.
References
----------
@ -185,17 +166,16 @@ class DensityField:
Parameters
----------
grid : int
The grid size.
Grid size.
smooth_scale : float, optional
Scale to smoothen the density field, in units matching
`self.boxsize`. By default `None`, i.e. no smoothing is applied.
verbose : float, optional
A verbosity flag. By default `True`.
`self.boxsize`. By default no smoothing is applied.
verbose : bool
Verbosity flag.
Returns
-------
overdensity : 3-dimensional array of shape `(grid, grid, grid)`.
Overdensity field.
"""
# Get the overdensity
delta = self.density_field(grid, smooth_scale, verbose)
@ -210,17 +190,16 @@ class DensityField:
Parameters
----------
grid : int
The grid size.
Grid size.
smooth_scale : float, optional
Scale to smoothen the original density field, in units matching
`self.boxsize`. By default `None`, i.e. no smoothing is applied.
verbose : float, optional
A verbosity flag. By default `True`.
Scale to smoothen the density field, in units matching
`self.boxsize`. By default no smoothing is applied.
verbose : bool
Verbosity flag.
Returns
-------
potential : 3-dimensional array of shape `(grid, grid, grid)`.
Potential field.
"""
delta = self.overdensity_field(grid, smooth_scale, verbose)
if verbose:
@ -236,12 +215,12 @@ class DensityField:
Parameters
----------
grid : int
The grid size.
Grid size.
smooth_scale : float, optional
Scale to smoothen the original density field, in units matching
`self.boxsize`. By default `None`, i.e. no smoothing is applied.
verbose : float, optional
A verbosity flag. By default `True`.
Scale to smoothen the density field, in units matching
`self.boxsize`. By default no smoothing is applied.
verbose : bool
Verbosity flag.
Returns
-------
@ -260,12 +239,12 @@ class DensityField:
Parameters
----------
grid : int
The grid size.
Grid size.
smooth_scale : float, optional
Scale to smoothen the original density field, in units matching
`self.boxsize`. By default `None`, i.e. no smoothing is applied.
verbose : float, optional
A verbosity flag. By default `True`.
Scale to smoothen the density field, in units matching
`self.boxsize`. By default no smoothing is applied.
verbose : bool, optional
A verbosity flag.
Returns
-------
@ -284,24 +263,23 @@ class DensityField:
Parameters
----------
grid : int
The grid size.
Grid size.
smooth_scale : float, optional
Scale to smoothen the original density field, in units matching
`self.boxsize`. By default `None`, i.e. no smoothing is applied.
verbose : float, optional
A verbosity flag. By default `True`.
Scale to smoothen the density field, in units matching
`self.boxsize`. By default no smoothing is applied.
verbose : bool, optional
Verbosity flag.
Returns
-------
pk : py:class`Pk_library.Pk`
Power spectrum object.
"""
delta = self.overdensity_field(grid, smooth_scale, verbose)
return PKL.Pk(
delta, self.boxsize, axis=1, MAS=self.MAS, threads=1,
verbose=verbose)
def smooth_field(self, field, scale, threads=1):
def smooth_field(self, field, smooth_scale, threads=1):
"""
Smooth a field with a Gaussian filter.
@ -309,9 +287,9 @@ class DensityField:
----------
field : 3-dimensional array of shape `(grid, grid, grid)`
The field to be smoothed.
scale : float
The smoothing scale of the Gaussian filter. Units must match that
of `self.boxsize`.
smooth_scale : float, optional
Scale to smoothen the density field, in units matching
`self.boxsize`. By default no smoothing is applied.
threads : int, optional
Number of threads. By default 1.
@ -322,17 +300,17 @@ class DensityField:
Filter = "Gaussian"
grid = field.shape[0]
# FFT of the filter
W_k = SL.FT_filter(self.boxsize, scale, grid, Filter, threads)
W_k = SL.FT_filter(self.boxsize, smooth_scale, grid, Filter, threads)
return SL.field_smoothing(field, W_k, threads)
def evaluate_field(self, *field, pos):
"""
Evaluate the field at Cartesian coordinates.
Evaluate the field at Cartesian coordinates using CIC interpolation.
Parameters
----------
field : (list of) 3-dimensional array of shape `(grid, grid, grid)`
The density field that is to be interpolated.
Density field that is to be interpolated.
pos : 2-dimensional array of shape `(n_samples, 3)`
Positions to evaluate the density field. The coordinates span range
of [0, boxsize].
@ -340,7 +318,6 @@ class DensityField:
Returns
-------
interp_field : (list of) 1-dimensional array of shape `(n_samples,).
Interpolated fields at `pos`.
"""
self._force_f32(pos, "pos")
@ -353,12 +330,13 @@ class DensityField:
def evaluate_sky(self, *field, pos, isdeg=True):
"""
Evaluate the field at given distance, right ascension and declination.
Assumes that the observed is in the centre of the box.
Assumes that the observed is in the centre of the box and uses CIC
interpolation.
Parameters
----------
field : (list of) 3-dimensional array of shape `(grid, grid, grid)`
The density field that is to be interpolated. Assumed to be defined
Density field that is to be interpolated. Assumed to be defined
on a Cartesian grid.
pos : 2-dimensional array of shape `(n_samples, 3)`
Spherical coordinates to evaluate the field. Should be distance,
@ -369,7 +347,6 @@ class DensityField:
Returns
-------
interp_field : (list of) 1-dimensional array of shape `(n_samples,).
Interpolated fields at `pos`.
"""
self._force_f32(pos, "pos")
X = numpy.vstack(
@ -387,12 +364,11 @@ class DensityField:
Parameters
----------
gx, gy, gz : 1-dimensional arrays of shape `(n_samples,)`
Gravitational field components.
Gravitational field Cartesian components.
Returns
-------
g : 1-dimensional array of shape `(n_samples,)`
Gravitational field norm.
"""
return numpy.sqrt(gx * gx + gy * gy + gz * gz)
@ -410,7 +386,6 @@ class DensityField:
Returns
-------
eigvals : 2-dimensional array of shape `(n_samples, 3)`
Eigenvalues of each sample.
"""
n_samples = T00.size
# Fill array of shape `(n_samples, 3, 3)` to calculate eigvals
@ -446,14 +421,14 @@ class DensityField:
Directions to evaluate the field. Assumes `dec` is in [-90, 90]
degrees (or equivalently in radians).
field : 3-dimensional array of shape `(grid, grid, grid)`
The density field that is to be interpolated. Assumed to be defined
Density field that is to be interpolated. Assumed to be defined
on a Cartesian grid `[0, self.boxsize]^3`.
dist_marg : 1-dimensional array
Radial distances to evaluate the field.
isdeg : bool, optional
Whether `ra` and `dec` are in degres. By default `True`.
verbose : bool, optional
Verbosity flag. By default `True`.
Verbosity flag.
Returns
-------

1
csiborgtools/fits/__init__.py
View file

@ -12,7 +12,6 @@
# 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.
from .haloprofile import (NFWProfile, NFWPosterior) # noqa
from .halo import (distribute_halos, clump_with_particles, # noqa
dump_split_particles, load_split_particles, # noqa

253
csiborgtools/fits/halo.py
View file

@ -15,12 +15,10 @@
"""
Tools for splitting the particles and a clump object.
"""
import numpy
from os import remove
from warnings import warn
from os.path import join
import numpy
from tqdm import trange
from ..read import ParticleReader
@ -34,11 +32,11 @@ def clump_with_particles(particle_clumps, clumps):
particle_clumps : 1-dimensional array
Array of particles' clump IDs.
clumps : structured array
The clumps array.
Clumps array.
Returns
-------
with_particles : 1-dimensional array
with_particles : 1-dimensional boolean array
Array of whether a clump has any particles.
"""
return numpy.isin(clumps["index"], particle_clumps)
@ -54,13 +52,12 @@ def distribute_halos(n_splits, clumps):
n_splits : int
Number of splits.
clumps : structured array
The clumps array.
Clumps array.
Returns
-------
splits : 2-dimensional array
Array of starting and ending indices of each CPU of shape
`(njobs, 2)`.
splits : 2-dimensional array of shape `(njobs, 2)`
Array of starting and ending indices of each CPU.
"""
# Make sure these are unique IDs
indxs = clumps["index"]
@ -75,7 +72,7 @@ def distribute_halos(n_splits, clumps):
def dump_split_particles(particles, particle_clumps, clumps, n_splits,
paths, verbose=True):
nsnap, nsim, paths, verbose=True):
"""
Save the data needed for each split so that a process does not have to load
everything.
@ -90,6 +87,10 @@ def dump_split_particles(particles, particle_clumps, clumps, n_splits,
The clumps array.
n_splits : int
Number of times to split the clumps.
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
paths : py:class`csiborgtools.read.CSiBORGPaths`
CSiBORG paths-handling object with set `n_sim` and `n_snap`.
verbose : bool, optional
@ -112,9 +113,8 @@ def dump_split_particles(particles, particle_clumps, clumps, n_splits,
splits = distribute_halos(n_splits, clumps)
fname = join(paths.temp_dumpdir, "out_{}_snap_{}_{}.npz")
iters = trange(n_splits) if verbose else range(n_splits)
tot = 0
for n in iters:
for n in trange(n_splits) if verbose else range(n_splits):
# Lower and upper array index of the clumps array
i, j = splits[n, :]
# Clump indices in this split
@ -130,7 +130,7 @@ def dump_split_particles(particles, particle_clumps, clumps, n_splits,
"with no particles.".format(n, indxs.size, npart_unique))
# Dump it!
tot += mask.sum()
fout = fname.format(paths.n_sim, paths.n_snap, n)
fout = fname.format(nsim, nsnap, n)
numpy.savez(fout, particles[mask], particle_clumps[mask], clumps[i:j])
# There are particles whose clump ID is > 1 and have no counterpart in the
@ -166,7 +166,7 @@ def split_jobs(njobs, ncpu):
return jobs
def load_split_particles(n_split, paths, remove_split=False):
def load_split_particles(nsplit, nsnap, nsim, paths, remove_split=False):
"""
Load particles of a split saved by `dump_split_particles`.
@ -174,6 +174,10 @@ def load_split_particles(n_split, paths, remove_split=False):
--------
n_split : int
Split index.
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
paths : py:class`csiborgtools.read.CSiBORGPaths`
CSiBORG paths-handling object with set `n_sim` and `n_snap`.
remove_split : bool, optional
@ -188,9 +192,8 @@ def load_split_particles(n_split, paths, remove_split=False):
clumps : 1-dimensional array
Clumps belonging to this split.
"""
fname = join(
paths.temp_dumpdir, "out_{}_snap_{}_{}.npz".format(
paths.n_sim, paths.n_snap, n_split))
fname = join(paths.temp_dumpdir,
"out_{}_snap_{}_{}.npz".format(nsim, nsnap, nsplit))
file = numpy.load(fname)
particles, clump_indxs, clumps = (file[f] for f in file.files)
if remove_split:
@ -228,6 +231,11 @@ def pick_single_clump(n, particles, particle_clumps, clumps):
return particles[mask], clumps[n]
###############################################################################
# Clump object #
###############################################################################
class Clump:
r"""
A clump (halo) object to handle the particles and their clump's data.
@ -264,28 +272,31 @@ class Clump:
G : float, optional
The gravitational constant :math:`G` in box units. By default not set.
"""
_r = None
_rmin = None
_rmax = None
_pos = None
_clump_pos = None
_clump_mass = None
_vel = None
_Npart = None
_index = None
_rhoc = None
_G = None
def __init__(self, x, y, z, m, x0, y0, z0, clump_mass=None,
vx=None, vy=None, vz=None, index=None, rhoc=None, G=None):
self.pos = (x, y, z, x0, y0, z0)
self.clump_pos = (x0, y0, z0)
self.clump_mass = clump_mass
self.vel = (vx, vy, vz)
self.m = m
self.index = index
self.rhoc = rhoc
self.G = G
self._pos = numpy.vstack([x - x0, y - y0, z - z0]).T
self._clump_pos = numpy.asarray((x0, y0, z0))
assert clump_mass is None or isinstance(clump_mass, float)
self._clump_mass = clump_mass
if all(v is not None for v in (vx, vy, vz)):
self._vel = numpy.vstack([vx, vy, vz]).T
assert self._vel.shape == self.pos.shape
assert m.ndim == 1 and m.size == self.Npart
self._m = m
assert index is None or (isinstance(index, (int, numpy.int64)) and index >= 0) # noqa
self._index = index
assert rhoc is None or rhoc > 0
self._rhoc = rhoc
assert G is None or G > 0
self._G = G
@property
def pos(self):
@ -298,15 +309,16 @@ class Clump:
"""
return self._pos
@pos.setter
def pos(self, X):
"""Sets `pos` and calculates radial distance."""
x, y, z, x0, y0, z0 = X
self._pos = numpy.vstack([x - x0, y - y0, z - z0]).T
self._r = numpy.sum(self.pos**2, axis=1)**0.5
self._rmin = numpy.min(self._r)
self._rmax = numpy.max(self._r)
self._Npart = self._r.size
@property
def Npart(self):
"""
Number of particles associated with this clump.
Returns
-------
Npart : int
"""
return self.pos.shape[0]
@property
def r(self):
@ -317,7 +329,7 @@ class Clump:
-------
r : 1-dimensional array of shape `(n_particles, )`.
"""
return self._r
return numpy.sum(self.pos**2, axis=1)**0.5
@property
def rmin(self):
@ -328,7 +340,7 @@ class Clump:
-------
rmin : float
"""
return self._rmin
return numpy.min(self.r)
@property
def rmax(self):
@ -339,23 +351,12 @@ class Clump:
-------
rmin : float
"""
return self._rmax
@property
def Npart(self):
"""
Number of particles associated with this clump.
Returns
-------
Npart : int
"""
return self._Npart
return numpy.max(self.r)
@property
def clump_pos(self):
"""
Cartesian clump coordinates.
Cartesian position components of the clump.
Returns
-------
@ -363,14 +364,6 @@ class Clump:
"""
return self._clump_pos
@clump_pos.setter
def clump_pos(self, pos):
"""Sets `clump_pos`. Makes sure it is the correct shape."""
pos = numpy.asarray(pos)
if pos.shape != (3,):
raise TypeError("Invalid clump position `{}`".format(pos.shape))
self._clump_pos = pos
@property
def clump_mass(self):
"""
@ -384,17 +377,10 @@ class Clump:
raise ValueError("Clump mass `clump_mass` has not been set.")
return self._clump_mass
@clump_mass.setter
def clump_mass(self, mass):
"""Sets `clump_mass`, making sure it is a float."""
if mass is not None and not isinstance(mass, float):
raise ValueError("`clump_mass` must be a float.")
self._clump_mass = mass
@property
def vel(self):
"""
Cartesian particle velocities. Throws an error if they are not set.
Cartesian velocity components of the clump.
Returns
-------
@ -404,16 +390,6 @@ class Clump:
raise ValueError("Velocities `vel` have not been set.")
return self._vel
@vel.setter
def vel(self, V):
"""Sets the particle velocities, making sure the shape is OK."""
if any(v is None for v in V):
return
vx, vy, vz = V
self._vel = numpy.vstack([vx, vy, vz]).T
if self.pos.shape != self.vel.shape:
raise ValueError("Different `pos` and `vel` arrays!")
@property
def m(self):
"""
@ -425,18 +401,11 @@ class Clump:
"""
return self._m
@m.setter
def m(self, m):
"""Sets particle masses `m`, ensuring it is the right size."""
if not isinstance(m, numpy.ndarray) and m.size != self.r.size:
raise TypeError("`r` and `m` must be equal size 1-dim arrays.")
self._m = m
@property
def center_mass(self):
"""
Clump center of mass. Note that this is already in a frame centered at
the clump's potential minimum.
Cartesian position components of the clump centre of mass. Note that
this is already in a frame centered at the clump's potential minimum.
Returns
-------
@ -444,10 +413,43 @@ class Clump:
"""
return numpy.average(self.pos, axis=0, weights=self.m)
@property
def angular_momentum(self):
"""
Clump angular momentum in the box coordinates.
Returns
-------
J : 1-dimensional array or shape `(3, )`
"""
J = numpy.cross(self.pos - self.center_mass, self.vel)
return numpy.einsum("i,ij->j", self.m, J)
@property
def lambda200c(self):
r"""
Clump Bullock spin, see Eq. 5 in [1], in a radius of
:math:`R_{\rm 200c}`.
References
----------
[1] A Universal Angular Momentum Profile for Galactic Halos; 2001;
Bullock, J. S.; Dekel, A.; Kolatt, T. S.; Kravtsov, A. V.;
Klypin, A. A.; Porciani, C.; Primack, J. R.
Returns
-------
lambda200c : float
"""
J = self.angular_momentum
R, M = self.spherical_overdensity_mass(200)
V = numpy.sqrt(self.G * M / R)
return numpy.linalg.norm(J) / (numpy.sqrt(2) * M * V * R)
@property
def index(self):
"""
The halo finder clump index.
Halo finder clump index.
Returns
-------
@ -457,17 +459,10 @@ class Clump:
raise ValueError("Halo index `hindex` has not been set.")
return self._index
@index.setter
def index(self, n):
"""Sets the halo index, making sure it is an integer."""
if n is not None and not (isinstance(n, (int, numpy.int64)) and n > 0):
raise ValueError("Halo index `index` must be an integer > 0.")
self._index = n
@property
def rhoc(self):
r"""
The critical density :math:`\rho_c` at this snapshot in box units.
Critical density :math:`\rho_c` at this snapshot in box units.
Returns
-------
@ -477,17 +472,10 @@ class Clump:
raise ValueError("The critical density `rhoc` has not been set.")
return self._rhoc
@rhoc.setter
def rhoc(self, rhoc):
"""Sets the critical density. Makes sure it is > 0."""
if rhoc is not None and not rhoc > 0:
raise ValueError("Critical density `rho_c` must be > 0.")
self._rhoc = rhoc
@property
def G(self):
r"""
The gravitational constant :math:`G` in box units.
Gravitational constant :math:`G` in box units.
Returns
-------
@ -497,13 +485,6 @@ class Clump:
raise ValueError("The grav. constant `G` has not been set.")
return self._G
@G.setter
def G(self, G):
"""Sets the gravitational constant. Makes sure it is > 0."""
if G is not None and not G > 0:
raise ValueError("Gravitational constant `G` must be > 0.")
self._G = G
@property
def total_particle_mass(self):
"""
@ -528,8 +509,7 @@ class Clump:
def enclosed_spherical_mass(self, rmax, rmin=0):
"""
The enclosed spherical mass between two radii. All quantities remain
in the box units.
Enclosed spherical mass between two radii in box units.
Parameters
----------
@ -547,14 +527,14 @@ class Clump:
def enclosed_spherical_volume(self, rmax, rmin=0):
"""
Enclosed volume within two radii.
Enclosed spherical volume within two radii in box units.
Parameters
----------
rmax : float
The maximum radial distance.
Maximum radial distance.
rmin : float, optional
The minimum radial distance. By default 0.
Minimum radial distance. By default 0.
Returns
-------
@ -581,15 +561,15 @@ class Clump:
The :math:`\delta_{\rm x}` parameters where :math:`\mathrm{x}` is
the overdensity multiple.
n_particles_min : int
The minimum number of enclosed particles for a radius to be
Minimum number of enclosed particles for a radius to be
considered trustworthy.
Returns
-------
rx : float
The radius where the enclosed density reaches required value.
Radius where the enclosed density reaches required value.
mx : float
The corresponding spherical enclosed mass.
Corresponding spherical enclosed mass.
"""
# If single `delta` turn to list
delta = [delta] if isinstance(delta, (float, int)) else delta
@ -636,39 +616,6 @@ class Clump:
return rfound, mfound
@property
def angular_momentum(self):
"""
The clump angular momentum in the box coordinates.
Returns
-------
J : 1-dimensional array or shape `(3, )`
"""
J = numpy.cross(self.pos - self.center_mass, self.vel)
return numpy.einsum("i,ij->j", self.m, J)
@property
def lambda200c(self):
r"""
The clump Bullock spin, see Eq. 5 in [1], in a radius of
:math:`R_{\rm 200c}`.
References
----------
[1] A Universal Angular Momentum Profile for Galactic Halos; 2001;
Bullock, J. S.; Dekel, A.; Kolatt, T. S.; Kravtsov, A. V.;
Klypin, A. A.; Porciani, C.; Primack, J. R.
Returns
-------
lambda200c : float
"""
J = self.angular_momentum
R, M = self.spherical_overdensity_mass(200)
V = numpy.sqrt(self.G * M / R)
return numpy.linalg.norm(J) / (numpy.sqrt(2) * M * V * R)
@classmethod
def from_arrays(cls, particles, clump, rhoc=None, G=None):
r"""

80
csiborgtools/fits/haloprofile.py
View file

@ -15,37 +15,14 @@
"""
Halo profiles functions and posteriors.
"""
import numpy
from jax import numpy as jnumpy
from jax import grad
import numpy
from scipy.optimize import minimize_scalar
from scipy.stats import uniform
from .halo import Clump
def jax_switch(f1, f2, use_jax):
"""
Return function `f1` if `use_jax` else return `f2`.
Parameters
----------
f1 : function
The JAX function.
f2 : function
The non-JAX function.
use_jax : bool
Whether to use the JAX function.
Returns
-------
chosen_func : function
The chosen function.
"""
return f1 if use_jax else f2
class NFWProfile:
r"""
The Navarro-Frenk-White (NFW) density profile defined as
@ -62,10 +39,6 @@ class NFWProfile:
rho0 : float
NFW density parameter :math:`\rho_0`.
"""
def __init__(self):
pass
@staticmethod
def profile(r, Rs, rho0):
r"""
@ -109,7 +82,7 @@ class NFWProfile:
logdensity : float or 1-dimensional array
Logarithmic density of the NFW profile at :math:`r`.
"""
log = jax_switch(jnumpy.log, numpy.log, use_jax)
log = jnumpy.log if use_jax else numpy.log
x = r / Rs
return log(rho0) - log(x) - 2 * log(1 + x)
@ -134,8 +107,9 @@ class NFWProfile:
M : float or 1-dimensional array
The enclosed mass.
"""
log = jnumpy.log if use_jax else numpy.log
x = r / Rs
out = jax_switch(jnumpy.log, numpy.log, use_jax)(1 + x) - x / (1 + x)
out = log(1 + x) - x / (1 + x)
return 4 * numpy.pi * rho0 * Rs**3 * out
def bounded_enclosed_mass(self, rmin, rmax, Rs, rho0, use_jax=False):
@ -145,9 +119,9 @@ class NFWProfile:
Parameters
----------
rmin : float
The minimum radius.
Minimum radius.
rmax : float
The maximum radius.
Maximum radius.
Rs : float
Scale radius :math:`R_s`.
rho0 : float
@ -158,7 +132,7 @@ class NFWProfile:
Returns
-------
M : float
The enclosed mass within the radial range.
Enclosed mass within the radial range.
"""
return (self.enclosed_mass(rmax, Rs, rho0, use_jax)
- self.enclosed_mass(rmin, Rs, rho0, use_jax))
@ -182,9 +156,9 @@ class NFWProfile:
Rs : float
Scale radius :math:`R_s`.
rmin : float
The minimum radius.
Minimum radius.
rmax : float
The maximum radius.
Maximum radius.
Returns
-------
@ -202,9 +176,9 @@ class NFWProfile:
Parameters
----------
rmin : float
The minimum radius.
Minimum radius.
rmax : float
The maximum radius.
Maximum radius.
Rs : float
Scale radius :math:`R_s`.
N : int, optional
@ -271,7 +245,6 @@ class NFWPosterior(NFWProfile):
Returns
-------
clump : `Clump`
The clump object.
"""
return self._clump
@ -284,7 +257,6 @@ class NFWPosterior(NFWProfile):
Returns
-------
r : 1-dimensional array
Array of shape `(n_particles, )`.
"""
return self._r
@ -297,7 +269,6 @@ class NFWPosterior(NFWProfile):
Returns
-------
Npart : int
Number of particles.
"""
return self._Npart
@ -310,7 +281,6 @@ class NFWPosterior(NFWProfile):
Returns
-------
r : 1-dimensional array
Array of shape `(n_particles, )`.
"""
return self._m
@ -322,7 +292,6 @@ class NFWPosterior(NFWProfile):
Returns
-------
rmin : float
The minimum distance.
"""
return self._rmin
@ -335,7 +304,6 @@ class NFWPosterior(NFWProfile):
Returns
-------
rmax : float
The R200c radius.
"""
return self._rmax
@ -384,7 +352,6 @@ class NFWPosterior(NFWProfile):
Returns
-------
rho0: float
The NFW density parameter.
"""
Mtot = numpy.exp(self._logMtot)
Mnfw_norm = self.bounded_enclosed_mass(self.rmin, self.rmax, Rs, 1)
@ -401,8 +368,7 @@ class NFWPosterior(NFWProfile):
Returns
-------
ll : float
The logarithmic prior.
lp : float
"""
if not self._logrmin < logRs < self._logrmax:
return - numpy.infty
@ -422,16 +388,14 @@ class NFWPosterior(NFWProfile):
Returns
-------
ll : float
The logarithmic likelihood.
"""
Rs = 10**logRs
log = jnumpy.log if use_jax else numpy.log
# Expected enclosed mass from a NFW
Mnfw = self.bounded_enclosed_mass(self.rmin, self.rmax,
Rs, 1, use_jax)
ll = jax_switch(jnumpy.sum, numpy.sum, use_jax)(self.logprofile(
self.r, Rs, 1, use_jax))
ll += self._ll0
fsum = jnumpy.sum if use_jax else numpy.sum
ll = fsum(self.logprofile(self.r, Rs, 1, use_jax)) + self._ll0
return ll - self.Npart * log(Mnfw)
@property
@ -444,7 +408,6 @@ class NFWPosterior(NFWProfile):
Returns
-------
initlogRs : float
The initial guess of :math:`\log R_{\rm s}`.
"""
bins = numpy.linspace(self.rmin, self.rmax,
self._binsguess)
@ -465,7 +428,6 @@ class NFWPosterior(NFWProfile):
Returns
-------
lpost : float
The logarithmic posterior.
"""
lp = self.logprior(logRs)
if not numpy.isfinite(lp):
@ -475,7 +437,7 @@ class NFWPosterior(NFWProfile):
def uncertainty_at_maxpost(self, logRs_max):
r"""
Calculate Gaussian approximation of the uncertainty at `logRs_max`, the
maximum a-posteriori esitmate. This is the square root of the negative
maximum a-posteriori estimate. This is the square root of the negative
inverse 2nd derivate of the logarithimic posterior with respect to the
logarithm of the scale factor. This is only valid `logRs_max` is the
maximum of the posterior!
@ -492,7 +454,6 @@ class NFWPosterior(NFWProfile):
Returns
-------
uncertainty : float
The uncertainty calculated as outlined above.
"""
def f(x):
return self(x, use_jax=True)
@ -506,9 +467,10 @@ class NFWPosterior(NFWProfile):
def maxpost_logRs(self, calc_err=False, eps=1e-4):
r"""
Maximum a-posterio estimate of the scale radius :math:`\log R_{\rm s}`.
Returns the scale radius if the fit converged, otherwise `numpy.nan`.
Checks whether :math:`log r_{\rm max} / R_{\rm s} > \epsilon`, where
Maximum a-posteriori estimate of the scale radius
:math:`\log R_{\rm s}`. Returns the scale radius if the fit converged,
otherwise `numpy.nan`. Checks whether
:math:`log r_{\rm max} / R_{\rm s} > \epsilon`, where
to ensure that the scale radius is not too close to the boundary which
occurs if the fit fails.
@ -521,9 +483,9 @@ class NFWPosterior(NFWProfile):
Returns
-------
logRs: float
The log scale radius.
Log scale radius.
uncertainty : float
The uncertainty on the scale radius. Calculated following
Uncertainty on the scale radius. Calculated following
`self.uncertainty_at_maxpost`.
"""
# Loss function to optimize

2
csiborgtools/match/__init__.py
View file

@ -13,7 +13,7 @@
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
from .match import (brute_spatial_separation, RealisationsMatcher, cosine_similarity, # noqa
from .match import (RealisationsMatcher, cosine_similarity, # noqa
ParticleOverlap, get_clumplims, fill_delta, fill_delta_indxs, # noqa
calculate_overlap, calculate_overlap_indxs, # noqa
dist_centmass, dist_percentile) # noqa

1
csiborgtools/match/correlation.py
View file

@ -12,7 +12,6 @@
# 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.
import numpy
from Corrfunc.mocks import DDtheta_mocks
from Corrfunc.utils import convert_3d_counts_to_cf

1
csiborgtools/match/knn.py
View file

@ -258,7 +258,6 @@ class kNN_CDF:
range=(rmin, rmax))
cdf1[k, :] += _counts
joint_cdf = numpy.cumsum(joint_cdf, axis=-1)
cdf0 = numpy.cumsum(cdf0, axis=-1)
cdf1 = numpy.cumsum(cdf1, axis=-1)

156
csiborgtools/match/match.py
View file

@ -12,64 +12,21 @@
# 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.
"""
Support for matching halos between CSiBORG IC realisations.
"""
from gc import collect
import numpy
from scipy.ndimage import gaussian_filter
from tqdm import (tqdm, trange)
from astropy.coordinates import SkyCoord
from numba import jit
from gc import collect
from tqdm import (tqdm, trange)
from ..read import concatenate_clumps
from ..utils import now
def brute_spatial_separation(c1, c2, angular=False, N=None, verbose=False):
"""
Calculate for each point in `c1` the `N` closest points in `c2`.
Parameters
----------
c1 : `astropy.coordinates.SkyCoord`
Coordinates of the first set of points.
c2 : `astropy.coordinates.SkyCoord`
Coordinates of the second set of points.
angular : bool, optional
Whether to calculate angular separation or 3D separation. By default
`False` and 3D separation is calculated.
N : int, optional
Number of closest points in `c2` to each object in `c1` to return.
verbose : bool, optional
Verbosity flag. By default `False`.
Returns
-------
sep : 1-dimensional array
Separation of each object in `c1` to `N` closest objects in `c2`. The
array shape is `(c1.size, N)`. Separation is in units of `c1`.
indxs : 1-dimensional array
Indexes of the closest objects in `c2` for each object in `c1`. The
array shape is `(c1.size, N)`.
"""
if not (isinstance(c1, SkyCoord) and isinstance(c2, SkyCoord)):
raise TypeError("`c1` & `c2` must be `astropy.coordinates.SkyCoord`.")
N1 = c1.size
N2 = c2.size if N is None else N
# Pre-allocate arrays
sep = numpy.full((N1, N2), numpy.nan)
indxs = numpy.full((N1, N2), numpy.nan, dtype=int)
iters = tqdm(range(N1)) if verbose else range(N1)
for i in iters:
if angular:
dist = c1[i].separation(c2).value
else:
dist = c1[i].separation_3d(c2).value
# Sort the distances
sort = numpy.argsort(dist)[:N2]
indxs[i, :] = sort
sep[i, :] = dist[sort]
return sep, indxs
###############################################################################
# Realisations matcher for calculating overlaps #
###############################################################################
class RealisationsMatcher:
@ -101,9 +58,12 @@ class RealisationsMatcher:
def __init__(self, nmult=1., dlogmass=2., mass_kind="totpartmass",
overlapper_kwargs={}):
self.nmult = nmult
self.dlogmass = dlogmass
self.mass_kind = mass_kind
assert nmult > 0
assert dlogmass > 0
assert isinstance(mass_kind, str)
self._nmult = nmult
self._dlogmass = dlogmass
self._mass_kind = mass_kind
self._overlapper = ParticleOverlap(**overlapper_kwargs)
@property
@ -118,12 +78,6 @@ class RealisationsMatcher:
"""
return self._nmult
@nmult.setter
def nmult(self, nmult):
"""Set `nmult`."""
assert nmult > 0
self._nmult = nmult
@property
def dlogmass(self):
"""
@ -136,16 +90,10 @@ class RealisationsMatcher:
"""
return self._dlogmass
@dlogmass.setter
def dlogmass(self, dlogmass):
"""Set `dlogmass`."""
assert dlogmass > 0
self._dlogmass = dlogmass
@property
def mass_kind(self):
"""
The mass kind whose similarity is to be checked.
Mass kind whose similarity is to be checked.
Returns
-------
@ -153,12 +101,6 @@ class RealisationsMatcher:
"""
return self._mass_kind
@mass_kind.setter
def mass_kind(self, mass_kind):
"""Set `mass_kind`."""
assert isinstance(mass_kind, str)
self._mass_kind = mass_kind
@property
def overlapper(self):
"""
@ -213,7 +155,7 @@ class RealisationsMatcher:
# Query the KNN
verbose and print("{}: querying the KNN.".format(now()), flush=True)
match_indxs = radius_neighbours(
catx.knn(select_initial=True), cat0.positions0,
catx.knn(select_initial=True), cat0.positions(in_initial=True),
radiusX=cat0["lagpatch"], radiusKNN=catx["lagpatch"],
nmult=self.nmult, enforce_in32=True, verbose=verbose)
@ -381,10 +323,6 @@ class ParticleOverlap:
the nearest grid position particle assignment scheme, with optional
Gaussian smoothing.
"""
_inv_clength = None
_clength = None
_nshift = None
def __init__(self):
# Inverse cell length in box units. By default :math:`2^11`, which
# matches the initial RAMSES grid resolution.
@ -679,6 +617,11 @@ class ParticleOverlap:
nonzero, mass1, mass2)
###############################################################################
# Halo matching supplementary functions #
###############################################################################
@jit(nopython=True)
def fill_delta(delta, xcell, ycell, zcell, xmin, ymin, zmin, weights):
"""
@ -757,7 +700,7 @@ def get_clumplims(clumps, ncells, nshift=None):
Returns
-------
mins, maxs : 2-dimensional arrays of shape `(n_samples, 3)`
The minimum and maximum along each axis.
Minimum and maximum along each axis.
"""
dtype = clumps[0][0]['x'].dtype # dtype of the first clump's 'x'
# Check that for real positions we cannot apply nshift
@ -983,3 +926,58 @@ def radius_neighbours(knn, X, radiusX, radiusKNN, nmult=1., enforce_in32=False,
indxs[i] = indxs[i].astype(numpy.int32)
return numpy.asarray(indxs, dtype=object)
###############################################################################
# Sky mathing #
###############################################################################
# def brute_spatial_separation(c1, c2, angular=False, N=None, verbose=False):
# """
# Calculate for each point in `c1` the `N` closest points in `c2`.
# Parameters
# ----------
# c1 : `astropy.coordinates.SkyCoord`
# Coordinates of the first set of points.
# c2 : `astropy.coordinates.SkyCoord`
# Coordinates of the second set of points.
# angular : bool, optional
# Whether to calculate angular separation or 3D separation. By default
# `False` and 3D separation is calculated.
# N : int, optional
# Number of closest points in `c2` to each object in `c1` to return.
# verbose : bool, optional
# Verbosity flag. By default `False`.
# Returns
# -------
# sep : 1-dimensional array
# Separation of each object in `c1` to `N` closest objects in `c2`. The
# array shape is `(c1.size, N)`. Separation is in units of `c1`.
# indxs : 1-dimensional array
# Indexes of the closest objects in `c2` for each object in `c1`. The
# array shape is `(c1.size, N)`.
# """
# if not (isinstance(c1, SkyCoord) and isinstance(c2, SkyCoord)):
# raise TypeError(
# "`c1` & `c2` must be `astropy.coordinates.SkyCoord`.")
# N1 = c1.size
# N2 = c2.size if N is None else N
# # Pre-allocate arrays
# sep = numpy.full((N1, N2), numpy.nan)
# indxs = numpy.full((N1, N2), numpy.nan, dtype=int)
# iters = tqdm(range(N1)) if verbose else range(N1)
# for i in iters:
# if angular:
# dist = c1[i].separation(c2).value
# else:
# dist = c1[i].separation_3d(c2).value
# # Sort the distances
# sort = numpy.argsort(dist)[:N2]
# indxs[i, :] = sort
# sep[i, :] = dist[sort]
# return sep, indxs

1
csiborgtools/match/num_density.py
View file

@ -15,7 +15,6 @@
"""
Calculation of number density functions.
"""
import numpy

4
csiborgtools/read/__init__.py
View file

@ -17,6 +17,6 @@ from .readsim import (CSiBORGPaths, ParticleReader, read_mmain, read_initcm, hal
from .make_cat import (HaloCatalogue, concatenate_clumps) # noqa
from .readobs import (PlanckClusters, MCXCClusters, TwoMPPGalaxies, # noqa
TwoMPPGroups, SDSS) # noqa
from .outsim import (dump_split, combine_splits, make_ascii_powmes) # noqa
from .summaries import (PKReader, kNNCDFReader, PairOverlap, NPairsOverlap,
from .outsim import (dump_split, combine_splits) # noqa
from .summaries import (PKReader, kNNCDFReader, PairOverlap, NPairsOverlap, # noqa
binned_resample_mean) # noqa

209
csiborgtools/read/make_cat.py
View file

@ -29,8 +29,10 @@ class HaloCatalogue:
Parameters
----------
paths : py:class:`csiborgtools.read.CSiBORGPaths`
CSiBORG paths-handling object with set `n_sim` and `n_snap`.
nsim : int
IC realisation index.
paths : py:class`csiborgtools.read.CSiBORGPaths`
CSiBORG paths object.
min_mass : float, optional
The minimum :math:`M_{rm tot} / M_\odot` mass. By default no threshold.
max_dist : float, optional
@ -38,20 +40,49 @@ class HaloCatalogue:
load_init : bool, optional
Whether to load the initial snapshot information. By default False.
"""
_box = None
_paths = None
_nsim = None
_data = None
_selmask = None
def __init__(self, nsim, min_mass=None, max_dist=None, load_init=False):
# Set up paths
paths = CSiBORGPaths(n_sim=nsim)
paths.n_snap = paths.get_maximum_snapshot()
self._paths = paths
self._box = BoxUnits(paths)
def __init__(self, nsim, paths, min_mass=None, max_dist=None,
load_init=False):
assert isinstance(paths, CSiBORGPaths)
self._nsim = nsim
self._paths = paths
self._set_data(min_mass, max_dist, load_init)
@property
def nsim(self):
return self._nsim
@property
def paths(self):
return self._paths
@property
def nsnap(self):
"""
Catalogue's snapshot index corresponding to the maximum simulation
snapshot index.
Returns
-------
nsnap : int
"""
return max(self.paths.get_snapshots(self.nsim))
@property
def box(self):
"""
CSiBORG box object handling unit conversion.
Returns
-------
box : :py:class:`csiborgtools.units.BoxUnits`
"""
return BoxUnits(self.nsnap, self.nsim, self.paths)
@property
def data(self):
"""
@ -61,81 +92,19 @@ class HaloCatalogue:
-------
cat : structured array
"""
if self._data is None:
raise ValueError("`data` is not set!")
return self._data
@property
def box(self):
"""
Box object, useful for change of units.
Returns
-------
box : :py:class:`csiborgtools.units.BoxUnits`
"""
return self._box
@property
def paths(self):
"""
The paths-handling object.
Returns
-------
paths : :py:class:`csiborgtools.read.CSiBORGPaths`
"""
return self._paths
@property
def n_snap(self):
"""
The snapshot ID.
Returns
-------
n_snap : int
"""
return self.paths.n_snap
@property
def n_sim(self):
"""
The initial condition (IC) realisation ID.
Returns
-------
n_sim : int
"""
return self.paths.n_sim
def knn(self, select_initial):
"""
kNN object of all halo positions.
Parameters
----------
select_initial : bool
Whether to define the KNN on the initial or final snapshot.
Returns
-------
knn : :py:class:`sklearn.neighbors.NearestNeighbors`
"""
knn = NearestNeighbors()
return knn.fit(self.positions0 if select_initial else self.positions)
def _set_data(self, min_mass, max_dist, load_init):
"""
Loads the data, merges with mmain, does various coordinate transforms.
"""
# Load the processed data
fname = "ramses_out_{}_{}.npy".format(
str(self.n_sim).zfill(5), str(self.n_snap).zfill(5))
str(self.nsim).zfill(5), str(self.nsnap).zfill(5))
data = numpy.load(join(self.paths.dumpdir, fname))
# Load the mmain file and add it to the data
mmain = read_mmain(self.n_sim, self.paths.mmain_path)
mmain = read_mmain(self.nsim, self.paths.mmain_path)
data = self.merge_mmain_to_clumps(data, mmain)
flip_cols(data, "peak_x", "peak_z")
@ -144,7 +113,7 @@ class HaloCatalogue:
# Now also load the initial positions
if load_init:
initcm = read_initcm(self.n_sim, self.paths.initmatch_path)
initcm = read_initcm(self.nsim, self.paths.initmatch_path)
if initcm is not None:
data = self.merge_initmatch_to_clumps(data, initcm)
flip_cols(data, "x0", "z0")
@ -185,7 +154,7 @@ class HaloCatalogue:
formats.append(numpy.float32)
dtype = numpy.dtype({"names": names, "formats": formats})
# Apply cuts on distance and min500 if any
# Apply cuts on distance and total particle mass if any
data = data[data["dist"] < max_dist] if max_dist is not None else data
data = (data[data["totpartmass"] > min_mass]
if min_mass is not None else data)
@ -247,76 +216,75 @@ class HaloCatalogue:
X[:, i] = initcat[p]
return add_columns(clumps, X, ["x0", "y0", "z0", "lagpatch"])
@property
def positions(self):
def positions(self, in_initial=False):
r"""
3D positions of halos in :math:`\mathrm{cMpc}`.
Cartesian position components of halos in :math:`\mathrm{cMpc}`.
Parameters
----------
in_initial : bool, optional
Whether to define the kNN on the initial or final snapshot.
Returns
-------
X : 2-dimensional array
Array of shape `(n_halos, 3)`, where the latter axis represents
`x`, `y` and `z`.
pos : 2-dimensional array of shape `(nhalos, 3)`
"""
return numpy.vstack(
[self._data["peak_{}".format(p)] for p in ("x", "y", "z")]).T
if in_initial:
ps = ["x0", "y0", "z0"]
else:
ps = ["peak_x", "peak_y", "peak_z"]
return numpy.vstack([self[p] for p in ps]).T
@property
def positions0(self):
r"""
3D positions of halos in the initial snapshot in :math:`\mathrm{cMpc}`.
Returns
-------
X : 2-dimensional array
Array of shape `(n_halos, 3)`, where the latter axis represents
`x`, `y` and `z`.
"""
try:
return numpy.vstack(
[self._data["{}".format(p)] for p in ("x0", "y0", "z0")]).T
except KeyError:
raise RuntimeError("Initial positions are not set!")
@property
def velocities(self):
"""
Cartesian velocities of halos.
Cartesian velocity components of halos. Likely in box units.
Returns
-------
vel : 2-dimensional array
Array of shape `(n_halos, 3)`.
vel : 2-dimensional array of shape `(nhalos, 3)`
"""
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.
Cartesian angular momentum components of halos in the box coordinate
system. Likely in box units.
Returns
-------
angmom : 2-dimensional array
Array of shape `(n_halos, 3)`.
angmom : 2-dimensional array of shape `(nhalos, 3)`
"""
return numpy.vstack([self["L{}".format(p)] for p in ("x", "y", "z")]).T
def radius_neigbours(self, X, radius, select_initial=True):
def knn(self, in_initial):
"""
kNN object of all halo positions.
Parameters
----------
in_initial : bool
Whether to define the kNN on the initial or final snapshot.
Returns
-------
knn : :py:class:`sklearn.neighbors.NearestNeighbors`
"""
knn = NearestNeighbors()
return knn.fit(self.positions(in_initial))
def radius_neigbours(self, X, radius, in_initial):
r"""
Return sorted nearest neigbours within `radius` of `X` in the initial
Sorted nearest neigbours within `radius` of `X` in the initial
or final snapshot.
Parameters
----------
X : 2-dimensional array
Array of shape `(n_queries, 3)`, where the latter axis represents
`x`, `y` and `z`.
X : 2-dimensional array of shape `(n_queries, 3)`
Cartesian query position components in :math:`\mathrm{cMpc}`.
radius : float
Limiting distance of neighbours.
select_initial : bool, optional
Whether to search for neighbours in the initial or final snapshot.
By default `True`, i.e. the final snapshot.
Limiting neighbour distance.
in_initial : bool
Whether to define the kNN on the initial or final snapshot.
Returns
-------
@ -329,7 +297,7 @@ class HaloCatalogue:
"""
if not (X.ndim == 2 and X.shape[1] == 3):
raise TypeError("`X` must be an array of shape `(n_samples, 3)`.")
knn = self.knn(select_initial)
knn = self.knn(in_initial)
return knn.radius_neighbors(X, radius, sort_results=True)
@property
@ -347,6 +315,11 @@ class HaloCatalogue:
return self.data.size
###############################################################################
# Useful functions #
###############################################################################
def concatenate_clumps(clumps):
"""
Concatenate an array of clumps to a single array containing all particles.

82
csiborgtools/read/outsim.py
View file

@ -15,19 +15,13 @@
"""
I/O functions for analysing the CSiBORG realisations.
"""
import numpy
from os.path import (join, isfile)
from os.path import join
from os import remove
from tqdm import trange
from astropy.table import Table
I64 = numpy.int64
F64 = numpy.float64
def dump_split(arr, n_split, paths):
def dump_split(arr, nsplit, nsnap, nsim, paths):
"""
Dump an array from a split.
@ -35,30 +29,26 @@ def dump_split(arr, n_split, paths):
----------
arr : n-dimensional or structured array
Array to be saved.
n_split: int
The split index.
nsplit : int
Split index.
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
paths : py:class`csiborgtools.read.CSiBORGPaths`
CSiBORG paths-handling object with set `n_sim` and `n_snap`.
n_sim : int
The CSiBORG realisation index.
n_snap : int
The index of a redshift snapshot.
outdir : string
Directory where to save the temporary files.
Returns
-------
None
"""
n_sim = str(paths.n_sim).zfill(5)
n_snap = str(paths.n_snap).zfill(5)
fname = join(paths.temp_dumpdir, "ramses_out_{}_{}_{}.npy"
.format(n_sim, n_snap, n_split))
.format(str(nsim).zfill(5), str(nsnap).zfill(5), nsplit))
numpy.save(fname, arr)
def combine_splits(n_splits, part_reader, cols_add, remove_splits=False,
verbose=True):
def combine_splits(nsplits, nsnap, nsim, part_reader, cols_add,
remove_splits=False, verbose=True):
"""
Combine results of many splits saved from `dump_split`. Identifies to which
clump the clumps in the split correspond to by matching their index.
@ -67,8 +57,12 @@ def combine_splits(n_splits, part_reader, cols_add, remove_splits=False,
Paramaters
----------
n_splits : int
The total number of clump splits.
nsplits : int
Total number of clump splits.
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
part_reader : py:class`csiborgtools.read.ParticleReadear`
CSiBORG particle reader.
cols_add : list of `(str, dtype)`
@ -84,26 +78,20 @@ def combine_splits(n_splits, part_reader, cols_add, remove_splits=False,
out : structured array
Clump array with appended results from the splits.
"""
# Load clumps to see how many there are and will add to this array
n_sim = part_reader.paths.n_sim
n_snap = part_reader.paths.n_snap
clumps = part_reader.read_clumps(cols=None)
clumps = part_reader.read_clumps(nsnap, nsim, cols=None)
# Get the old + new dtypes and create an empty array
descr = clumps.dtype.descr + cols_add
out = numpy.full(clumps.size, numpy.nan, dtype=descr)
# Now put the old values into the array
for par in clumps.dtype.names:
for par in clumps.dtype.names: # Now put the old values into the array
out[par] = clumps[par]
# Filename of splits data
froot = "ramses_out_{}_{}".format(
str(n_sim).zfill(5), str(n_snap).zfill(5))
froot = "ramses_out_{}_{}".format(str(nsim).zfill(5), str(nsnap).zfill(5))
fname = join(part_reader.paths.temp_dumpdir, froot + "_{}.npy")
# Iterate over splits and add to the output array
cols_add_names = [col[0] for col in cols_add]
iters = trange(n_splits) if verbose else range(n_splits)
iters = trange(nsplits) if verbose else range(nsplits)
for n in iters:
fnamesplit = fname.format(n)
arr = numpy.load(fnamesplit)
@ -121,31 +109,3 @@ def combine_splits(n_splits, part_reader, cols_add, remove_splits=False,
remove(fnamesplit)
return out
def make_ascii_powmes(particles, fout):
"""
Write an ASCII file with appropriate formatting for POWMES.
Parameters
----------
particles : structured array
Array of particles.
fout : str
File path to store the ASCII file.
Returns
-------
None
"""
out = Table()
for p in ('x', 'y', 'z', 'M'):
out[p] = particles[p]
Npart = particles.size
# If fout exists, remove
if isfile(fout):
remove(fout)
with open(fout, 'wb') as f:
numpy.savetxt(f, out, delimiter=',', header=str(Npart), comments='')

555
csiborgtools/read/readsim.py
View file

@ -15,23 +15,15 @@
"""
Functions to read in the particle and clump files.
"""
import numpy
from scipy.io import FortranFile
from os.path import (join, isfile, isdir)
from glob import glob
from tqdm import tqdm
from warnings import warn
import numpy
from scipy.io import FortranFile
from tqdm import tqdm
from ..utils import (cols_to_structured)
F16 = numpy.float16
F32 = numpy.float32
F64 = numpy.float64
I32 = numpy.int32
I64 = numpy.int64
###############################################################################
# Paths manager #
###############################################################################
@ -43,285 +35,117 @@ class CSiBORGPaths:
Parameters
----------
n_sim : int, optional
CSiBORG IC realisation index. By default not set.
n_snap : int, optional
Snapshot index. By default not set.
srcdir : str, optional
The file path to the folder where realisations of the ICs are stored.
By default `/mnt/extraspace/hdesmond/`.
Path to the folder where CSiBORG simulations are stored.
dumpdir : str, optional
Path to where files from `run_fit_halos` are stored. By default
`/mnt/extraspace/rstiskalek/csiborg/`.
Path to the folder where files from `run_fit_halos` are stored.
mmain_path : str, optional
Path to where mmain files are stored. By default
`/mnt/zfsusers/hdesmond/Mmain`.
Path to folder where mmain files are stored.
initmatch_path : str, optional
Path to where match between the first and final snapshot is stored. By
default `/mnt/extraspace/rstiskalek/csiborg/initmatch/`.
to_new : bool, optional
Whether the paths should point to `new` files, for example
`ramses_out_8452_new`.
Path to the folder where particle ID match between the first and final
snapshot is stored.
"""
_srcdir = None
_n_sim = None
_n_snap = None
_dumpdir = None
_mmain_path = None
_initmatch_path = None
_to_new = None
def __init__(self, n_sim=None, n_snap=None, srcdir=None, dumpdir=None,
mmain_path=None, initmatch_path=None, to_new=False):
if srcdir is None:
srcdir = "/mnt/extraspace/hdesmond/"
self.srcdir = srcdir
if dumpdir is None:
dumpdir = "/mnt/extraspace/rstiskalek/csiborg/"
self.dumpdir = dumpdir
if mmain_path is None:
mmain_path = "/mnt/zfsusers/hdesmond/Mmain"
self.mmain_path = mmain_path
if initmatch_path is None:
initmatch_path = "/mnt/extraspace/rstiskalek/csiborg/initmatch/"
self.initmatch_path = initmatch_path
self.to_new = to_new
if n_sim is not None and n_snap is not None:
self.set_info(n_sim, n_snap)
if n_sim is not None:
self.n_sim = n_sim
def __init__(self, srcdir="/mnt/extraspace/hdesmond/",
dumpdir="/mnt/extraspace/rstiskalek/csiborg/",
mmain_path="/mnt/zfsusers/hdesmond/Mmain",
initmatch_path="/mnt/extraspace/rstiskalek/csiborg/initmatch/"): # noqa
for path in [srcdir, dumpdir, mmain_path, initmatch_path]:
self._check_directory(path)
self._srcdir = srcdir
self._dumpdir = dumpdir
self._mmain_path = mmain_path
self._initmatch_path = initmatch_path
@staticmethod
def _check_directory(path):
if not isdir(path):
raise IOError("Invalid directory `{}`!".format(path))
@property
def srcdir(self):
"""
Folder where CSiBORG simulations are stored.
Path to the folder where CSiBORG simulations are stored.
Returns
-------
srcdir : int
path : str
"""
return self._srcdir
@srcdir.setter
def srcdir(self, srcdir):
"""
Set `srcdir`, check that the directory exists.
"""
if not isdir(srcdir):
raise IOError("Invalid directory `{}`!".format(srcdir))
self._srcdir = srcdir
@property
def dumpdir(self):
"""
Folder where files from `run_fit_halos` are stored.
Path to the folder where files from `run_fit_halos` are stored.
Returns
-------
dumpdir : str
path : str
"""
return self._dumpdir
@property
def temp_dumpdir(self):
"""
Temporary dumping directory.
Path to a temporary dumping folder.
Returns
-------
temp_dumpdir : str
path : str
"""
fpath = join(self.dumpdir, "temp")
if not isdir(fpath):
raise IOError("Invalid directory `{}`!".format(fpath))
raise IOError("Invalid directory `{}`.".format(fpath))
return fpath
@dumpdir.setter
def dumpdir(self, dumpdir):
"""
Set `dumpdir`, check that the directory exists.
"""
if not isdir(dumpdir):
raise IOError("Invalid directory `{}`!".format(dumpdir))
self._dumpdir = dumpdir
@property
def mmain_path(self):
"""
Path where mmain files are stored.
Path to the folder where mmain files are stored.
Returns
-------
mmain_path : str
path : str
"""
return self._mmain_path
@mmain_path.setter
def mmain_path(self, mmain_path):
"""
Set `mmain_path`, check that the directory exists.
"""
if not isdir(mmain_path):
raise IOError("Invalid directory `{}`!".format(mmain_path))
self._mmain_path = mmain_path
@property
def initmatch_path(self):
"""
Path to where match between the first and final snapshot is stored.
Path to the folder where the match between the first and final
snapshot is stored.
Returns
-------
initmach_path : str
path : str
"""
return self._initmatch_path
@initmatch_path.setter
def initmatch_path(self, initmatch_path):
def ic_ids(self, tonew):
"""
Set `initmatch_path`, check that the directory exists.
"""
if not isdir(initmatch_path):
raise IOError("Invalid directory `{}`!".format(initmatch_path))
self._initmatch_path = initmatch_path
@property
def to_new(self):
"""
Flag whether paths should point to `new` files, for example
`ramses_out_8452_new`.
Returns
-------
to_new : bool
"""
return self._to_new
@to_new.setter
def to_new(self, to_new):
"""Set `to_new`."""
if not isinstance(to_new, bool):
raise TypeError("`to_new` must be be a bool")
self._to_new = to_new
@property
def n_sim(self):
"""
The IC realisation index set by the user.
Returns
-------
n_sim : int
"""
if self._n_sim is None:
raise ValueError(
"`self.n_sim` is not set! Either provide a value directly "
"or set it using `self.set_info(...)`")
return self._n_sim
@n_sim.setter
def n_sim(self, n_sim):
"""Set `n_sim`, ensure it is a valid simulation index."""
if n_sim not in self.ic_ids:
raise ValueError(
"`{}` is not a valid IC realisation index.".format(n_sim))
self._n_sim = n_sim
@property
def n_snap(self):
"""
The snapshot index of a IC realisation set by the user.
Returns
-------
n_snap: int
"""
if self._n_snap is None:
raise ValueError(
"`self.n_sim` is not set! Either provide a value directly "
"or set it using `self.set_info(...)`")
return self._n_snap
@n_snap.setter
def n_snap(self, n_snap):
"""Set `n_snap`."""
self._n_snap = n_snap
def set_info(self, n_sim, n_snap):
"""
Convenience function for setting `n_sim` and `n_snap`.
CSiBORG IC realisation IDs from the list of folders in `self.srcdir`.
Parameters
----------
n_sim : int
CSiBORG IC realisation index.
n_snap : int
Snapshot index.
"""
self.n_sim = n_sim
if n_snap not in self.get_snapshots(n_sim):
raise ValueError(
"Invalid snapshot number `{}` for IC realisation `{}`."
.format(n_snap, n_sim))
self.n_snap = n_snap
def reset_info(self):
"""
Reset `self.n_sim` and `self.n_snap`.
"""
self._n_sim = None
self._n_snap = None
def get_n_sim(self, n_sim):
"""
Get `n_sim`. If `self.n_sim` return it, otherwise returns `n_sim`.
"""
if n_sim is None:
return self.n_sim
return n_sim
def get_n_snap(self, n_snap):
"""
Get `n_snap`. If `self.n_snap` return it, otherwise returns `n_snap`.
"""
if n_snap is None:
return self.n_snap
return n_snap
@property
def ic_ids(self):
"""
CSiBORG initial condition (IC) simulation IDs from the list of folders
in `self.srcdir`.
Returns
-------
ids : 1-dimensional array
Array of CSiBORG simulation IDs.
"""
if self.to_new:
return self._ic_ids_new
return self._ic_ids
@property
def _ic_ids(self):
"""
IC simulation IDs.
tonew : bool
If `True`, path to the '_new' ICs is returned.
Returns
-------
ids : 1-dimensional array
"""
files = glob(join(self.srcdir, "ramses_out*"))
# Select only file names
files = [f.split("/")[-1] for f in files]
# Remove files with inverted ICs
files = [f for f in files if "_inv" not in f]
# Remove the new files with z = 70 only
files = [f for f in files if "_new" not in f]
# Remove the filename with _old
files = [f for f in files if "OLD" not in f]
files = [f.split("/")[-1] for f in files] # Select only file names
if tonew:
files = [f for f in files if "_new" in f]
ids = [int(f.split("_")[2]) for f in files] # Take the IC IDs
else:
files = [f for f in files if "_inv" not in f] # Remove inv. ICs
files = [f for f in files if "_new" not in f] # Remove _new
files = [f for f in files if "OLD" not in f] # Remove _old
ids = [int(f.split("_")[-1]) for f in files]
try:
ids.remove(5511)
@ -329,111 +153,56 @@ class CSiBORGPaths:
pass
return numpy.sort(ids)
@property
def _ic_ids_new(self):
def ic_path(self, nsim, tonew=False):
"""
ICs simulation IDs denoted as `new` with recoved :math:`z = 70`
particle information.
Returns
-------
ids : 1-dimensional array
"""
files = glob(join(self.srcdir, "ramses_out*"))
# Select only file names
files = [f.split("/")[-1] for f in files]
# Only _new files
files = [f for f in files if "_new" in f]
# Take the ICs
ids = [int(f.split("_")[2]) for f in files]
return numpy.sort(ids)
def ic_path(self, n_sim=None):
"""
Path to `n_sim`th CSiBORG IC realisation.
Path to a CSiBORG IC realisation folder.
Parameters
----------
n_sim : int, optional
The index of the initial conditions (IC) realisation. By default
`None` and the set value is attempted to be used.
nsim : int
IC realisation index.
tonew : bool, optional
Whether to return the path to the '_new' IC realisation.
Returns
-------
path : str
"""
n_sim = self.get_n_sim(n_sim)
fname = "ramses_out_{}"
if self.to_new:
if tonew:
fname += "_new"
return join(self.srcdir, fname.format(n_sim))
return join(self.srcdir, fname.format(nsim))
def get_snapshots(self, n_sim=None):
def get_snapshots(self, nsim):
"""
List of snapshots for the `n_sim`th IC realisation.
List of available snapshots of a CSiBORG IC realisation.
Parameters
----------
n_sim : int
The index of the initial conditions (IC) realisation. By default
`None` and the set value is attempted to be used.
nsim : int
IC realisation index.
Returns
-------
snapshots : 1-dimensional array
"""
n_sim = self.get_n_sim(n_sim)
simpath = self.ic_path(n_sim)
simpath = self.ic_path(nsim, tonew=False)
# Get all files in simpath that start with output_
snaps = glob(join(simpath, "output_*"))
# Take just the last _00XXXX from each file and strip zeros
snaps = [int(snap.split('_')[-1].lstrip('0')) for snap in snaps]
return numpy.sort(snaps)
def get_maximum_snapshot(self, n_sim=None):
"""
Return the maximum snapshot of an IC realisation.
Parameters
----------
n_sim : int
The index of the initial conditions (IC) realisation. By default
`None` and the set value is attempted to be used.
Returns
-------
maxsnap : float
"""
n_sim = self.get_n_sim(n_sim)
return max(self.get_snapshots(n_sim))
def get_minimum_snapshot(self, n_sim=None):
"""
Return the maximum snapshot of an IC realisation.
Parameters
----------
n_sim : int
The index of the initial conditions (IC) realisation. By default
`None` and the set value is attempted to be used.
Returns
-------
minsnap : float
"""
n_sim = self.get_n_sim(n_sim)
return min(self.get_snapshots(n_sim))
def clump0_path(self, nsim):
"""
Path to a single dumped clump's particles. This is expected to point
to a dictonary whose keys are the clump indices and items structured
Path to a single dumped clump's particles. Expected to point to a
dictonary whose keys are the clump indices and items structured
arrays with the clump's particles in the initial snapshot.
Parameters
----------
nsim : int
Index of the initial conditions (IC) realisation.
IC realisation index.
Returns
-------
@ -442,27 +211,25 @@ class CSiBORGPaths:
cdir = join(self.dumpdir, "initmatch")
return join(cdir, "clump_{}_{}.npy".format(nsim, "particles"))
def snapshot_path(self, n_snap=None, n_sim=None):
def snapshot_path(self, nsnap, nsim, tonew=False):
"""
Path to a CSiBORG IC realisation snapshot.
Parameters
----------
n_snap : int
Snapshot index. By default `None` and the set value is attempted
to be used.
n_sim : str
Corresponding CSiBORG IC realisation index. By default `None` and
the set value is attempted to be used.
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
tonew : bool, optional
Whether to return the path to the '_new' IC realisation.
Returns
-------
snappath : str
"""
n_snap = self.get_n_snap(n_snap)
n_sim = self.get_n_sim(n_sim)
simpath = self.ic_path(n_sim)
return join(simpath, "output_{}".format(str(n_snap).zfill(5)))
simpath = self.ic_path(nsim, tonew=tonew)
return join(simpath, "output_{}".format(str(nsnap).zfill(5)))
###############################################################################
@ -472,56 +239,41 @@ class CSiBORGPaths:
class ParticleReader:
"""
Tools to read in particle files alon with their corresponding clumps.
Shortcut to read in particle files along with their corresponding clumps.
Parameters
----------
paths : py:class`csiborgtools.read.CSiBORGPaths`
CSiBORG paths-handling object with set `n_sim` and `n_snap`.
"""
_paths = None
def __init__(self, paths):
self.paths = paths
assert isinstance(paths, CSiBORGPaths)
self._paths = paths
@property
def paths(self):
"""
The paths-handling object.
Returns
-------
paths : :py:class:`csiborgtools.read.CSiBORGPaths`
"""
return self._paths
@paths.setter
def paths(self, paths):
"""
Set `paths`. Makes sure it is the right object and `n_sim` and `n_snap`
are both set.
"""
if not isinstance(paths, CSiBORGPaths):
raise TypeError("`paths` must be of type `CSiBORGPaths`.")
if paths.n_sim is None or paths.n_snap is None:
raise ValueError(
"`paths` must have set both `n_sim` and `n_snap`!")
self._paths = paths
def read_info(self):
def read_info(self, nsnap, nsim):
"""
Read CSiBORG simulation snapshot info.
Parameters
----------
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
Returns
-------
info : dict
Dictionary of info paramaters. Note that both keys and values are
strings.
Dictionary of information paramaters. Note that both keys and
values are strings.
"""
# Open the info file
n_snap = self.paths.n_snap
snappath = self.paths.snapshot_path()
filename = join(snappath, "info_{}.txt".format(str(n_snap).zfill(5)))
snappath = self.paths.snapshot_path(nsnap, nsim)
filename = join(snappath, "info_{}.txt".format(str(nsnap).zfill(5)))
with open(filename, "r") as f:
info = f.read().split()
# Throw anything below ordering line out
@ -533,12 +285,16 @@ class ParticleReader:
vals = info[eqs + 1]
return {key: val for key, val in zip(keys, vals)}
def open_particle(self, verbose=True):
def open_particle(self, nsnap, nsim, verbose=True):
"""
Open particle files to a given CSiBORG simulation.
Parameters
----------
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
verbose : bool, optional
Verbosity flag.
@ -549,28 +305,25 @@ class ParticleReader:
partfiles : list of `scipy.io.FortranFile`
Opened part files.
"""
n_snap = self.paths.n_snap
# Zeros filled snapshot number and the snapshot path
nout = str(n_snap).zfill(5)
snappath = self.paths.snapshot_path()
snappath = self.paths.snapshot_path(nsnap, nsim)
ncpu = int(self.read_info()["ncpu"])
nsnap = str(nsnap).zfill(5)
if verbose:
print("Reading in output `{}` with ncpu = `{}`."
.format(nout, ncpu))
.format(nsnap, ncpu))
# First read the headers. Reallocate arrays and fill them.
nparts = numpy.zeros(ncpu, dtype=int)
partfiles = [None] * ncpu
for cpu in range(ncpu):
cpu_str = str(cpu + 1).zfill(5)
fpath = join(snappath, "part_{}.out{}".format(nout, cpu_str))
fpath = join(snappath, "part_{}.out{}".format(nsnap, cpu_str))
f = FortranFile(fpath)
# Read in this order
ncpuloc = f.read_ints()
if ncpuloc != ncpu:
infopath = join(snappath, "info_{}.txt".format(nout))
infopath = join(snappath, "info_{}.txt".format(nsnap))
raise ValueError(
"`ncpu = {}` of `{}` disagrees with `ncpu = {}` "
"of `{}`.".format(ncpu, infopath, ncpuloc, fpath))
@ -590,8 +343,7 @@ class ParticleReader:
@staticmethod
def read_sp(dtype, partfile):
"""
Utility function to read a single particle file, depending on its
dtype.
Utility function to read a single particle file.
Parameters
----------
@ -609,9 +361,9 @@ class ParticleReader:
simpath : str
The complete path to the CSiBORG simulation.
"""
if dtype in [F16, F32, F64]:
if dtype in [numpy.float16, numpy.float32, numpy.float64]:
return partfile.read_reals('d')
elif dtype in [I32]:
elif dtype in [numpy.int32]:
return partfile.read_ints()
else:
raise TypeError("Unexpected dtype `{}`.".format(dtype))
@ -620,7 +372,8 @@ class ParticleReader:
def nparts_to_start_ind(nparts):
"""
Convert `nparts` array to starting indices in a pre-allocated array for
looping over the CPU number.
looping over the CPU number. The starting indices calculated as a
cumulative sum starting at 0.
Parameters
----------
@ -630,40 +383,40 @@ class ParticleReader:
Returns
-------
start_ind : 1-dimensional array
The starting indices calculated as a cumulative sum starting at 0.
"""
return numpy.hstack([[0], numpy.cumsum(nparts[:-1])])
def read_particle(self, pars_extract, verbose=True):
def read_particle(self, nsnap, nsim, pars_extract, verbose=True):
"""
Read particle files of a simulation at a given snapshot and return
values of `pars_extract`.
Parameters
----------
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
pars_extract : list of str
Parameters to be extacted.
Nsnap : int
The index of the redshift snapshot.
simpath : str
The complete path to the CSiBORG simulation.
verbose : bool, optional
Verbosity flag while for reading the CPU outputs.
Returns
-------
out : structured array
The data read from the particle file.
"""
# Open the particle files
nparts, partfiles = self.open_particle(verbose=verbose)
nparts, partfiles = self.open_particle(nsnap, nsim, verbose=verbose)
if verbose:
print("Opened {} particle files.".format(nparts.size))
ncpu = nparts.size
# Order in which the particles are written in the FortranFile
forder = [("x", F32), ("y", F32), ("z", F32),
("vx", F32), ("vy", F32), ("vz", F32),
("M", F32), ("ID", I32), ("level", I32)]
forder = [("x", numpy.float32), ("y", numpy.float32),
("z", numpy.float32), ("vx", numpy.float32),
("vy", numpy.float32), ("vz", numpy.float32),
("M", numpy.float32), ("ID", numpy.int32),
("level", numpy.int32)]
fnames = [fp[0] for fp in forder]
fdtypes = [fp[1] for fp in forder]
# Check there are no strange parameters
@ -677,7 +430,7 @@ class ParticleReader:
npart_tot = numpy.sum(nparts)
# A dummy array is necessary for reading the fortran files.
dum = numpy.full(npart_tot, numpy.nan, dtype=F16)
dum = numpy.full(npart_tot, numpy.nan, dtype=numpy.float16)
# These are the data we read along with types
dtype = {"names": pars_extract,
"formats": [forder[fnames.index(p)][1] for p in pars_extract]}
@ -699,13 +452,17 @@ class ParticleReader:
return out
def open_unbinding(self, cpu):
def open_unbinding(self, nsnap, nsim, cpu):
"""
Open particle files to a given CSiBORG simulation. Note that to be
consistent CPU is incremented by 1.
Parameters
----------
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
cpu : int
The CPU index.
@ -714,18 +471,23 @@ class ParticleReader:
unbinding : `scipy.io.FortranFile`
The opened unbinding FortranFile.
"""
nout = str(self.paths.n_snap).zfill(5)
nsnap = str(nsnap).zfill(5)
cpu = str(cpu + 1).zfill(5)
fpath = join(self.paths.ic_path(), "output_{}".format(nout),
"unbinding_{}.out{}".format(nout, cpu))
fpath = join(self.paths.ic_path(nsim, to_new=False),
"output_{}".format(nsnap),
"unbinding_{}.out{}".format(nsnap, cpu))
return FortranFile(fpath)
def read_clumpid(self, verbose=True):
def read_clumpid(self, nsnap, nsim, verbose=True):
"""
Read clump IDs of particles from unbinding files.
Parameters
----------
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
verbose : bool, optional
Verbosity flag while for reading the CPU outputs.
@ -734,16 +496,16 @@ class ParticleReader:
clumpid : 1-dimensional array
The array of clump IDs.
"""
nparts, __ = self.open_particle(verbose)
nparts, __ = self.open_particle(nsnap, nsim, verbose)
start_ind = self.nparts_to_start_ind(nparts)
ncpu = nparts.size
clumpid = numpy.full(numpy.sum(nparts), numpy.nan, dtype=I32)
clumpid = numpy.full(numpy.sum(nparts), numpy.nan, dtype=numpy.int32)
iters = tqdm(range(ncpu)) if verbose else range(ncpu)
for cpu in iters:
i = start_ind[cpu]
j = nparts[cpu]
ff = self.open_unbinding(cpu)
ff = self.open_unbinding(nsnap, nsim, cpu)
clumpid[i:i + j] = ff.read_ints()
# Close
ff.close()
@ -772,12 +534,17 @@ class ParticleReader:
mask = clump_ids != 0
return clump_ids[mask], particles[mask]
def read_clumps(self, cols=None):
def read_clumps(self, nsnap, nsim, cols=None):
"""
Read in a clump file `clump_Nsnap.dat`.
Parameters
----------
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
cols : list of str, optional.
Columns to extract. By default `None` and all columns are
extracted.
@ -787,9 +554,10 @@ class ParticleReader:
out : structured array
Structured array of the clumps.
"""
n_snap = str(self.paths.n_snap).zfill(5)
fname = join(self.paths.ic_path(), "output_{}".format(n_snap),
"clump_{}.dat".format(n_snap))
nsnap = str(nsnap).zfill(5)
fname = join(self.paths.ic_path(nsim, to_new=False),
"output_{}".format(nsnap),
"clump_{}.dat".format(nsnap))
# Check the file exists.
if not isfile(fname):
raise FileExistsError(
@ -797,10 +565,12 @@ class ParticleReader:
# Read in the clump array. This is how the columns must be written!
data = numpy.genfromtxt(fname)
clump_cols = [("index", I64), ("level", I64), ("parent", I64),
("ncell", F64), ("peak_x", F64), ("peak_y", F64),
("peak_z", F64), ("rho-", F64), ("rho+", F64),
("rho_av", F64), ("mass_cl", F64), ("relevance", F64)]
clump_cols = [("index", numpy.int64), ("level", numpy.int64),
("parent", numpy.int64), ("ncell", numpy.float64),
("peak_x", numpy.float64), ("peak_y", numpy.float64),
("peak_z", numpy.float64), ("rho-", numpy.float64),
("rho+", numpy.float64), ("rho_av", numpy.float64),
("mass_cl", numpy.float64), ("relevance", numpy.float64)]
out0 = cols_to_structured(data.shape[0], clump_cols)
for i, name in enumerate(out0.dtype.names):
out0[name] = data[:, i]
@ -824,19 +594,19 @@ class ParticleReader:
return out
def read_mmain(n, srcdir, fname="Mmain_{}.npy"):
def read_mmain(nsim, srcdir, fname="Mmain_{}.npy"):
"""
Read `mmain` numpy arrays of central halos whose mass contains their
substracture contribution.
Parameters
----------
n : int
The index of the initial conditions (IC) realisation.
nsim : int
IC realisation index.
srcdir : str
The path to the folder containing the files.
Path to the folder containing the files.
fname : str, optional
The file name convention. By default `Mmain_{}.npy`, where the
File name convention. By default `Mmain_{}.npy`, where the
substituted value is `n`.
Returns
@ -844,11 +614,12 @@ def read_mmain(n, srcdir, fname="Mmain_{}.npy"):
out : structured array
Array with the central halo information.
"""
fpath = join(srcdir, fname.format(n))
fpath = join(srcdir, fname.format(nsim))
arr = numpy.load(fpath)
cols = [("index", I64), ("peak_x", F64), ("peak_y", F64),
("peak_z", F64), ("mass_cl", F64), ("sub_frac", F64)]
cols = [("index", numpy.int64), ("peak_x", numpy.float64),
("peak_y", numpy.float64), ("peak_z", numpy.float64),
("mass_cl", numpy.float64), ("sub_frac", numpy.float64)]
out = cols_to_structured(arr.shape[0], cols)
for i, name in enumerate(out.dtype.names):
out[name] = arr[:, i]
@ -856,26 +627,26 @@ def read_mmain(n, srcdir, fname="Mmain_{}.npy"):
return out
def read_initcm(n, srcdir, fname="clump_{}_cm.npy"):
def read_initcm(nsim, srcdir, fname="clump_{}_cm.npy"):
"""
Read `clump_cm`, i.e. the center of mass of a clump at redshift z = 70.
If the file does not exist returns `None`.
Parameters
----------
n : int
The index of the initial conditions (IC) realisation.
nsim : int
IC realisation index.
srcdir : str
The path to the folder containing the files.
Path to the folder containing the files.
fname : str, optional
The file name convention. By default `clump_cm_{}.npy`, where the
substituted value is `n`.
File name convention. By default `clump_cm_{}.npy`, where the
substituted value is `nsim`.
Returns
-------
out : structured array
"""
fpath = join(srcdir, fname.format(n))
fpath = join(srcdir, fname.format(nsim))
try:
return numpy.load(fpath)
except FileNotFoundError:

19
csiborgtools/read/summaries.py
View file

@ -15,13 +15,18 @@
"""
Tools for summarising various results.
"""
import numpy
import joblib
from os.path import (join, isfile)
from glob import glob
import numpy
import joblib
from tqdm import tqdm
###############################################################################
# PKReader #
###############################################################################
class PKReader:
"""
A shortcut object for reading in the power spectrum files.
@ -170,6 +175,11 @@ class PKReader:
return ks, xpks
###############################################################################
# PKReader #
###############################################################################
class kNNCDFReader:
"""
Shortcut object to read in the kNN CDF data.
@ -304,6 +314,11 @@ class kNNCDFReader:
return [file for file in glob(join(folder, "*")) if str(ic) in file]
###############################################################################
# PKReader #
###############################################################################
class PairOverlap:
r"""
A shortcut object for reading in the results of matching two simulations.

2
csiborgtools/units/__init__.py
View file

@ -14,4 +14,4 @@
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
from .transforms import cartesian_to_radec, radec_to_cartesian # noqa
from .box_units import (BoxUnits) # noqa
from .box_units import BoxUnits # noqa

45
csiborgtools/units/box_units.py
View file

@ -15,7 +15,6 @@
"""
Simulation box unit transformations.
"""
import numpy
from scipy.interpolate import interp1d
from astropy.cosmology import LambdaCDM
@ -38,17 +37,21 @@ class BoxUnits:
Paramaters
----------
nsnap : int
Snapshot index.
nsim : int
IC realisation index.
paths : py:class`csiborgtools.read.CSiBORGPaths`
CSiBORG paths-handling object with set `n_sim` and `n_snap`.
CSiBORG paths object.
"""
_cosmo = None
def __init__(self, paths):
def __init__(self, nsnap, nsim, paths):
"""
Read in the snapshot info file and set the units from it.
"""
partreader = ParticleReader(paths)
info = partreader.read_info()
info = partreader.read_info(nsnap, nsim)
pars = ["boxlen", "time", "aexp", "H0",
"omega_m", "omega_l", "omega_k", "omega_b",
"unit_l", "unit_d", "unit_t"]
@ -68,7 +71,6 @@ class BoxUnits:
Returns
-------
cosmo : `astropy.cosmology.LambdaCDM`
The CSiBORG cosmology.
"""
return self._cosmo
@ -81,7 +83,6 @@ class BoxUnits:
Returns
-------
H0 : float
Hubble constant.
"""
return self._H0
@ -93,7 +94,6 @@ class BoxUnits:
Returns
-------
h : float
The little h
"""
return self._H0 / 100
@ -106,7 +106,6 @@ class BoxUnits:
Returns
-------
G : float
The gravitational constant.
"""
return constants.G.cgs.value * (self._unit_d * self._unit_t ** 2)
@ -118,7 +117,6 @@ class BoxUnits:
Returns
-------
H0 : float
The Hubble constant.
"""
return self.H0 * 1e5 / units.Mpc.to(units.cm) * self._unit_t
@ -130,7 +128,6 @@ class BoxUnits:
Returns
-------
c : float
The speed of light.
"""
return constants.c.cgs.value * self._unit_t / self._unit_l
@ -142,7 +139,6 @@ class BoxUnits:
Returns
-------
rhoc : float
The critical density.
"""
return 3 * self.box_H0 ** 2 / (8 * numpy.pi * self.box_G)
@ -150,8 +146,7 @@ class BoxUnits:
def box2kpc(self, length):
r"""
Convert length from box units to :math:`\mathrm{ckpc}` (with
:math:`h=0.705`). It appears that `self.unit_l` must be in
:math:`\mathrm{cm}`.
:math:`h=0.705`).
Parameters
----------
@ -224,12 +219,12 @@ class BoxUnits:
Parameters
----------
dist : float
Dist in box units.
Distance in box units.
Returns
-------
cosmo_redshift : foat
The cosmological redshift.
Cosmological redshift.
"""
x = numpy.linspace(0., 1., 5001)
y = self.cosmo.comoving_distance(x)
@ -244,17 +239,17 @@ class BoxUnits:
Parameters
----------
vx, vy, vz : 1-dimensional arrays
The Cartesian velocity components.
Cartesian velocity components.
px, py, pz : 1-dimensional arrays
The Cartesian position vectors components.
Cartesian position vectors components.
p0x, p0y, p0z : floats
The centre of the box. By default 0, in which it is assumed that
the coordinates are already centred.
Centre of the box coordinates. By default 0, in which it is assumed
that the coordinates are already centred.
Returns
-------
pec_redshift : 1-dimensional array
The peculiar redshift.
Peculiar redshift.
"""
# Peculiar velocity along the radial distance
r = numpy.vstack([px - p0x, py - p0y, pz - p0z]).T
@ -275,17 +270,17 @@ class BoxUnits:
Parameters
----------
vx, vy, vz : 1-dimensional arrays
The Cartesian velocity components.
Cartesian velocity components.
px, py, pz : 1-dimensional arrays
The Cartesian position vectors components.
Cartesian position vectors components.
p0x, p0y, p0z : floats
The centre of the box. By default 0, in which it is assumed that
the coordinates are already centred.
Centre of the box coordinates. By default 0, in which it is assumed
that the coordinates are already centred.
Returns
-------
obs_redshift : 1-dimensional array
The observed redshift.
Observed redshift.
"""
r = numpy.vstack([px - p0x, py - p0y, pz - p0z]).T
zcosmo = self.box2cosmoredshift(numpy.sum(r**2, axis=1)**0.5)

3
csiborgtools/units/transforms.py
View file

@ -15,7 +15,6 @@
"""
Various coordinate transformations.
"""
import numpy
@ -49,7 +48,7 @@ def radec_to_cartesian(dist, ra, dec, isdeg=True):
Parameters
----------
dist, ra, dec : 1-dimensional arrays
The spherical coordinates.
Spherical coordinates.
isdeg : bool, optional
Whether `ra` and `dec` are in degres. By default `True`.

18
csiborgtools/utils/recarray_manip.py
View file

@ -15,8 +15,6 @@
"""
Utilility functions for manipulation structured arrays.
"""
import numpy
@ -51,7 +49,7 @@ def add_columns(arr, X, cols):
Parameters
----------
arr : record array
The record array to add columns to.
Record array to add columns to.
X : (list of) 1-dimensional array(s) or 2-dimensional array
Columns to be added.
cols : str or list of str
@ -60,7 +58,6 @@ def add_columns(arr, X, cols):
Returns
-------
out : record array
The new record array with added values.
"""
# Make sure cols is a list of str and X a 2D array
cols = [cols] if isinstance(cols, str) else cols
@ -95,14 +92,13 @@ def rm_columns(arr, cols):
Parameters
----------
arr : record array
The record array to remove columns from.
Record array to remove columns from.
cols : str or list of str
Column names to be removed.
Returns
-------
out : record array
Record array with removed columns.
"""
# Check columns we wish to delete are in the array
cols = [cols] if isinstance(cols, str) else cols
@ -173,7 +169,7 @@ def array_to_structured(arr, cols):
Returns
-------
out : structured array
The output structured array.
Output structured array.
"""
cols = [cols] if isinstance(cols, str) else cols
if arr.ndim != 2 and arr.shape[1] != len(cols):
@ -196,15 +192,15 @@ def flip_cols(arr, col1, col2):
Parameters
----------
arr : structured array
The array whose columns are to be converted.
Array whose columns are to be converted.
col1 : str
The first column name.
First column name.
col2 : str
The second column name.
Second column name.
Returns
-------
nothing
None
"""
dum = numpy.copy(arr[col1])
arr[col1] = arr[col2]

4242
notebooks/knn.ipynb
View file

File diff suppressed because one or more lines are too long

68
scripts/run_asciipos.py
View file

@ -1,68 +0,0 @@
# 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 os.path import join, exists
from argparse import ArgumentParser
import numpy
from datetime import datetime
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!")

13
scripts/run_asciipos.sh
View file

@ -1,13 +0,0 @@
nthreads=1
memory=75
queue="berg"
env="/mnt/zfsusers/rstiskalek/csiborgtools/venv_galomatch/bin/python"
file="run_asciipos.py"
mode="dump"
cm="addqueue -q $queue -n $nthreads -m $memory $env $file --mode $mode"
echo "Submitting:"
echo $cm
echo
$cm

71
scripts/run_crossmatch.py
View file

@ -1,71 +0,0 @@
# 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.
"""
MPI script to run the CSiBORG realisations matcher.
TODO
----
- [ ] Update this script
"""
import numpy
from datetime import datetime
from mpi4py import MPI
from os.path import join
try:
import csiborgtools
except ModuleNotFoundError:
import sys
sys.path.append("../")
import csiborgtools
import utils
# Get MPI things
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
nproc = comm.Get_size()
# File paths
fperm = join(utils.dumpdir, "overlap", "cross_{}.npy")
# fperm = join(utils.dumpdir, "match", "cross_matches.npy")
# Set up the catalogue
paths = csiborgtools.read.CSiBORGPaths(to_new=False)
print("{}: started reading in the combined catalogue.".format(datetime.now()),
flush=True)
cat = csiborgtools.read.CombinedHaloCatalogue(
paths, min_m500=None, max_dist=None, verbose=False)
print("{}: finished reading in the combined catalogue with `{}`."
.format(datetime.now(), cat.n_sims), flush=True)
matcher = csiborgtools.match.RealisationsMatcher(cat)
for i in csiborgtools.fits.split_jobs(len(cat.n_sims), nproc)[rank]:
n = cat.n_sims[i]
print("{}: rank {} working on simulation `{}`."
.format(datetime.now(), rank, n), flush=True)
out = matcher.cross_knn_position_single(
i, nmult=15, dlogmass=2, init_dist=True, overlap=False, verbose=False,
overlapper_kwargs={"smooth_scale": 1})
# Dump the result
fout = fperm.format(n)
print("Saving results to `{}`.".format(fout))
with open(fout, "wb") as f:
numpy.save(fout, out)
comm.Barrier()
if rank == 0:
print("All finished.")

17
scripts/run_crossmatch.sh
View file

@ -1,17 +0,0 @@
nthreads=1
memory=32
queue="berg"
env="/mnt/zfsusers/rstiskalek/csiborgtools/venv_galomatch/bin/python"
file="run_crossmatch.py"
pythoncm="$env $file"
# echo "Submitting:"
# echo $pythoncm
# echo
# $pythoncm
cm="addqueue -q $queue -n $nthreads -m $memory $pythoncm"
echo "Submitting:"
echo $cm
echo
$cm

43
scripts/run_crosspk.py
View file

@ -16,14 +16,14 @@
MPI script to calculate the matter cross power spectrum between CSiBORG
IC realisations. Units are Mpc/h.
"""
from gc import collect
from argparse import ArgumentParser
from os import remove
from os.path import join
from itertools import combinations
from datetime import datetime
import numpy
import joblib
from datetime import datetime
from itertools import combinations
from os.path import join
from os import remove
from gc import collect
from mpi4py import MPI
import Pk_library as PKL
try:
@ -32,9 +32,9 @@ except ModuleNotFoundError:
import sys
sys.path.append("../")
import csiborgtools
import utils
dumpdir = "/mnt/extraspace/rstiskalek/csiborg/"
parser = ArgumentParser()
parser.add_argument("--grid", type=int)
parser.add_argument("--halfwidth", type=float, default=0.5)
@ -47,27 +47,24 @@ nproc = comm.Get_size()
MAS = "CIC" # mass asignment scheme
paths = csiborgtools.read.CSiBORGPaths()
ics = paths.ic_ids
n_sims = len(ics)
box = csiborgtools.units.BoxUnits(paths)
reader = csiborgtools.read.ParticleReader(paths)
ics = paths.ic_ids(tonew=False)
nsims = len(ics)
# File paths
ftemp = join(utils.dumpdir, "temp_crosspk",
ftemp = join(dumpdir, "temp_crosspk",
"out_{}_{}" + "_{}".format(args.halfwidth))
fout = join(utils.dumpdir, "crosspk",
fout = join(dumpdir, "crosspk",
"out_{}_{}" + "_{}.p".format(args.halfwidth))
jobs = csiborgtools.fits.split_jobs(n_sims, nproc)[rank]
jobs = csiborgtools.fits.split_jobs(nsims, nproc)[rank]
for n in jobs:
print("Rank {}@{}: saving {}th delta.".format(rank, datetime.now(), n))
# Set the paths
n_sim = ics[n]
paths.set_info(n_sim, paths.get_maximum_snapshot(n_sim))
# Set reader and the box
reader = csiborgtools.read.ParticleReader(paths)
box = csiborgtools.units.BoxUnits(paths)
# Read particles
particles = reader.read_particle(["x", "y", "z", "M"], verbose=False)
nsim = ics[n]
particles = reader.read_particle(max(paths.get_snapshots(nsim)), nsim,
["x", "y", "z", "M"], verbose=False)
# Halfwidth -- particle selection
if args.halfwidth < 0.5:
particles = csiborgtools.read.halfwidth_select(
@ -85,9 +82,9 @@ for n in jobs:
collect()
# Dump the results
with open(ftemp.format(n_sim, "delta") + ".npy", "wb") as f:
with open(ftemp.format(nsim, "delta") + ".npy", "wb") as f:
numpy.save(f, delta)
joblib.dump([aexp, length], ftemp.format(n_sim, "lengths") + ".p")
joblib.dump([aexp, length], ftemp.format(nsim, "lengths") + ".p")
# Try to clean up memory
del delta
@ -97,8 +94,8 @@ for n in jobs:
comm.Barrier()
# Get off-diagonal elements and append the diagoal
combs = [c for c in combinations(range(n_sims), 2)]
for i in range(n_sims):
combs = [c for c in combinations(range(nsims), 2)]
for i in range(nsims):
combs.append((i, i))
prev_delta = [-1, None, None, None] # i, delta, aexp, length

41
scripts/run_fieldprop.py
View file

@ -16,11 +16,11 @@
MPI script to evaluate field properties at the galaxy positions.
"""
from argparse import ArgumentParser
import numpy
from datetime import datetime
from mpi4py import MPI
from os.path import join
from os import remove
from datetime import datetime
import numpy
from mpi4py import MPI
try:
import csiborgtools
except ModuleNotFoundError:
@ -29,6 +29,7 @@ except ModuleNotFoundError:
import csiborgtools
import utils
dumpdir = "/mnt/extraspace/rstiskalek/csiborg/"
parser = ArgumentParser()
parser.add_argument("--survey", type=str, choices=["SDSS"])
parser.add_argument("--grid", type=int)
@ -52,30 +53,28 @@ pos = pos.astype(numpy.float32)
# File paths
fname = "out_{}_{}_{}_{}_{}".format(
survey.name, args.grid, args.MAS, args.halfwidth, args.smooth_scale)
ftemp = join(utils.dumpdir, "temp_fields", fname + "_{}.npy")
fperm = join(utils.dumpdir, "fields", fname + ".npy")
ftemp = join(dumpdir, "temp_fields", fname + "_{}.npy")
fperm = join(dumpdir, "fields", fname + ".npy")
# Edit depending on what is calculated
dtype = {"names": ["delta", "phi"], "formats": [numpy.float32] * 2}
# CSiBORG simulation paths
paths = csiborgtools.read.CSiBORGPaths()
ics = paths.ic_ids
n_sims = len(ics)
ics = paths.ic_ids(tonew=False)
nsims = len(ics)
for n in csiborgtools.fits.split_jobs(n_sims, nproc)[rank]:
for n in csiborgtools.fits.split_jobs(nsims, nproc)[rank]:
print("Rank {}@{}: working on {}th IC.".format(rank, datetime.now(), n),
flush=True)
# Set the paths
n_sim = ics[n]
paths.set_info(n_sim, paths.get_maximum_snapshot(n_sim))
# Set reader and the box
nsim = ics[n]
nsnap = max(paths.get_snapshots(nsim))
reader = csiborgtools.read.ParticleReader(paths)
box = csiborgtools.units.BoxUnits(paths)
box = csiborgtools.units.BoxUnits(nsnap, nsim, paths)
# Read particles and select a subset of them
particles = reader.read_particle(["x", "y", "z", "M"], verbose=False)
particles = reader.read_particle(nsnap, nsim, ["x", "y", "z", "M"],
verbose=False)
if args.halfwidth < 0.5:
particles = csiborgtools.read.halfwidth_select(
args.halfwidth, particles)
@ -101,7 +100,7 @@ for n in csiborgtools.fits.split_jobs(n_sims, nproc)[rank]:
# ...
# Dump the results
with open(ftemp.format(n_sim), "wb") as f:
with open(ftemp.format(nsim), "wb") as f:
numpy.save(f, out)
# Wait for all ranks to finish
@ -109,16 +108,16 @@ comm.Barrier()
if rank == 0:
print("Collecting files...", flush=True)
out = numpy.full((n_sims, pos.shape[0]), numpy.nan, dtype=dtype)
out = numpy.full((nsims, pos.shape[0]), numpy.nan, dtype=dtype)
for n in range(n_sims):
n_sim = ics[n]
with open(ftemp.format(n_sim), "rb") as f:
for n in range(nsims):
nsim = ics[n]
with open(ftemp.format(nsim), "rb") as f:
fin = numpy.load(f, allow_pickle=True)
for name in dtype["names"]:
out[name][n, ...] = fin[name]
# Remove the temporary file
remove(ftemp.format(n_sim))
remove(ftemp.format(nsim))
print("Saving results to `{}`.".format(fperm), flush=True)
with open(fperm, "wb") as f:

69
scripts/run_fit_halos.py
View file

@ -14,12 +14,11 @@
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
"""
A script to fit halos (concentration, ...). The particle array of each CSiBORG
realisation must have been split in advance by `run_split_halos`.
realisation must have been split in advance by `runsplit_halos`.
"""
import numpy
from datetime import datetime
from os.path import join
from datetime import datetime
import numpy
from mpi4py import MPI
try:
import csiborgtools
@ -29,47 +28,48 @@ except ModuleNotFoundError:
import csiborgtools
import utils
F64 = numpy.float64
I64 = numpy.int64
# Get MPI things
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
nproc = comm.Get_size()
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)]
paths = csiborgtools.read.CSiBORGPaths()
dumpdir = "/mnt/extraspace/rstiskalek/csiborg/"
loaddir = join(dumpdir, "temp")
cols_collect = [("npart", numpy.int64), ("totpartmass", numpy.float64),
("Rs", numpy.float64), ("vx", numpy.float64),
("vy", numpy.float64), ("vz", numpy.float64),
("Lx", numpy.float64), ("Ly", numpy.float64),
("Lz", numpy.float64), ("rho0", numpy.float64),
("conc", numpy.float64), ("rmin", numpy.float64),
("rmax", numpy.float64), ("r200", numpy.float64),
("r500", numpy.float64), ("m200", numpy.float64),
("m500", numpy.float64), ("lambda200c", numpy.float64)]
for i, n_sim in enumerate(paths.ic_ids):
for i, nsim in enumerate(paths.ic_ids(tonew=False)):
if rank == 0:
print("{}: calculating {}th simulation.".format(datetime.now(), i))
# Correctly set the paths!
n_snap = paths.get_maximum_snapshot(n_sim)
paths.set_info(n_sim, n_snap)
box = csiborgtools.units.BoxUnits(paths)
nsnap = max(paths.get_snapshots(nsim))
box = csiborgtools.units.BoxUnits(nsnap, nsim, paths)
jobs = csiborgtools.fits.split_jobs(utils.Nsplits, nproc)[rank]
for n_split in jobs:
for nsplit in jobs:
parts, part_clumps, clumps = csiborgtools.fits.load_split_particles(
n_split, paths, remove_split=False)
nsplit, nsnap, nsim, paths, remove_split=False)
N = clumps.size
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)]
cols = [("index", numpy.int64), ("npart", numpy.int64),
("totpartmass", numpy.float64), ("Rs", numpy.float64),
("rho0", numpy.float64), ("conc", numpy.float64),
("lambda200c", numpy.float64), ("vx", numpy.float64),
("vy", numpy.float64), ("vz", numpy.float64),
("Lx", numpy.float64), ("Ly", numpy.float64),
("Lz", numpy.float64), ("rmin", numpy.float64),
("rmax", numpy.float64), ("r200", numpy.float64),
("r500", numpy.float64), ("m200", numpy.float64),
("m500", numpy.float64)]
out = csiborgtools.utils.cols_to_structured(N, cols)
out["index"] = clumps["index"]
@ -106,7 +106,7 @@ for i, n_sim in enumerate(paths.ic_ids):
out["rho0"][n] = nfwpost.rho0_from_Rs(Rs)
out["conc"][n] = out["r200"][n] / Rs
csiborgtools.read.dump_split(out, n_split, paths)
csiborgtools.read.dump_split(out, nsplit, nsnap, nsim, paths)
# Wait until all jobs finished before moving to another simulation
comm.Barrier()
@ -116,11 +116,10 @@ for i, n_sim in enumerate(paths.ic_ids):
print("Collecting results!")
partreader = csiborgtools.read.ParticleReader(paths)
out_collected = csiborgtools.read.combine_splits(
utils.Nsplits, partreader, cols_collect, remove_splits=True,
verbose=False)
utils.Nsplits, nsnap, nsim, partreader, cols_collect,
remove_splits=True, verbose=False)
fname = join(paths.dumpdir, "ramses_out_{}_{}.npy"
.format(str(paths.n_sim).zfill(5),
str(paths.n_snap).zfill(5)))
.format(str(nsim).zfill(5), str(nsnap).zfill(5)))
print("Saving results to `{}`.".format(fname))
numpy.save(fname, out_collected)

33
scripts/run_initmatch.py
View file

@ -19,14 +19,14 @@ are grouped in a clump at present redshift.
Optionally also dumps the clumps information, however watch out as this will
eat up a lot of memory.
"""
import numpy
from argparse import ArgumentParser
from distutils.util import strtobool
from datetime import datetime
from mpi4py import MPI
from gc import collect
from os.path import join
from os import remove
from gc import collect
from argparse import ArgumentParser
from datetime import datetime
from distutils.util import strtobool
import numpy
from mpi4py import MPI
try:
import csiborgtools
except ModuleNotFoundError:
@ -44,9 +44,8 @@ parser = ArgumentParser()
parser.add_argument("--dump_clumps", type=lambda x: bool(strtobool(x)))
args = parser.parse_args()
init_paths = csiborgtools.read.CSiBORGPaths(to_new=True)
fin_paths = csiborgtools.read.CSiBORGPaths(to_new=False)
nsims = init_paths.ic_ids
paths = csiborgtools.read.CSiBORGPaths()
nsims = paths.ic_ids(tonew=True)
# Output files
dumpdir = "/mnt/extraspace/rstiskalek/csiborg/"
@ -58,22 +57,18 @@ for nsim in nsims:
if rank == 0:
print("{}: reading simulation {}.".format(datetime.now(), nsim),
flush=True)
# Set the snapshot numbers
init_paths.set_info(nsim, init_paths.get_minimum_snapshot(nsim))
fin_paths.set_info(nsim, fin_paths.get_maximum_snapshot(nsim))
# Set the readers
init_reader = csiborgtools.read.ParticleReader(init_paths)
fin_reader = csiborgtools.read.ParticleReader(fin_paths)
nsnap_min = min(paths.get_snapshots(nsim))
nsnap_max = max(paths.get_snapshots(nsim))
reader = csiborgtools.read.ParticleReader(paths)
# Read and sort the initial particle files by their particle IDs
part0 = init_reader.read_particle(["x", "y", "z", "M", "ID"],
part0 = reader.read_particle(nsnap_min, nsim, ["x", "y", "z", "M", "ID"],
verbose=False)
part0 = part0[numpy.argsort(part0["ID"])]
# Order the final snapshot clump IDs by the particle IDs
pid = fin_reader.read_particle(["ID"], verbose=False)["ID"]
clump_ids = fin_reader.read_clumpid(verbose=False)
pid = reader.read_particle(nsnap_max, nsim, ["ID"], verbose=False)["ID"]
clump_ids = reader.read_clumpid(nsnap_max, nsim, verbose=False)
clump_ids = clump_ids[numpy.argsort(pid)]
del pid

13
scripts/run_knn.py
View file

@ -62,6 +62,7 @@ ics = [7444, 7468, 7492, 7516, 7540, 7564, 7588, 7612, 7636, 7660, 7684,
dumpdir = "/mnt/extraspace/rstiskalek/csiborg/knn"
fout_auto = join(dumpdir, "auto", "knncdf_{}.p")
fout_cross = join(dumpdir, "cross", "knncdf_{}_{}.p")
paths = csiborgtools.read.CSiBORGPaths()
###############################################################################
@ -72,11 +73,11 @@ knncdf = csiborgtools.match.kNN_CDF()
def do_auto(ic):
out = {}
cat = csiborgtools.read.HaloCatalogue(ic, max_dist=Rmax)
cat = csiborgtools.read.HaloCatalogue(ic, paths, max_dist=Rmax)
for i, mmin in enumerate(mass_threshold):
knn = NearestNeighbors()
knn.fit(cat.positions[cat["totpartmass"] > mmin, ...])
knn.fit(cat.positions(False)[cat["totpartmass"] > mmin, ...])
rs, cdf = knncdf(knn, nneighbours=args.nneighbours, Rmax=Rmax,
rmin=args.rmin, rmax=args.rmax, nsamples=args.nsamples,
@ -90,15 +91,15 @@ def do_auto(ic):
def do_cross(ics):
out = {}
cat1 = csiborgtools.read.HaloCatalogue(ics[0], max_dist=Rmax)
cat2 = csiborgtools.read.HaloCatalogue(ics[1], max_dist=Rmax)
cat1 = csiborgtools.read.HaloCatalogue(ics[0], paths, max_dist=Rmax)
cat2 = csiborgtools.read.HaloCatalogue(ics[1], paths, max_dist=Rmax)
for i, mmin in enumerate(mass_threshold):
knn1 = NearestNeighbors()
knn1.fit(cat1.positions[cat1["totpartmass"] > mmin, ...])
knn1.fit(cat1.positions()[cat1["totpartmass"] > mmin, ...])
knn2 = NearestNeighbors()
knn2.fit(cat2.positions[cat2["totpartmass"] > mmin, ...])
knn2.fit(cat2.positions()[cat2["totpartmass"] > mmin, ...])
rs, cdf0, cdf1, joint_cdf = knncdf.joint(
knn1, knn2, nneighbours=args.nneighbours, Rmax=Rmax,

9
scripts/run_singlematch.py
View file

@ -35,6 +35,7 @@ parser.add_argument("--sigma", type=float)
args = parser.parse_args()
# File paths
paths = csiborgtools.read.CSiBORGPaths()
fout = join(utils.dumpdir, "overlap",
"cross_{}_{}.npz".format(args.nsim0, args.nsimx))
smooth_kwargs = {"sigma": args.sigma, "mode": "constant", "cval": 0.0}
@ -43,18 +44,18 @@ overlapper = csiborgtools.match.ParticleOverlap()
# Load catalogues
print("{}: loading catalogues {} and {}."
.format(datetime.now(), args.nsim0, args.nsimx), flush=True)
cat0 = csiborgtools.read.HaloCatalogue(args.nsim0)
catx = csiborgtools.read.HaloCatalogue(args.nsimx)
cat0 = csiborgtools.read.HaloCatalogue(args.nsim0, paths)
catx = csiborgtools.read.HaloCatalogue(args.nsimx, paths)
print("{}: loading simulation {} and converting positions to cell numbers."
.format(datetime.now(), args.nsim0), flush=True)
with open(cat0.paths.clump0_path(args.nsim0), "rb") as f:
with open(paths.clump0_path(args.nsim0), "rb") as f:
clumps0 = numpy.load(f, allow_pickle=True)
overlapper.clumps_pos2cell(clumps0)
print("{}: loading simulation {} and converting positions to cell numbers."
.format(datetime.now(), args.nsimx), flush=True)
with open(catx.paths.clump0_path(args.nsimx), 'rb') as f:
with open(paths.clump0_path(args.nsimx), 'rb') as f:
clumpsx = numpy.load(f, allow_pickle=True)
overlapper.clumps_pos2cell(clumpsx)

28
scripts/run_split_halos.py
View file

@ -15,9 +15,8 @@
"""
Script to split particles into smaller files according to their clump
membership for faster manipulation. Currently does this for the maximum
snapshot of each simulation. Running this will require a lot of memory.
snapshot of each simulation. Running this requires a lot of memory.
"""
from mpi4py import MPI
from datetime import datetime
try:
@ -34,27 +33,26 @@ rank = comm.Get_rank()
nproc = comm.Get_size()
paths = csiborgtools.read.CSiBORGPaths()
n_sims = paths.ic_ids[:1]
sims = paths.ic_ids(False)
partcols = ["x", "y", "z", "vx", "vy", "vz", "M", "level"]
jobs = csiborgtools.fits.split_jobs(len(n_sims), nproc)[rank]
jobs = csiborgtools.fits.split_jobs(len(sims), nproc)[rank]
for icount, sim_index in enumerate(jobs):
print("{}: rank {} working {} / {} jobs.".format(datetime.now(), rank,
icount + 1, len(jobs)))
n_sim = n_sims[sim_index]
n_snap = paths.get_maximum_snapshot(n_sim)
# Set paths and inifitalise a particle reader
paths.set_info(n_sim, n_snap)
print("{}: rank {} working {} / {} jobs."
.format(datetime.now(), rank, icount + 1, len(jobs)), flush=True)
nsim = sims[sim_index]
nsnap = max(paths.get_snapshots(nsim))
partreader = csiborgtools.read.ParticleReader(paths)
# Load the clumps, particles' clump IDs and particles.
clumps = partreader.read_clumps()
particle_clumps = partreader.read_clumpid(verbose=False)
particles = partreader.read_particle(partcols, verbose=False)
clumps = partreader.read_clumps(nsnap, nsim)
particle_clumps = partreader.read_clumpid(nsnap, nsim, verbose=False)
particles = partreader.read_particle(nsnap, nsim, partcols, verbose=False)
# Drop all particles whose clump index is 0 (not assigned to any halo)
particle_clumps, particles = partreader.drop_zero_indx(
particle_clumps, particles)
# Dump it!
csiborgtools.fits.dump_split_particles(particles, particle_clumps, clumps,
utils.Nsplits, paths, verbose=False)
utils.Nsplits, nsnap, nsim, paths,
verbose=False)
print("All finished!")
print("All finished!", flush=True)