From d578c71b83b34d49425f9d0a057bf3fb2ac01479 Mon Sep 17 00:00:00 2001 From: Richard Stiskalek Date: Wed, 14 Aug 2024 13:02:38 +0200 Subject: [PATCH] LSS projected basics (#140) * Move files * Move files * Add galactic to RA/dec * Update sky maps * Add projected fields * Remove old import * Quick update * Add IO * Add imports * Update imports * Add basic file --- csiborgtools/__init__.py | 2 +- csiborgtools/cmb_correlation/__init__.py | 16 ++ csiborgtools/cmb_correlation/io.py | 59 +++++++ csiborgtools/field/__init__.py | 4 +- csiborgtools/field/interp.py | 96 +++-------- csiborgtools/read/paths.py | 21 +++ csiborgtools/summary/__init__.py | 1 + csiborgtools/summary/cmb_correlation.py | 80 +++++++++ csiborgtools/utils.py | 6 + scripts/{ => field_prop}/field_bulk.py | 0 scripts/{ => field_prop}/field_bulk.sh | 0 scripts/{ => field_prop}/field_los.py | 0 scripts/{ => field_prop}/field_los.sh | 0 scripts/field_prop/field_projected.py | 127 ++++++++++++++ scripts/field_prop/field_projected.sh | 26 +++ scripts/{ => field_prop}/field_prop.py | 0 scripts/{ => field_prop}/field_prop.sh | 0 scripts/{ => field_prop}/field_sample.py | 0 scripts/{ => field_prop}/field_sample.sh | 0 scripts/{ => flow}/flow_validation.py | 0 scripts/{ => flow}/flow_validation.sh | 0 scripts/{ => flow}/post_upglade.py | 0 scripts/{ => flow}/post_upglade.sh | 0 scripts/{ => flow}/quijote_bulkflow.py | 0 scripts/{ => flow}/quijote_bulkflow.sh | 0 scripts/{ => flow}/quijote_pecvel_covmat.py | 0 scripts/{ => flow}/quijote_pecvel_covmat.sh | 0 scripts/{ => halo_overlap}/fit_init.py | 0 scripts/{ => halo_overlap}/fit_init.sh | 0 .../{ => halo_overlap}/match_overlap_all.py | 0 .../{ => halo_overlap}/match_overlap_all.sh | 0 .../match_overlap_single.py | 0 .../match_overlap_single.sh | 0 scripts/{ => mah}/mah_random.py | 0 scripts/{ => mah}/mah_random.sh | 0 scripts/vrad_covariance.py | 158 ------------------ 36 files changed, 365 insertions(+), 231 deletions(-) create mode 100644 csiborgtools/cmb_correlation/__init__.py create mode 100644 csiborgtools/cmb_correlation/io.py create mode 100644 csiborgtools/summary/cmb_correlation.py rename scripts/{ => field_prop}/field_bulk.py (100%) rename scripts/{ => field_prop}/field_bulk.sh (100%) rename scripts/{ => field_prop}/field_los.py (100%) rename scripts/{ => field_prop}/field_los.sh (100%) create mode 100644 scripts/field_prop/field_projected.py create mode 100755 scripts/field_prop/field_projected.sh rename scripts/{ => field_prop}/field_prop.py (100%) rename scripts/{ => field_prop}/field_prop.sh (100%) rename scripts/{ => field_prop}/field_sample.py (100%) rename scripts/{ => field_prop}/field_sample.sh (100%) rename scripts/{ => flow}/flow_validation.py (100%) rename scripts/{ => flow}/flow_validation.sh (100%) rename scripts/{ => flow}/post_upglade.py (100%) rename scripts/{ => flow}/post_upglade.sh (100%) rename scripts/{ => flow}/quijote_bulkflow.py (100%) rename scripts/{ => flow}/quijote_bulkflow.sh (100%) rename scripts/{ => flow}/quijote_pecvel_covmat.py (100%) rename scripts/{ => flow}/quijote_pecvel_covmat.sh (100%) rename scripts/{ => halo_overlap}/fit_init.py (100%) rename scripts/{ => halo_overlap}/fit_init.sh (100%) rename scripts/{ => halo_overlap}/match_overlap_all.py (100%) rename scripts/{ => halo_overlap}/match_overlap_all.sh (100%) rename scripts/{ => halo_overlap}/match_overlap_single.py (100%) rename scripts/{ => halo_overlap}/match_overlap_single.sh (100%) rename scripts/{ => mah}/mah_random.py (100%) rename scripts/{ => mah}/mah_random.sh (100%) delete mode 100644 scripts/vrad_covariance.py diff --git a/csiborgtools/__init__.py b/csiborgtools/__init__.py index 2c9fa70..e4b7afc 100644 --- a/csiborgtools/__init__.py +++ b/csiborgtools/__init__.py @@ -22,7 +22,7 @@ from .utils import (center_of_mass, delta2ncells, number_counts, radec_to_cartesian, cartesian_to_radec, # noqa thin_samples_by_acl, BIC_AIC, radec_to_galactic, # noqa heliocentric_to_cmb, calculate_acl, harmonic_evidence, # noqa - dict_samples_to_array) # noqa + dict_samples_to_array, galactic_to_radec) # noqa from .params import (paths_glamdring, simname2boxsize, simname2Omega_m, # noqa snap2redshift) # noqa diff --git a/csiborgtools/cmb_correlation/__init__.py b/csiborgtools/cmb_correlation/__init__.py new file mode 100644 index 0000000..8e56a48 --- /dev/null +++ b/csiborgtools/cmb_correlation/__init__.py @@ -0,0 +1,16 @@ +# Copyright (C) 2024 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. + +from .io import read_CMB_temperature # noqa \ No newline at end of file diff --git a/csiborgtools/cmb_correlation/io.py b/csiborgtools/cmb_correlation/io.py new file mode 100644 index 0000000..1b854f4 --- /dev/null +++ b/csiborgtools/cmb_correlation/io.py @@ -0,0 +1,59 @@ +# 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. +""" +Various I/O functions for reading and writing data. +""" +from warnings import warn + +import healpy as hp +from astropy.io import fits + + +def read_CMB_temperature(fname=None, nside_out=None, + convert_to_ring_ordering=True, normalize=False, + verbose=True): + """ + Read the CMB temperature map from a FITS file. + """ + if fname is None: + warn("Using the glamdrnig path to the default temperature map.", + UserWarning) + fname = "/mnt/extraspace/rstiskalek/catalogs/CMB/COM_CMB_IQU-smica_2048_R3.00_full.fits" # noqa + + f = fits.open(fname) + if verbose: + print(f"Reading CMB temperature map from `{fname}`.") + + skymap = f[1].data["I_STOKES"] + mask = f[1].data["TMASK"] + f.close() + + if nside_out is not None: + if verbose: + print(f"Moving to nside = {nside_out}...") + skymap = hp.pixelfunc.ud_grade(skymap, nside_out, order_in="NESTED") + mask = hp.pixelfunc.ud_grade(mask, nside_out, order_in="NESTED") + + if convert_to_ring_ordering: + if verbose: + print("Converting to RING ordering...") + skymap = hp.reorder(skymap, n2r=True) + mask = hp.reorder(mask, n2r=True) + + if normalize: + skymap -= skymap.mean() + skymap /= skymap.std() + + return skymap, mask diff --git a/csiborgtools/field/__init__.py b/csiborgtools/field/__init__.py index 53e41a5..1570096 100644 --- a/csiborgtools/field/__init__.py +++ b/csiborgtools/field/__init__.py @@ -19,8 +19,8 @@ try: from .density import (DensityField, PotentialField, TidalTensorField, # noqa VelocityField, radial_velocity, power_spectrum, # noqa overdensity_field) # noqa - from .interp import (evaluate_cartesian_cic, evaluate_sky, evaluate_los, # noqa - field2rsp, fill_outside, make_sky, # noqa + from .interp import (evaluate_cartesian_cic, evaluate_los, field2rsp, # noqa + fill_outside, make_sky, # noqa observer_peculiar_velocity, smoothen_field, # noqa field_at_distance) # noqa except ImportError: diff --git a/csiborgtools/field/interp.py b/csiborgtools/field/interp.py index 3311718..e27bea5 100644 --- a/csiborgtools/field/interp.py +++ b/csiborgtools/field/interp.py @@ -20,7 +20,7 @@ import numpy as np import smoothing_library as SL from numba import jit from scipy.interpolate import RegularGridInterpolator -from tqdm import tqdm, trange +from tqdm import tqdm from ..utils import periodic_wrap_grid, radec_to_cartesian from .utils import divide_nonzero, force_single_precision, nside2radec @@ -303,56 +303,13 @@ def evaluate_los(*fields, sky_pos, boxsize, rmax, dr, smooth_scales=None, ############################################################################### -def evaluate_sky(*fields, pos, mpc2box, smooth_scales=None, verbose=False): - """ - Evaluate a scalar field(s) at radial distance `Mpc / h`, right ascensions - [0, 360) deg and declinations [-90, 90] deg. Uses CIC interpolation. - - Parameters - ---------- - fields : (list of) 3-dimensional array of shape `(grid, grid, grid)` - Field to be interpolated. - pos : 2-dimensional array of shape `(n_samples, 3)` - Query spherical coordinates. - mpc2box : float - Conversion factor to multiply the radial distance by to get box units. - smooth_scales : (list of) float, optional - Smoothing scales in `Mpc / h`. If `None`, no smoothing is performed. - verbose : bool, optional - Smoothing verbosity flag. - - Returns - ------- - (list of) 1-dimensional array of shape `(n_samples, len(smooth_scales))` - """ - # Make a copy of the positions to avoid modifying the input. - pos = np.copy(pos) - - pos = force_single_precision(pos) - pos[:, 0] *= mpc2box - - cart_pos = radec_to_cartesian(pos) + 0.5 - - if smooth_scales is not None: - if isinstance(smooth_scales, (int, float)): - smooth_scales = [smooth_scales] - - if isinstance(smooth_scales, list): - smooth_scales = np.array(smooth_scales, dtype=np.float32) - - smooth_scales *= mpc2box - - return evaluate_cartesian_cic(*fields, pos=cart_pos, - smooth_scales=smooth_scales, - verbose=verbose) - - -def make_sky(field, angpos, dist, boxsize, verbose=True): +def make_sky(field, angpos, rmax, dr, boxsize, interpolation_method="cic", + return_full=False, verbose=True): r""" Make a sky map of a scalar field. The observer is in the centre of the box the field is evaluated along directions `angpos` (RA [0, 360) deg, - dec [-90, 90] deg). Along each direction, the field is evaluated distances - `dist` (`Mpc / h`) and summed. Uses CIC interpolation. + dec [-90, 90] deg). The field is evaluated up to `rmax` with a linear + spacing of `dr` in `Mpc / h`. Parameters ---------- @@ -360,39 +317,38 @@ def make_sky(field, angpos, dist, boxsize, verbose=True): Field to be interpolated angpos : 2-dimensional arrays of shape `(ndir, 2)` Directions to evaluate the field. - dist : 1-dimensional array - Uniformly spaced radial distances to evaluate the field in `Mpc / h`. + rmax : float + Maximum radial distance in `Mpc / h`. + dr : float + Radial distance step in `Mpc / h`. boxsize : float Box size in `Mpc / h`. + interpolation_method : str, optional + Interpolation method. Must be one of `cic` or one of the methods of + `scipy.interpolate.RegularGridInterpolator`. + return_full : bool, optional + If `True`, return the full interpolated field instead of the average + field at each radial distance. verbose : bool, optional Verbosity flag. Returns ------- - interp_field : 1-dimensional array of shape `(n_pos, )`. + interp_field : 1-dimensional array of shape `(n_pos, )` """ - dx = dist[1] - dist[0] - assert np.allclose(dist[1:] - dist[:-1], dx) - assert angpos.ndim == 2 and dist.ndim == 1 + rdist, finterp = evaluate_los( + field, sky_pos=angpos, boxsize=boxsize, rmax=rmax, dr=dr, + smooth_scales=None, verbose=verbose, + interpolation_method=interpolation_method) - # We loop over the angular directions, at each step evaluating a vector - # of distances. We pre-allocate arrays for speed. - dir_loop = np.full((dist.size, 3), np.nan, dtype=np.float32) + if return_full: + return rdist, finterp - ndir = angpos.shape[0] - out = np.full(ndir, np.nan, dtype=np.float32) - for i in trange(ndir) if verbose else range(ndir): - dir_loop[:, 0] = dist - dir_loop[:, 1] = angpos[i, 0] - dir_loop[:, 2] = angpos[i, 1] + finterp *= rdist**2 + finterp = np.trapz(finterp, x=rdist, axis=1) + finterp /= np.trapz(rdist**2, x=rdist) - out[i] = np.sum( - dist**2 * evaluate_sky(field, pos=dir_loop, mpc2box=1 / boxsize)) - - # Assuming the field is in h^2 Msun / kpc**3, we need to convert Mpc / h - # to kpc / h and multiply by the pixel area. - out *= dx * 1e9 * 4 * np.pi / len(angpos) - return out + return finterp ############################################################################### diff --git a/csiborgtools/read/paths.py b/csiborgtools/read/paths.py index 0cdb331..a66c1a3 100644 --- a/csiborgtools/read/paths.py +++ b/csiborgtools/read/paths.py @@ -634,3 +634,24 @@ class Paths: fdir = join(self.postdir, "field_los") try_create_directory(fdir) return join(fdir, f"los_{catalogue}_{simnname}.hdf5") + + def field_projected(self, simname, kind): + """ + Path to the files containing the projected fields on the sky. + + Parameters + ---------- + simname : str + Simulation name. + kind : str + Field type. + + Returns + ------- + str + """ + fdir = join(self.postdir, "field_projected") + try_create_directory(fdir) + return join(fdir, f"{simname}_{kind}_volume_weighted.hdf5") + + diff --git a/csiborgtools/summary/__init__.py b/csiborgtools/summary/__init__.py index a08e165..3c85d95 100644 --- a/csiborgtools/summary/__init__.py +++ b/csiborgtools/summary/__init__.py @@ -13,6 +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 .cmb_correlation import read_projected_matter # noqa from .knn_summary import kNNCDFReader # noqa from .nearest_neighbour_summary import NearestNeighbourReader # noqa from .overlap_summary import weighted_stats # noqa diff --git a/csiborgtools/summary/cmb_correlation.py b/csiborgtools/summary/cmb_correlation.py new file mode 100644 index 0000000..4c2f8e6 --- /dev/null +++ b/csiborgtools/summary/cmb_correlation.py @@ -0,0 +1,80 @@ +# Copyright (C) 2024 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. +""" +Useful functions for getting summary data for CMB x MATTER cross-correlation. +""" +import numpy as np +from h5py import File +from healpy.sphtfunc import smoothing +import healpy as hp +from tqdm import tqdm + + +def read_projected_matter(simname, paths, fwhm_deg=None, remove_monopole=False, + remove_dipole=False, normalize=False): + """ + Read the projected matter density field for a given simulation. + + Parameters + ---------- + simname : str + The name of the simulation. + paths : csiborgtools.read.Paths + Paths object. + fwhm_deg : float, optional + The full-width at half-maximum of the smoothing kernel in degrees. + remove_monopole : bool, optional + Whether to remove the monopole from the field. + remove_dipole : bool, optional + Whether to remove the dipole from the field. + normalize : bool, optional + Whether to apply standard normalization to the field. + + Returns + ------- + dist_ranges : 2-dimensional array of shape (n_dist_ranges, 2) + The distance ranges for the field. + data : 3-dimensional array of shape (n_sims, n_dist_ranges, npix) + The projected matter density field. + """ + kind = "density" + nsims = paths.get_ics(simname) + + fname = paths.field_projected(simname, kind) + with File(fname, "r") as f: + dist_ranges = f["dist_ranges"][...] + + npix = len(f[f"nsim_{nsims[0]}/dist_range_0"]) + + data = np.zeros((len(nsims), len(dist_ranges), npix)) + for i, nsim in enumerate(tqdm(nsims, desc="Simulations")): + for j in range(len(dist_ranges)): + skymap = f[f"nsim_{nsim}/dist_range_{j}"][...] + + if fwhm_deg is not None: + skymap = smoothing(skymap, fwhm=fwhm_deg * np.pi / 180.0) + + if remove_monopole: + hp.pixelfunc.remove_monopole(skymap, copy=False) + if remove_dipole: + hp.pixelfunc.remove_dipole(skymap, copy=False) + + if normalize: + skymap -= np.mean(skymap) + skymap /= np.std(skymap) + + data[i, j] = skymap + + return dist_ranges, data diff --git a/csiborgtools/utils.py b/csiborgtools/utils.py index 85d2492..c986457 100644 --- a/csiborgtools/utils.py +++ b/csiborgtools/utils.py @@ -163,6 +163,12 @@ def radec_to_galactic(ra, dec): return c.galactic.l.degree, c.galactic.b.degree +def galactic_to_radec(l, b): # noqa + """Convert galactic coordinates to right ascension and declination.""" + c = SkyCoord(l=l*u.degree, b=b*u.degree, frame='galactic') + return c.icrs.ra.degree, c.icrs.dec.degree + + def radec_to_supergalactic(ra, dec): """Convert right ascension and declination to supergalactic coordinates.""" c = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs') diff --git a/scripts/field_bulk.py b/scripts/field_prop/field_bulk.py similarity index 100% rename from scripts/field_bulk.py rename to scripts/field_prop/field_bulk.py diff --git a/scripts/field_bulk.sh b/scripts/field_prop/field_bulk.sh similarity index 100% rename from scripts/field_bulk.sh rename to scripts/field_prop/field_bulk.sh diff --git a/scripts/field_los.py b/scripts/field_prop/field_los.py similarity index 100% rename from scripts/field_los.py rename to scripts/field_prop/field_los.py diff --git a/scripts/field_los.sh b/scripts/field_prop/field_los.sh similarity index 100% rename from scripts/field_los.sh rename to scripts/field_prop/field_los.sh diff --git a/scripts/field_prop/field_projected.py b/scripts/field_prop/field_projected.py new file mode 100644 index 0000000..72731a7 --- /dev/null +++ b/scripts/field_prop/field_projected.py @@ -0,0 +1,127 @@ +# Copyright (C) 2023 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 calculate the projected density field for a given simulation as a +sky map. The script is not parallelized in any way. The generated fields are +converted to galactic coordinates to match the CMB maps. +""" +from argparse import ArgumentParser +from os import remove +from os.path import exists, join + +import csiborgtools +import numpy as np +from h5py import File + + +def get_field(simname, nsim, field_kind, MAS, grid): + """Get the appropriate field reader for the simulation.""" + if simname == "csiborg1": + reader = csiborgtools.read.CSiBORG1Field(nsim) + elif "csiborg2" in simname: + kind = simname.split("_")[-1] + reader = csiborgtools.read.CSiBORG2Field(nsim, kind) + elif simname == "csiborg2X": + reader = csiborgtools.read.CSiBORG2XField(nsim) + elif "quijote" in simname: + reader = csiborgtools.read.QuijoteField(nsim) + else: + raise ValueError(f"Unknown simname: `{simname}`.") + + if field_kind == "density": + return reader.density_field(MAS, grid) + else: + raise ValueError(f"Unknown field kind: `{field_kind}`.") + + +def main(simname, nsims, field_kind, nside, dist_ranges, MAS, grid, + volume_weight, folder, normalize_to_overdensity=True): + boxsize = csiborgtools.simname2boxsize(simname) + Om0 = csiborgtools.simname2Omega_m(simname) + matter_density = Om0 * 277.53662724583074 # Msun / kpc^3 + + fname = join(folder, f"{simname}_{field_kind}.hdf5") + if volume_weight: + fname = fname.replace(".hdf5", "_volume_weighted.hdf5") + print(f"Writing to `{fname}`...") + if exists(fname): + remove(fname) + + with File(fname, "w") as f: + f.create_dataset("dist_ranges", data=np.asarray(dist_ranges)) + f.create_dataset("nsims", data=nsims) + + # These are at first generated in RA/dec but we can assume it is galactic + # and convert it to RA/dec. + pixel_angpos = csiborgtools.field.nside2radec(nside) + RA, dec = csiborgtools.galactic_to_radec(*pixel_angpos.T) + pixel_angpos = np.vstack([RA, dec]).T + npix = len(pixel_angpos) + + Rmax = np.asanyarray(dist_ranges).reshape(-1).max() + dr = 0.5 * boxsize / grid + print(f"{'R_max:':<20} {Rmax} Mpc / h", flush=True) + print(f"{'dr:':<20} {dr} Mpc / h", flush=True) + + for nsim in nsims: + print(f"Interpolating at {npix} pixel for simulation {nsim}...", + flush=True) + + field = get_field(simname, nsim, field_kind, MAS, grid) + rdist, finterp = csiborgtools.field.make_sky( + field, pixel_angpos, Rmax, dr, boxsize, return_full=True, + interpolation_method="linear") + + with File(fname, "a") as f: + grp = f.create_group(f"nsim_{nsim}") + + for n in range(len(dist_ranges)): + dmin, dmax = dist_ranges[n] + k_start = np.searchsorted(rdist, dmin) + k_end = np.searchsorted(rdist, dmax) + + r = rdist[k_start:k_end + 1] + y = r**2 * finterp[:, k_start:k_end + 1] + skymap = np.trapz(y, r, axis=-1) / np.trapz(r**2, r) + + if normalize_to_overdensity: + skymap /= matter_density + skymap -= 1 + + grp.create_dataset(f"dist_range_{n}", data=skymap) + + +if __name__ == "__main__": + parser = ArgumentParser() + parser.add_argument("--simname", type=str, help="Simulation name.") + args = parser.parse_args() + + fdir = "/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_projected" + dx = 25 + dist_ranges = [[0, n * dx] for n in range(1, 5)] + dist_ranges += [[n * dx, (n + 1) * dx] for n in range(0, 5)] + + paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring) + nsims = paths.get_ics(args.simname) + print(f"{'Num. sims:':<20} {len(nsims)}", flush=True) + + MAS = "SPH" + grid = 1024 + nside = 128 + field_kind = "density" + volume_weight = True + + main(args.simname, nsims, field_kind, nside, dist_ranges, MAS, grid, + volume_weight, fdir) diff --git a/scripts/field_prop/field_projected.sh b/scripts/field_prop/field_projected.sh new file mode 100755 index 0000000..1192e00 --- /dev/null +++ b/scripts/field_prop/field_projected.sh @@ -0,0 +1,26 @@ +nthreads=1 +memory=64 +on_login=${1} +simname=${2} +queue="berg" +env="/mnt/users/rstiskalek/csiborgtools/venv_csiborg/bin/python" +file="field_projected.py" + + +if [ "$on_login" != "1" ] && [ "$on_login" != "0" ] +then + echo "'on_login' (1) must be either 0 or 1." + exit 1 +fi + +pythoncm="$env $file --simname $simname" +if [ $on_login -eq 1 ]; then + echo $pythoncm + $pythoncm +else + cm="addqueue -q $queue -n $nthreads -m $memory $pythoncm" + echo "Submitting:" + echo $cm + echo + eval $cm +fi diff --git a/scripts/field_prop.py b/scripts/field_prop/field_prop.py similarity index 100% rename from scripts/field_prop.py rename to scripts/field_prop/field_prop.py diff --git a/scripts/field_prop.sh b/scripts/field_prop/field_prop.sh similarity index 100% rename from scripts/field_prop.sh rename to scripts/field_prop/field_prop.sh diff --git a/scripts/field_sample.py b/scripts/field_prop/field_sample.py similarity index 100% rename from scripts/field_sample.py rename to scripts/field_prop/field_sample.py diff --git a/scripts/field_sample.sh b/scripts/field_prop/field_sample.sh similarity index 100% rename from scripts/field_sample.sh rename to scripts/field_prop/field_sample.sh diff --git a/scripts/flow_validation.py b/scripts/flow/flow_validation.py similarity index 100% rename from scripts/flow_validation.py rename to scripts/flow/flow_validation.py diff --git a/scripts/flow_validation.sh b/scripts/flow/flow_validation.sh similarity index 100% rename from scripts/flow_validation.sh rename to scripts/flow/flow_validation.sh diff --git a/scripts/post_upglade.py b/scripts/flow/post_upglade.py similarity index 100% rename from scripts/post_upglade.py rename to scripts/flow/post_upglade.py diff --git a/scripts/post_upglade.sh b/scripts/flow/post_upglade.sh similarity index 100% rename from scripts/post_upglade.sh rename to scripts/flow/post_upglade.sh diff --git a/scripts/quijote_bulkflow.py b/scripts/flow/quijote_bulkflow.py similarity index 100% rename from scripts/quijote_bulkflow.py rename to scripts/flow/quijote_bulkflow.py diff --git a/scripts/quijote_bulkflow.sh b/scripts/flow/quijote_bulkflow.sh similarity index 100% rename from scripts/quijote_bulkflow.sh rename to scripts/flow/quijote_bulkflow.sh diff --git a/scripts/quijote_pecvel_covmat.py b/scripts/flow/quijote_pecvel_covmat.py similarity index 100% rename from scripts/quijote_pecvel_covmat.py rename to scripts/flow/quijote_pecvel_covmat.py diff --git a/scripts/quijote_pecvel_covmat.sh b/scripts/flow/quijote_pecvel_covmat.sh similarity index 100% rename from scripts/quijote_pecvel_covmat.sh rename to scripts/flow/quijote_pecvel_covmat.sh diff --git a/scripts/fit_init.py b/scripts/halo_overlap/fit_init.py similarity index 100% rename from scripts/fit_init.py rename to scripts/halo_overlap/fit_init.py diff --git a/scripts/fit_init.sh b/scripts/halo_overlap/fit_init.sh similarity index 100% rename from scripts/fit_init.sh rename to scripts/halo_overlap/fit_init.sh diff --git a/scripts/match_overlap_all.py b/scripts/halo_overlap/match_overlap_all.py similarity index 100% rename from scripts/match_overlap_all.py rename to scripts/halo_overlap/match_overlap_all.py diff --git a/scripts/match_overlap_all.sh b/scripts/halo_overlap/match_overlap_all.sh similarity index 100% rename from scripts/match_overlap_all.sh rename to scripts/halo_overlap/match_overlap_all.sh diff --git a/scripts/match_overlap_single.py b/scripts/halo_overlap/match_overlap_single.py similarity index 100% rename from scripts/match_overlap_single.py rename to scripts/halo_overlap/match_overlap_single.py diff --git a/scripts/match_overlap_single.sh b/scripts/halo_overlap/match_overlap_single.sh similarity index 100% rename from scripts/match_overlap_single.sh rename to scripts/halo_overlap/match_overlap_single.sh diff --git a/scripts/mah_random.py b/scripts/mah/mah_random.py similarity index 100% rename from scripts/mah_random.py rename to scripts/mah/mah_random.py diff --git a/scripts/mah_random.sh b/scripts/mah/mah_random.sh similarity index 100% rename from scripts/mah_random.sh rename to scripts/mah/mah_random.sh diff --git a/scripts/vrad_covariance.py b/scripts/vrad_covariance.py deleted file mode 100644 index c8ad263..0000000 --- a/scripts/vrad_covariance.py +++ /dev/null @@ -1,158 +0,0 @@ -# Copyright (C) 2024 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. -""" -NOTE: The script is far from finished or written well. -""" -from os.path import join -import numpy as np - -import csiborgtools -from h5py import File - - - -from sklearn.neighbors import NearestNeighbors -from numba import jit -from scipy.stats import binned_statistic - - -### - -def find_indxs(rtree, r, radius): - """ - Find the indices of points that are within a given radius of a given - point `r`. - """ - if isinstance(r, (int, float)): - r = np.array(r) - - return rtree.radius_neighbors( - r.reshape(-1, 1), radius=radius, return_distance=False, )[0] - -@jit(nopython=True) -def dot_product_norm(x1_norm, x2_norm): - """Dot product of two normalised 1D vectors.""" - return x1_norm[0] * x2_norm[0] + x1_norm[1] * x2_norm[1] + x1_norm[2] * x2_norm[2] # noqa - - - -@jit(nopython=True) -def get_angdist_vrad_product(i_comb, j_comb, pos_norm, vrad, r, rbin_i, rbin_j): - # TODO: Add itself? - len_i = len(i_comb) - len_j = len(j_comb) - - cos_angdist_values = np.full(len_i * len_j, np.nan) - vrad_product = np.full(len_i * len_j, np.nan) - weights = np.full(len_i * len_j, np.nan) - - k = 0 - for i in range(len_i): - pos_norm_i = pos_norm[i_comb[i]] - vrad_i = vrad[i_comb[i]] - w_i = (rbin_i / r[i_comb[i]])**2 - for j in range(len_j): - cos_angdist_values[k] = dot_product_norm(pos_norm_i, pos_norm[j_comb[j]]) - - # Product of the peculiar velocities - vrad_product[k] = vrad_i * vrad[j_comb[j]] - # Weight the product - w = w_i * (rbin_j / r[j_comb[j]])**2 - vrad_product[k] *= w - - weights[k] = w - - k += 1 - - return cos_angdist_values, vrad_product, weights - - -def main(out_summed_product, out_weights, ri, rj, angular_bins, pos, vel, observer, rmax): - # Centre the positions at the observer - pos = np.copy(pos) - observer - r = np.linalg.norm(pos, axis=1) - mask = r < rmax - - # Select only the haloes within the radial range - pos = pos[mask] - vel = vel[mask] - r = r[mask] - - # Create a KDTree for the radial positions - rtree = NearestNeighbors().fit(r.reshape(-1, 1)) - - # Calculate the radial velocity and the normalised position vector - pos_norm = pos / r[:, None] - vrad = np.sum(vel * pos_norm, axis=1) - - # TODO: eventually here loop over the radii - # for .... - dr = 2.5 - i_indxs = find_indxs(rtree, ri, radius=dr) - j_indxs = find_indxs(rtree, rj, radius=dr) - - # Calculate the cosine of the angular distance and the product of the - # radial velocities for each pair of points. - cos_angdist, vrad_product, weights = get_angdist_vrad_product( - i_indxs, j_indxs, pos_norm, vrad, r, ri, rj) - - out_summed_product += binned_statistic( - cos_angdist, vrad_product, bins=angular_bins, statistic="sum")[0] - - out_weights += binned_statistic( - cos_angdist, weights, bins=angular_bins, statistic="sum")[0] - - return out_summed_product, out_weights - - -if __name__ == "__main__": - fdir = "/mnt/extraspace/rstiskalek/BBF/Quijote_C_ij" - - rmax = 150 - nradial = 20 - nangular = 40 - ncatalogue = 1 - - # NOTE check this - radial_bins = np.linspace(0, rmax, nradial + 1) - angular_bins = np.linspace(-1, 1, nangular + 1) - - summed_product = np.zeros(nangular) - weights = np.zeros(nangular) - - fiducial_observers = csiborgtools.read.fiducial_observers(1000, rmax) - - ri = 100 - rj = 120 - - # for i in trange(ncatalogue, desc="Catalogues"): - for i in [30]: - cat = csiborgtools.read.QuijoteCatalogue(i) - - for j in range(len(fiducial_observers)): - # Loop over all the fiducial observers in this simulation - observer = fiducial_observers[j] - summed_product, weights = main( - summed_product, weights, ri, rj, angular_bins, - cat["cartesian_pos"], cat["cartesian_vel"], observer, rmax) - - - # TODO: Save - fname = join(fdir, "test.h5") - print(f"Saving to `{fname}`.") - with File(fname, 'w') as f: - f.create_dataset("summed_product", data=summed_product) - f.create_dataset("weights", data=weights) - f.create_dataset("angular_bins", data=angular_bins)