More void (#144)

* Add more void

* Update script

* Add check

* Add loop over LG observers

* Remove never used code

* Update docs

csiborgtools/flow/flow_model.py

@@ -22,7 +22,6 @@ References
 [1] https://arxiv.org/abs/1912.09383.
 """
 from abc import ABC, abstractmethod
-from os.path import join
 
 import numpy as np
 from astropy import units as u
@@ -202,14 +201,6 @@ class DataLoader:
             fpath = paths.field_los(simname, "Pantheon+")
         elif "CF4_TFR" in catalogue:
             fpath = paths.field_los(simname, "CF4_TFR")
-        elif "IndranilVoid" in catalogue:
-            fdir = "/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_los"  # noqa
-            if "exp" in catalogue:
-                fpath = join(fdir, "v_pec_EXP_IndranilVoid.dat")
-            elif "gauss" in catalogue:
-                fpath = join(fdir, "v_pec_GAUSS_IndranilVoid.dat")
-            else:
-                raise ValueError("Unknown `IndranilVoid` catalogue.")
         else:
             fpath = paths.field_los(simname, catalogue)
 

csiborgtools/read/paths.py

@@ -125,8 +125,21 @@ class Paths:
             files = [int(file) for file in files]
             # Downsample to only 20 realisations
             files = files[::5]
-        elif simname in ["Carrick2015", "Lilow2024", "no_field"] or "IndranilVoid" in simname:  # noqa
+        elif simname in ["Carrick2015", "Lilow2024", "no_field"]:
             files = [0]
+        elif "IndranilVoid" in simname:
+            kind = simname.split("_")[-1]
+            if kind not in ["exp", "gauss"]:
+                raise ValueError(f"Unknown void kind `{simname}`.")
+
+            kind = kind.upper()
+            fdir = join(self.aux_cat_dir, "IndranilVoid", f"{kind}profile")
+            files = glob(join(fdir, "v_pec_*.dat"))
+
+            files = [
+                search(rf'v_pec_{kind}profile_rLG_(\d+)\.dat', file).group(1)
+                for file in files]
+            files = [int(file) for file in files]
         else:
             raise ValueError(f"Unknown simulation name `{simname}`.")
 

scripts/field_prop/field_los_indranil_void.py

@@ -11,10 +11,7 @@
 # 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.
-"""
+"""Script to interpolate the Indranil void profiles for lines of sight."""
-MPI script to interpolate the density and velocity fields along the line of
-sight.
-"""
 from os.path import join
 
 import csiborgtools
@@ -23,53 +20,61 @@ from astropy.coordinates import SkyCoord, angular_separation
 from h5py import File
 from mpi4py import MPI
 from scipy.interpolate import RegularGridInterpolator
+from tqdm import trange
 
 from field_los import get_los
 
 
-def interpolate_indranil_void(kind, RA, dec, rmax, dr, dump_folder, catalogue):
+def interpolate_indranil_void(kind, nsims, RA, dec, rmax, dr, dump_folder,
-    fdir = "/mnt/extraspace/rstiskalek/catalogs"
+                              catalogue):
-    if kind == "exp":
+    if kind not in ["exp", "gauss"]:
-        fname = join(fdir, "v_pec_EXP_IndranilVoid.dat")
+        raise ValueError(f"Unknown void kind: `{kind}`.")
-    elif kind == "gauss":
+
-        fname = join(fdir, "v_pec_GAUSS_IndranilVoid.dat")
+    kind = kind.upper()
-    else:
+    fdir = join("/mnt/extraspace/rstiskalek/catalogs", "IndranilVoid",
-        raise ValueError("Invalid void kind.")
+                f"{kind}profile")
 
-    # These are only velocities.
-    data = np.loadtxt(fname)
     fname_out = join(dump_folder, f"los_{catalogue}_IndranilVoid_{kind}.hdf5")
-
-    r_grid = np.arange(0, 251)
-    phi_grid = np.arange(0, 181)
-    r_eval = np.arange(0, rmax, dr).astype(float) / 0.674
-
-    model_axis = SkyCoord(l=117, b=4, frame='galactic', unit='deg').icrs
-    coords = SkyCoord(ra=RA, dec=dec, unit='deg').icrs
-
-    # Get angular separation in degrees
-    phi = angular_separation(coords.ra.rad, coords.dec.rad,
-                             model_axis.ra.rad, model_axis.dec.rad)
-    phi *= 180 / np.pi
-
-    # Get the interpolator
-    f = RegularGridInterpolator((r_grid, phi_grid), data.T)
-    # Get the dummy x-values to evaluate for each LOS
-    x_dummy = np.ones((len(r_eval), 2))
-    x_dummy[:, 0] = r_eval
-
-    result = np.full((len(RA), len(r_eval)), np.nan)
-    for i in range(len(RA)):
-        x_dummy[:, 1] = phi[i]
-        result[i] = f(x_dummy)
-
-    # Write the output, homogenous density.
-    density = np.ones_like(result)
     print(f"Writing to `{fname_out}`.")
-    with File(fname_out, 'w') as f_out:
+    for k in trange(len(nsims), desc="LG observers"):
-        f_out.create_dataset("rdist_0", data=r_eval * 0.674)
+        nsim = nsims[k]
-        f_out.create_dataset("density_0", data=density)
+        # These are only velocities.
-        f_out.create_dataset("velocity_0", data=result)
+        fname = join(fdir, f"v_pec_{kind}profile_rLG_{nsim}.dat")
+        data = np.loadtxt(fname)
+
+        # The grid is in Mpc
+        r_grid = np.arange(0, 251)
+        # NOTE: The shape of the files is no longer (181, 251). It is now
+        # (180, 251), asked Sergij about this. He will produce new files.
+        phi_grid = np.arange(0, len(data))
+        # The input is in Mpc/h, so we need to convert to Mpc
+        r_eval = np.arange(0, rmax, dr).astype(float) / 0.674
+
+        model_axis = SkyCoord(l=117, b=4, frame='galactic', unit='deg').icrs
+        coords = SkyCoord(ra=RA, dec=dec, unit='deg').icrs
+
+        # Get angular separation in degrees
+        phi = angular_separation(coords.ra.rad, coords.dec.rad,
+                                 model_axis.ra.rad, model_axis.dec.rad)
+        phi *= 180 / np.pi
+
+        # Get the interpolator
+        f = RegularGridInterpolator((r_grid, phi_grid), data.T)
+        # Get the dummy x-values to evaluate for each LOS
+        x_dummy = np.ones((len(r_eval), 2))
+        x_dummy[:, 0] = r_eval
+
+        result = np.full((len(RA), len(r_eval)), np.nan)
+        for i in range(len(RA)):
+            x_dummy[:, 1] = phi[i]
+            result[i] = f(x_dummy)
+
+        # Write the output, homogenous density.
+        density = np.ones_like(result)
+        with File(fname_out, 'w') as f_out:
+            f_out.create_dataset(f"rdist_{k}", data=r_eval * 0.674)
+            f_out.create_dataset(f"density_{k}", data=density)
+            f_out.create_dataset(f"velocity_{k}", data=result)
 
 
 ###############################################################################
@@ -84,6 +89,7 @@ if __name__ == "__main__":
 
     comm = MPI.COMM_WORLD
     paths = csiborgtools.read.Paths(**csiborgtools.paths_glamdring)
+    nsims = paths.get_ics(f"IndranilVoid_{kind}")
 
     out_folder = "/mnt/extraspace/rstiskalek/csiborg_postprocessing/field_los"
 
@@ -92,4 +98,4 @@ if __name__ == "__main__":
 
         RA, dec = get_los(catalogue, "", comm).T
         interpolate_indranil_void(
-            kind, RA, dec, rmax, dr, out_folder, catalogue)
+            kind, nsims, RA, dec, rmax, dr, out_folder, catalogue)

scripts/flow/flow_validation.py

@@ -319,6 +319,9 @@ if __name__ == "__main__":
     if mag_selection and inference_method != "bayes":
         raise ValueError("Magnitude selection is only supported with `bayes` inference.")   # noqa
 
+    if "IndranilVoid" in ARGS.simname and ARGS.ksim is None:
+        raise ValueError("`IndranilVoid` must be run only per specific realization.")       # noqa
+
     if inference_method != "bayes":
         mag_selection = [None] * len(ARGS.catalogue)
     elif mag_selection is None or mag_selection:

scripts/flow/flow_validation.sh

@@ -37,24 +37,24 @@ else
 fi
 
 
-for simname in "Carrick2015" "csiborg2_main"; do
+for simname in "IndranilVoid_gauss"; do
 # for simname in "csiborg2_main" "csiborg2X" ; do
 # for simname in "Carrick2015" "Lilow2024" "csiborg2_main" "csiborg2X" "CF4"; do
 # for simname in "Carrick2015" "csiborg2X" "csiborg2_main"; do
 # for simname in "Carrick2015"; do
     # for catalogue in "LOSS" "Foundation" "2MTF" "SFI_gals" "CF4_TFR_i" "CF4_TFR_w1"; do
-    for catalogue in "2MTF" "SFI_gals" "CF4_TFR_i"; do
+    for catalogue in "LOSS"; do
     # for catalogue in "2MTF" "SFI" "CF4_TFR_not2MTForSFI_i"; do
     # for catalogue in "2MTF" "SFI_gals" "CF4_TFR_i"; do
     # for catalogue in "CF4_TFR_w1"; do
     # for catalogue in "CF4_GroupAll"; do
         for ksim in "none"; do
-        for ksmooth in 1 2 3 4; do
+        # for ksim in {0..500}; do
             pythoncm="$env $file --catalogue $catalogue --simname $simname --ksim $ksim --ksmooth $ksmooth --ndevice $ndevice --device $device"
 
             if [ "$on_login" == "1" ]; then
                 echo $pythoncm
-                eval $pythoncm
+                # eval $pythoncm
             else
                 if [ "$device" == "gpu" ]; then
                     cm="addqueue -q $queue -s -m $memory --gpus 1 --gputype $gputype $pythoncm"
@@ -63,14 +63,12 @@ for simname in "Carrick2015" "csiborg2_main"; do
                 fi
                 echo "Submitting:"
                 echo $cm
-                eval $cm
+                # eval $cm
             fi
 
             echo
-            sleep 0.001
+            sleep 0.01
 
         done
     done
-done
-
 done