Examples¶
Four examples, each a C++ program and a Python script that do the same thing: the same inputs, parameters and seed, and the same output file, byte for byte. Each folder in examples is a CMake project of its own, built against an installed otswap.
The data¶
The examples read the files of examples/data, which come with a clone of the repository; they are not in the Python package or the source distribution. They are provided for running the examples, may be redistributed unchanged with them, and may not be used for any other purpose: examples/data/README.md gives the terms and the content of each file.
Building and running¶
For an example <name>, with otswap installed under <prefix>:
cmake -S examples/<name> -B build/examples/<name> -DCMAKE_PREFIX_PATH=<prefix>
cmake --build build/examples/<name>
build/examples/<name>/<name> [data_dir]
python examples/<name>/<name>.py [data_dir]
Both read the examples’ data folder unless data_dir is given, and write to
examples/<name>/output/, created if missing. Where cfitsio is not in a
system path, give its prefix with -DCFITSIO_ROOT=... (and -DZLIB_ROOT=...
for a static cfitsio).
Box¶
Reconstructs the displacement field of halos in a cube, in real space, with box geometry: Cartesian coordinates and a mean particle separation computed from the tracers. Eight realizations, each on its own subset of the randoms.
Reads:
Writes output/displacement_box.dat, with otswap.io.write_displacements(); its columns
and keywords are on the file formats page.
# SPDX-License-Identifier: GPL-2.0-or-later
# Copyright (C) 2026 by Simone Sartori, simone.sartori@inaf.it
# =================================================================
# Example code: how to run the OT reconstruction with box geometry
# =================================================================
"""This example explains how to reconstruct the displacement field of a
tracer catalogue with box geometry. The box geometry can be used with any
catalogue in comoving Cartesian coordinates, including non-cubic ones,
provided that the mean particle separation is approximately constant
throughout the catalogue.
The reconstruction matches the tracers to randoms by optimal transport: in
each realization every tracer is paired with one random, and its
displacement points from the tracer (the Eulerian, observed position) to
that random (the Lagrangian position). Independent realizations are then
averaged.
The box is not periodic: nothing flows through its faces, so modes on the
scale of the box itself are not reconstructed, and the displacements are
least reliable near the faces. A box cut from a periodic simulation is
treated as a plain sub-volume.
Usage: python box.py [data_dir]. The catalogues are read from data_dir, by
default the examples' data folder; the output is written to the folder
output/ next to this file, created if missing.
"""
import sys
from pathlib import Path
import numpy as np
import otswap
def main():
# -------------------------------------------------------------
# ------------ Locate the input and output folders ------------
# -------------------------------------------------------------
# The data folder is the examples' one, next to this example's folder,
# unless a first command-line argument replaces it
here = Path(__file__).resolve().parent
data = Path(sys.argv[1]) if len(sys.argv) > 1 else here.parent / "data"
output = here / "output"
output.mkdir(parents=True, exist_ok=True)
# ---------------------------------------------------------------
# ------------ Read the tracer and random catalogues ------------
# ---------------------------------------------------------------
# otswap.io.read reads the listed columns of a table: for ASCII files
# 0-based indices, or the names of the "###" line of a file otswap
# wrote; for FITS files the column names. Table.values has shape (N, 3),
# x, y, z for each object; coordinates are comoving, in Mpc/h
tracers = otswap.io.read(data / "box_halos_real_space.dat", [0, 1, 2]).values
randoms = otswap.io.read(data / "box_randoms.dat", [0, 1, 2]).values
print(f"Tracers: {len(tracers)}")
print(f"Randoms: {len(randoms)}, {len(randoms) / len(tracers):g} per tracer")
# ------------------------------------------------------------------
# ------------ Set the parameters of the reconstruction ------------
# ------------------------------------------------------------------
# n_realizations: independent realizations. Each one consumes N
# randoms, disjoint from those of the others, so the random catalogue
# must hold at least n_realizations x N objects: the catalogue read here
# has exactly 8 per tracer. Omitting the randoms lets otswap draw them
# itself, uniformly in the tracers' bounding box.
#
# convergence: the sweeps stop once the fraction of successful swaps in a
# sweep, swaps per tracer visited, no longer exceeds this threshold.
# Larger values are faster and less accurate; the default is 1e-3.
#
# seed: a fixed seed makes the run reproducible; 0 draws a new seed at
# each run, and Result.seed records the one drawn. The number of threads
# follows OMP_NUM_THREADS, and with a fixed seed the result is the same
# whatever the number of threads.
#
# cell_size: the grid cell, in units of mps, affects the speed only,
# never the result; the default is 4.
#
# verbosity: what the reconstruction prints; "detailed" adds, before the
# line of the call, a line with the mean particle separation and its
# source
options = dict(n_realizations=8, convergence=1e-2, seed=12345, cell_size=4.0,
verbosity="detailed")
# ------------------------------------------------
# ------------ Run the reconstruction ------------
# ------------------------------------------------
# In box geometry the mean particle separation (mps) is a single number.
# Without one, as here, otswap takes (V/N)^(1/3), with V the volume of
# the tracers' bounding box and N their number, and records it in
# Result.mps; the keyword mps gives it instead. It sets the scale of the
# search for swap partners, so it must describe the whole catalogue:
# with a strongly varying density use the lightcone geometry, which
# measures it as a function of redshift
result = otswap.reconstruct_box(tracers, randoms, **options)
# ----------------------------------------------
# ------------ Summarize the result ------------
# ----------------------------------------------
# In a box result every realization is valid for every tracer, and
# mean_displacement is the average over all of them. Its average over
# the tracers is the offset between the centres of mass of the randoms
# used and of the tracers
mean = result.mean_displacement
all_valid = int(np.count_nonzero(result.valid_realizations == result.n_realizations))
print(f"Realizations: {result.n_realizations}, all valid for {all_valid} tracers")
print(f"Mean length of the mean displacement: {np.linalg.norm(mean, axis=1).mean():g} Mpc/h")
print("Average of the mean displacement over the tracers: ({:g}, {:g}, {:g}) Mpc/h".format(*mean.mean(axis=0)))
print("Displacement of the first tracer in the first realization: ({:g}, {:g}, {:g}) Mpc/h"
.format(*result.displacement[0, 0]))
# ------------------------------------------
# ------------ Write the output ------------
# ------------------------------------------
# One row per tracer, in the order of the input: its position, its mean
# Lagrangian position (the position plus the mean displacement), the
# mean displacement and the number of valid realizations, with a header
# giving the units and the parameters of the run. The writer is the C++
# library's, so the file is the C++ example's, byte for byte;
# otswap.io.write_displacement_field writes every realization instead,
# losslessly
otswap.io.write_displacements(output / "displacement_box.dat", result)
print("Written: output/displacement_box.dat")
if __name__ == "__main__":
main()
// SPDX-License-Identifier: GPL-2.0-or-later
// Copyright (C) 2026 by Simone Sartori, simone.sartori@inaf.it
// =================================================================
// Example code: how to run the OT reconstruction with box geometry
// =================================================================
/*
This example explains how to reconstruct the displacement field of a
tracer catalogue with box geometry. The box geometry can be used with
any catalogue in comoving Cartesian coordinates, including non-cubic
ones, provided that the mean particle separation is approximately
constant throughout the catalogue.
The reconstruction matches the tracers to randoms by optimal
transport: in each realization every tracer is paired with one
random, and its displacement points from the tracer (the Eulerian,
observed position) to that random (the Lagrangian position).
Independent realizations are then averaged.
The box is not periodic: nothing flows through its faces, so modes on
the scale of the box itself are not reconstructed, and the
displacements are least reliable near the faces. A box cut from a
periodic simulation is treated as a plain sub-volume.
Usage: box [data_dir]. The catalogues are read from data_dir, by
default the examples' data folder; the output is written to the
folder output/ next to this file, created if missing.
*/
#include <cmath>
#include <filesystem>
#include <iostream>
#include <string>
#include <vector>
#include "otswap/OT.h"
// ==========================================================================
int main (int argc, char** argv)
{
try {
// -------------------------------------------------------------
// ------------ Locate the input and output folders ------------
// -------------------------------------------------------------
// OTSWAP_EXAMPLES_DATA and OTSWAP_EXAMPLE_OUTPUT are set by this
// example's CMakeLists.txt to the data folder and to output/ next to
// this file. A first command-line argument replaces the data folder
const std::string data = (argc > 1) ? argv[1] : OTSWAP_EXAMPLES_DATA;
const std::string output = OTSWAP_EXAMPLE_OUTPUT;
std::filesystem::create_directories(output);
// ---------------------------------------------------------------
// ------------ Read the tracer and random catalogues ------------
// ---------------------------------------------------------------
// io::read reads the listed columns of a table: for ASCII files
// 0-based indices, or the names of the "###" line of a file otswap
// wrote; for FITS files (.fits, .fit, .fits.gz) the column names.
// Table::values holds the rows one after the other, x, y, z for each
// object: the flat layout every otswap function takes. Coordinates
// are comoving, in Mpc/h
const std::vector<double> tracers = otswap::io::read(data + "/box_halos_real_space.dat", {"0", "1", "2"}).values;
const std::vector<double> randoms = otswap::io::read(data + "/box_randoms.dat", {"0", "1", "2"}).values;
const std::size_t nTracers = tracers.size() / 3;
const std::size_t nRandoms = randoms.size() / 3;
std::cout << "Tracers: " << nTracers << std::endl;
std::cout << "Randoms: " << nRandoms << ", " << (double)nRandoms / (double)nTracers << " per tracer" << std::endl;
// ------------------------------------------------------------------
// ------------ Set the parameters of the reconstruction ------------
// ------------------------------------------------------------------
otswap::Config config;
// Independent realizations. Each one consumes N randoms, disjoint
// from those of the others, so the random catalogue must hold at least
// nRealizations x N objects: the catalogue read here has exactly 8 per
// tracer. An empty random array lets otswap draw the randoms itself,
// uniformly in the tracers' bounding box
config.nRealizations = 8;
// The sweeps stop once the fraction of successful swaps in a sweep,
// swaps per tracer visited, no longer exceeds this threshold. Larger
// values are faster and less accurate; the default is 1e-3
config.convergence = 1.e-2;
// A fixed seed makes the run reproducible; 0 draws a new seed at each
// run, and Result::config.seed records the one drawn. The number of
// threads follows OMP_NUM_THREADS, and with a fixed seed the result is
// the same whatever the number of threads
config.seed = 12345;
// The grid cell, in units of mps, affects the speed only, never the
// result; the default is 4
config.cellSize = 4.;
// What the reconstruction reports, and where: Detailed adds, before the
// line of the call, a line with the mean particle separation and its
// source. The report goes to std::clog by default; here to std::cout, in order
// with this program's own lines
config.verbosity = otswap::Verbosity::Detailed;
config.log = &std::cout;
// ------------------------------------------------
// ------------ Run the reconstruction ------------
// ------------------------------------------------
// In box geometry the mean particle separation (mps) is a single
// number. Without one, as here, otswap takes (V/N)^(1/3), with V the
// volume of the tracers' bounding box and N their number, and records
// it in Result::mps; an overload takes it from the caller. It sets the
// scale of the search for swap partners, so it must describe the whole
// catalogue: with a strongly varying density use the lightcone
// geometry, which measures it as a function of redshift
const otswap::Result result = otswap::reconstructBox(tracers, randoms, config);
// ----------------------------------------------
// ------------ Summarize the result ------------
// ----------------------------------------------
// In a box result every realization is valid for every tracer, and
// meanDisplacement is the average over all of them. Its average over
// the tracers is the offset between the centres of mass of the
// randoms used and of the tracers
std::size_t allValid = 0;
for (std::size_t i = 0; i < nTracers; ++i)
allValid += result.validRealizations[i] == result.nRealizations ? 1 : 0;
std::cout << "Realizations: " << result.nRealizations << ", all valid for " << allValid << " tracers" << std::endl;
double length = 0., average[3] = {0., 0., 0.};
for (std::size_t i = 0; i < nTracers; ++i) {
const double* d = &result.meanDisplacement[3*i];
length += std::sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
for (int c = 0; c < 3; ++c) average[c] += d[c];
}
std::cout << "Mean length of the mean displacement: " << length / (double)nTracers << " Mpc/h" << std::endl;
std::cout << "Average of the mean displacement over the tracers: (" << average[0] / (double)nTracers << ", "
<< average[1] / (double)nTracers << ", " << average[2] / (double)nTracers << ") Mpc/h" << std::endl;
std::cout << "Displacement of the first tracer in the first realization: (" << result.displacement[0] << ", "
<< result.displacement[1] << ", " << result.displacement[2] << ") Mpc/h" << std::endl;
// ------------------------------------------
// ------------ Write the output ------------
// ------------------------------------------
// One row per tracer, in the order of the input: its position, its
// mean Lagrangian position (the position plus the mean displacement),
// the mean displacement and the number of valid realizations, with a
// header giving the units and the parameters of the run. The format
// follows the extension, ASCII here; io::writeDisplacementField writes
// every realization instead, losslessly
otswap::io::writeDisplacements(output + "/displacement_box.dat", result);
std::cout << "Written: output/displacement_box.dat" << std::endl;
}
catch (const otswap::Error& e) { std::cerr << e.what() << std::endl; return 1; }
return 0;
}
Lightcone¶
Reconstructs the displacement field of a survey-like catalogue with lightcone geometry: sky coordinates, the mean particle separation measured as a function of redshift, a redshift cut and the survey’s HEALPix mask, with the displacements that cross unobserved pixels rejected.
Reads:
Writes output/displacement_lightcone.dat, with otswap.io.write_displacements(); its columns
and keywords are on the file formats page.
# SPDX-License-Identifier: GPL-2.0-or-later
# Copyright (C) 2026 by Simone Sartori, simone.sartori@inaf.it
# =======================================================================
# Example code: how to run the OT reconstruction with lightcone geometry
# =======================================================================
"""This example explains how to reconstruct the displacement field of a
survey catalogue with lightcone geometry. The lightcone geometry takes
observed coordinates (right ascension, declination, redshift), a mean
particle separation that changes with redshift, which otswap measures from
the tracers, and the survey's footprint, here a HEALPix mask.
The reconstruction matches the tracers to randoms by optimal transport: in
each realization every tracer is paired with one random, and its
displacement points from the tracer (the Eulerian, observed position) to
that random (the Lagrangian position). Independent realizations are then
averaged.
The randoms carry the survey's geometry and selection function, which
otswap does not model: use the survey's own random catalogue.
Usage: python lightcone.py [data_dir]. The catalogues are read from
data_dir, by default the examples' data folder; the output is written to the
folder output/ next to this file, created if missing.
"""
import sys
from pathlib import Path
import numpy as np
import otswap
def main():
# -------------------------------------------------------------
# ------------ Locate the input and output folders ------------
# -------------------------------------------------------------
# The data folder is the examples' one, next to this example's folder,
# unless a first command-line argument replaces it
here = Path(__file__).resolve().parent
data = Path(sys.argv[1]) if len(sys.argv) > 1 else here.parent / "data"
output = here / "output"
output.mkdir(parents=True, exist_ok=True)
# ---------------------------------------------------------------
# ------------ Read the tracer and random catalogues ------------
# ---------------------------------------------------------------
# otswap.io.read reads the listed columns of a table: for ASCII files
# 0-based indices, or the names of the "###" line of a file otswap
# wrote; for FITS files the column names. Table.values has shape (N, 3):
# right ascension and declination in degrees, and redshift. They are
# passed as they are, with angle_unit="deg"
tracers_sky = otswap.io.read(data / "lightcone_tracers.dat", [0, 1, 2]).values
randoms_sky = otswap.io.read(data / "lightcone_randoms.dat", [0, 1, 2]).values
print(f"Tracers read: {len(tracers_sky)}")
print(f"Randoms read: {len(randoms_sky)}, {len(randoms_sky) / len(tracers_sky):g} per tracer")
# ---------------------------------------
# ------------ Read the mask ------------
# ---------------------------------------
# A HEALPix map in RING or NESTED ordering, read from a FITS file: a
# pixel is observed when its value is greater than 0. A map already in
# memory can be passed to otswap.Mask.from_array instead. When a mask is
# available, give it to the reconstruction rather than a sky area:
# otswap then takes the area from it, and leaves out the objects on
# unobserved pixels
mask = otswap.Mask(str(data / "lightcone_mask.fits"))
print(f"Mask: NSIDE {mask.nside}, {mask.sky_area_deg2:g} deg^2")
# ----------------------------------------------------
# ------------ Set the cosmological model ------------
# ----------------------------------------------------
# The comoving distances, from Omega_m, h, w0 and wa of a flat
# cosmology, tabulated by default from z = 0 to 10. The table must cover
# the redshifts of tracers and randoms
distances = otswap.DistanceTable.flat(0.286, 0.7, w0=-1.0, wa=0.0)
# ------------------------------------------------------------------
# ------------ Set the parameters of the reconstruction ------------
# ------------------------------------------------------------------
# n_realizations: independent realizations. Each one consumes as many
# randoms as there are tracers, disjoint from those of the others, so
# the randoms kept must number at least n_realizations times the
# tracers kept.
#
# convergence: the sweeps stop once the fraction of successful swaps in a
# sweep, swaps per tracer visited, no longer exceeds this threshold.
# Larger values are faster and less accurate; the default is 1e-3.
#
# seed: a fixed seed makes the run reproducible; 0 draws a new seed at
# each run, and result.seed records the one drawn. The number of threads
# follows OMP_NUM_THREADS, and with a fixed seed the result is the same
# whatever the number of threads.
#
# cell_size: the grid cell, in units of the mean particle separation,
# affects the speed only, never the result; the default is 4.
#
# reject_crossings, max_unobserved_pixels_crossed: after the
# reconstruction, a displacement whose path on the sky (the great-circle
# arc from the tracer to its random) crosses more than this many
# unobserved pixels of the mask is marked invalid. The count is of
# pixels, so the same threshold means a different angle at a different
# NSIDE. These are the defaults.
#
# verbosity: what the reconstruction prints. With "normal", the
# default, how many tracers and randoms the redshift cut and the mask
# left out, how many displacements the crossing filter rejected, and
# then the time of the call and the tracers left without a valid
# displacement. "detailed" adds the mps(z) profile, "silent" prints
# nothing; the counts stay in result.selection either way.
#
# n_bins: the mean particle separation is measured in n_bins redshift
# bins of equal width. Each bin needs at least 278 tracers: if one is
# refused, the message says how many bins the catalogue supports.
#
# redshift_cut: tracers and randoms outside this redshift range are left
# out of the reconstruction. The range here spans the whole catalogue:
# narrow it to reconstruct a slice. The default, None, cuts nothing
options = dict(n_realizations=8, convergence=1e-2, seed=12345, cell_size=4.0,
reject_crossings=True, max_unobserved_pixels_crossed=0, verbosity="normal",
n_bins=30, redshift_cut=(0.885, 1.10))
# ------------------------------------------------
# ------------ Run the reconstruction ------------
# ------------------------------------------------
# Tracers left out by the cut or the mask keep their row in the result,
# flagged in outside_redshift_cut or outside_mask, with NaN
# displacements. The result also holds the tracers' Cartesian positions
# (tracers), their mean Lagrangian positions (lagrangian, and
# lagrangian_sky on the sky, in the angle unit of the call) and the
# measured mps(z) (mps_profile)
result = otswap.reconstruct_lightcone(tracers_sky, randoms_sky, mask=mask, distances=distances,
angle_unit="deg", **options)
# ----------------------------------------------
# ------------ Summarize the result ------------
# ----------------------------------------------
# valid holds, for each realization and tracer, whether the displacement
# survived the crossing filter; valid_realizations counts them per
# tracer, and mean_displacement averages the valid ones, NaN where there
# are none
n = len(tracers_sky)
for rec in range(result.n_realizations):
print(f"Realization {rec}: {np.count_nonzero(result.valid[rec])} valid displacements of {n}")
counts = np.bincount(result.valid_realizations, minlength=result.n_realizations + 1)
print("Tracers by number of valid realizations: " + ", ".join(f"{k}: {c}" for k, c in enumerate(counts)))
left_out = np.count_nonzero(result.outside_redshift_cut | result.outside_mask)
print(f"Tracers left out (redshift cut or mask): {left_out}")
mean = result.mean_displacement
no_mean = np.isnan(mean[:, 0])
print(f"Tracers with a NaN mean displacement: {np.count_nonzero(no_mean)}")
print(f"Mean length of the mean displacement: {np.linalg.norm(mean[~no_mean], axis=1).mean():g} Mpc/h")
# ------------------------------------------
# ------------ Write the output ------------
# ------------------------------------------
# One row per tracer, in the order of the input: its sky coordinates
# and those of its mean Lagrangian position, in degrees, its Cartesian
# position, its mean Lagrangian position and its mean displacement, in
# Mpc/h, its number of valid realizations and the two selection flags,
# with a header giving the units and the parameters of the run. A tracer
# left out, or without a valid realization, has NaN Lagrangian
# coordinates and displacement. The writer is the C++ library's, so the
# file is the C++ example's, byte for byte.
# otswap.io.write_displacement_field writes every realization instead,
# losslessly, and otswap.io.write_mps_profile the mps(z)
otswap.io.write_displacements(output / "displacement_lightcone.dat", result)
print("Written: output/displacement_lightcone.dat")
if __name__ == "__main__":
main()
// SPDX-License-Identifier: GPL-2.0-or-later
// Copyright (C) 2026 by Simone Sartori, simone.sartori@inaf.it
// =======================================================================
// Example code: how to run the OT reconstruction with lightcone geometry
// =======================================================================
/*
This example explains how to reconstruct the displacement field of a
survey catalogue with lightcone geometry. The lightcone geometry takes
observed coordinates (right ascension, declination, redshift), a mean
particle separation that changes with redshift, which otswap measures
from the tracers, and the survey's footprint, here a HEALPix mask.
The reconstruction matches the tracers to randoms by optimal
transport: in each realization every tracer is paired with one
random, and its displacement points from the tracer (the Eulerian,
observed position) to that random (the Lagrangian position).
Independent realizations are then averaged.
The randoms carry the survey's geometry and selection function, which
otswap does not model: use the survey's own random catalogue.
Usage: lightcone [data_dir]. The catalogues are read from data_dir, by
default the examples' data folder; the output is written to the
folder output/ next to this file, created if missing.
*/
#include <cmath>
#include <filesystem>
#include <iostream>
#include <string>
#include <vector>
#include "otswap/OT.h"
// ==========================================================================
int main (int argc, char** argv)
{
try {
// -------------------------------------------------------------
// ------------ Locate the input and output folders ------------
// -------------------------------------------------------------
// OTSWAP_EXAMPLES_DATA and OTSWAP_EXAMPLE_OUTPUT are set by this
// example's CMakeLists.txt to the data folder and to output/ next to
// this file. A first command-line argument replaces the data folder
const std::string data = (argc > 1) ? argv[1] : OTSWAP_EXAMPLES_DATA;
const std::string output = OTSWAP_EXAMPLE_OUTPUT;
std::filesystem::create_directories(output);
// ---------------------------------------------------------------
// ------------ Read the tracer and random catalogues ------------
// ---------------------------------------------------------------
// io::read reads the listed columns of a table: for ASCII files
// 0-based indices, or the names of the "###" line of a file otswap
// wrote; for FITS files (.fits, .fit, .fits.gz) the column names.
// Table::values holds the rows one after the other: right ascension,
// declination and redshift for each object, the flat layout every
// otswap function takes. The catalogues give the angles in degrees,
// and otswap takes radians: skyToRadians converts them with the factor
// numpy.deg2rad and otswap's Python interface use
const std::vector<double> tracersSky =
otswap::skyToRadians(otswap::io::read(data + "/lightcone_tracers.dat", {"0", "1", "2"}).values);
const std::vector<double> randomsSky =
otswap::skyToRadians(otswap::io::read(data + "/lightcone_randoms.dat", {"0", "1", "2"}).values);
const std::size_t nTracers = tracersSky.size() / 3;
const std::size_t nRandoms = randomsSky.size() / 3;
std::cout << "Tracers read: " << nTracers << std::endl;
std::cout << "Randoms read: " << nRandoms << ", " << (double)nRandoms / (double)nTracers << " per tracer" << std::endl;
// ---------------------------------------
// ------------ Read the mask ------------
// ---------------------------------------
// A HEALPix map in RING or NESTED ordering, read from a FITS file: a
// pixel is observed when its value is greater than 0. A map already in
// memory can be passed to the constructor taking an array instead.
// When a mask is available, give it to the reconstruction rather than
// a sky area: otswap then takes the area from it, and leaves out the
// objects on unobserved pixels
const otswap::Mask mask(data + "/lightcone_mask.fits");
std::cout << "Mask: NSIDE " << mask.nside() << ", " << mask.skyAreaDeg2() << " deg^2" << std::endl;
// ----------------------------------------------------
// ------------ Set the cosmological model ------------
// ----------------------------------------------------
// The comoving distances, from Omega_m, h, w0 and wa of a flat
// cosmology, tabulated by default from z = 0 to 10. The table must
// cover the redshifts of tracers and randoms
const otswap::DistanceTable distances(0.286, 0.7, -1., 0.);
// ------------------------------------------------------------------
// ------------ Set the parameters of the reconstruction ------------
// ------------------------------------------------------------------
otswap::Config config;
// Independent realizations. Each one consumes as many randoms as there
// are tracers, disjoint from those of the others, so the randoms kept
// must number at least nRealizations times the tracers kept
config.nRealizations = 8;
// The sweeps stop once the fraction of successful swaps in a sweep,
// swaps per tracer visited, no longer exceeds this threshold. Larger
// values are faster and less accurate; the default is 1e-3
config.convergence = 1.e-2;
// A fixed seed makes the run reproducible; 0 draws a new seed at each
// run, and Result::config.seed records the one drawn. The number of
// threads follows OMP_NUM_THREADS, and with a fixed seed the result is
// the same whatever the number of threads
config.seed = 12345;
// The grid cell, in units of the mean particle separation, affects the
// speed only, never the result; the default is 4
config.cellSize = 4.;
// After the reconstruction, a displacement whose path on the sky (the
// great-circle arc from the tracer to its random) crosses more than
// maxUnobservedPixelsCrossed unobserved pixels of the mask is marked
// invalid. The count is of pixels, so the same threshold means a
// different angle at a different NSIDE. These are the defaults
config.rejectCrossings = true;
config.maxUnobservedPixelsCrossed = 0;
// What the reconstruction reports: with Normal, the default, how many
// tracers and randoms the redshift cut and the mask left out, how many
// displacements the crossing filter rejected, and then the time of the
// call and the tracers left without a valid displacement. Detailed adds
// the mps(z) profile, Silent prints nothing; the counts stay in
// Result::selection either way. The report goes to std::clog by
// default; here to std::cout, in order with this program's own lines
config.verbosity = otswap::Verbosity::Normal;
config.log = &std::cout;
// The mean particle separation is measured in nBins redshift bins of
// equal width. Each bin needs at least 278 tracers: if one is refused,
// the message says how many bins the catalogue supports
const unsigned nBins = 30;
// Tracers and randoms outside this redshift range are left out of the
// reconstruction. The range here spans the whole catalogue: narrow it
// to reconstruct a slice. The default, {}, cuts nothing
const otswap::RedshiftCut cut {0.885, 1.10};
// ------------------------------------------------
// ------------ Run the reconstruction ------------
// ------------------------------------------------
// Tracers left out by the cut or the mask keep their row in the
// result, flagged in outsideRedshiftCut or outsideMask, with NaN
// displacements. The result also holds the tracers' Cartesian
// positions (tracers), their mean Lagrangian positions (lagrangian, and
// lagrangianSky on the sky) and the measured mps(z) (mpsProfile)
const otswap::Result result = otswap::reconstructLightcone(tracersSky, randomsSky, mask, nBins,
distances, config, cut);
// ----------------------------------------------
// ------------ Summarize the result ------------
// ----------------------------------------------
// valid holds, for each realization and tracer, whether the
// displacement survived the crossing filter; validRealizations counts
// them per tracer, and meanDisplacement averages the valid ones, NaN
// where there are none
for (unsigned rec = 0; rec < result.nRealizations; ++rec) {
std::size_t n = 0;
for (std::size_t i = 0; i < nTracers; ++i) n += result.valid[rec * nTracers + i];
std::cout << "Realization " << rec << ": " << n << " valid displacements of " << nTracers << std::endl;
}
std::vector<std::size_t> byCount(result.nRealizations + 1, 0);
std::size_t leftOut = 0, noMean = 0;
double length = 0.;
for (std::size_t i = 0; i < nTracers; ++i) {
++byCount[result.validRealizations[i]];
leftOut += (result.outsideRedshiftCut[i] || result.outsideMask[i]) ? 1 : 0;
const double* d = &result.meanDisplacement[3*i];
if (std::isnan(d[0])) ++noMean;
else length += std::sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
}
std::cout << "Tracers by number of valid realizations:";
for (std::size_t n = 0; n < byCount.size(); ++n)
std::cout << (n ? ", " : " ") << n << ": " << byCount[n];
std::cout << std::endl;
std::cout << "Tracers left out (redshift cut or mask): " << leftOut << std::endl;
std::cout << "Tracers with a NaN mean displacement: " << noMean << std::endl;
std::cout << "Mean length of the mean displacement: " << length / (double)(nTracers - noMean) << " Mpc/h" << std::endl;
// ------------------------------------------
// ------------ Write the output ------------
// ------------------------------------------
// One row per tracer, in the order of the input: its sky coordinates
// and those of its mean Lagrangian position, in degrees, its Cartesian
// position, its mean Lagrangian position and its mean displacement, in
// Mpc/h, its number of valid realizations and the two selection flags,
// with a header giving the units and the parameters of the run. A
// tracer left out, or without a valid realization, has NaN Lagrangian
// coordinates and displacement. io::writeDisplacementField writes every
// realization instead, losslessly, and io::writeMpsProfile the mps(z)
otswap::io::writeDisplacements(output + "/displacement_lightcone.dat", result);
std::cout << "Written: output/displacement_lightcone.dat" << std::endl;
}
catch (const otswap::Error& e) { std::cerr << e.what() << std::endl; return 1; }
return 0;
}
Redshift-space correction in a box¶
Runs the box reconstruction on the same halos in redshift space, distorted along one axis, then moves each halo along that axis to its real-space position with the factor f/(b + 3f/5).
Reads:
Writes output/reconstructed_catalogue_box.dat, with otswap.io.write_real_space_catalog(); its columns
and keywords are on the file formats page.
# SPDX-License-Identifier: GPL-2.0-or-later
# Copyright (C) 2026 by Simone Sartori, simone.sartori@inaf.it
# =====================================================================
# Example code: how to move a box catalogue from redshift to real space
# =====================================================================
"""This example explains how to correct a catalogue in box geometry for
redshift-space distortions with the OT reconstruction. The tracers are given
in redshift space, distorted along one Cartesian axis, the line of sight.
The example runs the reconstruction of the box example on them, then moves
each tracer along that axis to the position it would have in real space.
The correction of tracer i is a shift along the axis,
s_i = f(z) / (b + 3 f(z) / 5) * <Psi_axis>_i ,
with f the linear growth rate at the box's redshift z, b the tracers' linear
bias, and <Psi_axis>_i the component along the axis of the reconstructed
mean displacement, averaged with a gaussian of width sigma over the tracers
around i.
The box is not periodic, here or in the reconstruction: nothing flows
through its faces, so modes on the scale of the box itself are not
reconstructed, and the corrected positions are not wrapped. A box cut from a
periodic simulation is treated as a plain sub-volume.
Usage: python full_recons_box.py [data_dir]. The catalogues are read from
data_dir, by default the examples' data folder; the output is written to the
folder output/ next to this file, created if missing.
"""
import sys
from pathlib import Path
import numpy as np
import otswap
def main():
# -------------------------------------------------------------
# ------------ Locate the input and output folders ------------
# -------------------------------------------------------------
# The data folder is the examples' one, next to this example's folder,
# unless a first command-line argument replaces it
here = Path(__file__).resolve().parent
data = Path(sys.argv[1]) if len(sys.argv) > 1 else here.parent / "data"
output = here / "output"
output.mkdir(parents=True, exist_ok=True)
# ---------------------------------------------------------------
# ------------ Read the tracer and random catalogues ------------
# ---------------------------------------------------------------
# Comoving Cartesian coordinates in Mpc/h, shape (N, 3). The tracers are
# in redshift space, with the line of sight along y; the randoms are
# uniform in the volume
tracers = otswap.io.read(data / "box_halos_redshift_space.dat", [0, 1, 2]).values
randoms = otswap.io.read(data / "box_randoms.dat", [0, 1, 2]).values
print(f"Tracers: {len(tracers)}")
print(f"Randoms: {len(randoms)}, {len(randoms) / len(tracers):g} per tracer")
# ------------------------------------------------
# ------------ Run the reconstruction ------------
# ------------------------------------------------
# As in the box example: 8 realizations, a convergence threshold of
# 1e-2, a fixed seed, and the mean particle separation of the tracers'
# bounding box, computed by otswap and printed with "detailed"
result = otswap.reconstruct_box(tracers, randoms, n_realizations=8, convergence=1e-2,
seed=12345, verbosity="detailed")
# ----------------------------------------------------
# ------------ Set the cosmological model ------------
# ----------------------------------------------------
# Only the linear growth rate f at the box's redshift is needed. The
# table computes it from Omega_m, h, w0 and wa of a flat cosmology
distances = otswap.DistanceTable.flat(0.3186, 0.67, w0=-1.0, wa=0.0)
# ----------------------------------------------------------------
# ------------ Correct the redshift-space distortions ------------
# ----------------------------------------------------------------
# The box's redshift, the tracers' linear bias, and the line of sight:
# axis 0, 1 or 2 for x, y or z
redshift, bias, axis = 1.0, 1.91255, 1
# sigma, the width of the gaussian average in Mpc/h: 10 Mpc/h is a
# reasonable starting value, but the best value depends on the sample
# and should be checked in each analysis. With
# weight_by_realizations=True, each neighbour would be weighted by its
# number of valid realizations.
#
# The correction reads the tracers from the result. A tracer with no
# valid neighbour within 3 sigma, which in a box happens only for an
# isolated one, is left uncorrected, with NaN coordinates;
# catalogue.status says which
catalogue = otswap.real_space_box(result, axis=axis, redshift=redshift, distances=distances,
bias=bias, sigma=10.0, weight_by_realizations=False)
# ----------------------------------------------
# ------------ Summarize the result ------------
# ----------------------------------------------
# The single factor f/(b + 3f/5) of the box, and the shift applied along
# the line of sight to each corrected tracer, both kept in the catalogue
print(f"RSD factor f/(b + 3f/5): {catalogue.factor[0]:g}")
shift = catalogue.shift[~np.isnan(catalogue.shift)]
print(f"Shift along the line of sight over the {len(shift)} corrected tracers: mean {shift.mean():g}, "
f"rms {np.sqrt(np.mean(shift ** 2)):g} Mpc/h")
# ------------------------------------------
# ------------ Write the output ------------
# ------------------------------------------
# One row per tracer, in the order of the input: the corrected position,
# NaN for an uncorrected tracer, then the number of valid realizations,
# the two diagnostics of the average and the status of the tracer. The
# writer is the C++ library's, so the file is the C++ example's, byte
# for byte
otswap.io.write_real_space_catalog(output / "reconstructed_catalogue_box.dat", catalogue)
print("Written: output/reconstructed_catalogue_box.dat")
if __name__ == "__main__":
main()
// SPDX-License-Identifier: GPL-2.0-or-later
// Copyright (C) 2026 by Simone Sartori, simone.sartori@inaf.it
// =====================================================================
// Example code: how to move a box catalogue from redshift to real space
// =====================================================================
/*
This example explains how to correct a catalogue in box geometry for
redshift-space distortions with the OT reconstruction. The tracers
are given in redshift space, distorted along one Cartesian axis, the
line of sight. The example runs the reconstruction of the box example
on them, then moves each tracer along that axis to the position it
would have in real space.
The correction of tracer i is a shift along the axis,
s_i = f(z) / (b + 3 f(z) / 5) * <Psi_axis>_i ,
with f the linear growth rate at the box's redshift z, b the tracers'
linear bias, and <Psi_axis>_i the component along the axis of the
reconstructed mean displacement, averaged with a gaussian of width
sigma over the tracers around i.
The box is not periodic, here or in the reconstruction: nothing flows
through its faces, so modes on the scale of the box itself are not
reconstructed, and the corrected positions are not wrapped. A box cut
from a periodic simulation is treated as a plain sub-volume.
Usage: full_recons_box [data_dir]. The catalogues are read from
data_dir, by default the examples' data folder; the output is written
to the folder output/ next to this file, created if missing.
*/
#include <cmath>
#include <filesystem>
#include <iostream>
#include <string>
#include <vector>
#include "otswap/OT.h"
// ==========================================================================
int main (int argc, char** argv)
{
try {
// -------------------------------------------------------------
// ------------ Locate the input and output folders ------------
// -------------------------------------------------------------
// OTSWAP_EXAMPLES_DATA and OTSWAP_EXAMPLE_OUTPUT are set by this
// example's CMakeLists.txt to the data folder and to output/ next to
// this file. A first command-line argument replaces the data folder
const std::string data = (argc > 1) ? argv[1] : OTSWAP_EXAMPLES_DATA;
const std::string output = OTSWAP_EXAMPLE_OUTPUT;
std::filesystem::create_directories(output);
// ---------------------------------------------------------------
// ------------ Read the tracer and random catalogues ------------
// ---------------------------------------------------------------
// Comoving Cartesian coordinates in Mpc/h, x, y, z for each object,
// one row after the other. The tracers are in redshift space, with
// the line of sight along y; the randoms are uniform in the volume
const std::vector<double> tracers = otswap::io::read(data + "/box_halos_redshift_space.dat", {"0", "1", "2"}).values;
const std::vector<double> randoms = otswap::io::read(data + "/box_randoms.dat", {"0", "1", "2"}).values;
const std::size_t nTracers = tracers.size() / 3;
const std::size_t nRandoms = randoms.size() / 3;
std::cout << "Tracers: " << nTracers << std::endl;
std::cout << "Randoms: " << nRandoms << ", " << (double)nRandoms / (double)nTracers << " per tracer" << std::endl;
// ------------------------------------------------
// ------------ Run the reconstruction ------------
// ------------------------------------------------
// As in the box example: 8 realizations, a convergence threshold of
// 1e-2, a fixed seed, and the mean particle separation of the tracers'
// bounding box, computed by otswap and reported with Detailed. The
// report goes to std::cout, in order with this program's lines
otswap::Config config;
config.nRealizations = 8;
config.convergence = 1.e-2;
config.seed = 12345;
config.verbosity = otswap::Verbosity::Detailed;
config.log = &std::cout;
const otswap::Result result = otswap::reconstructBox(tracers, randoms, config);
// ----------------------------------------------------
// ------------ Set the cosmological model ------------
// ----------------------------------------------------
// Only the linear growth rate f at the box's redshift is needed. The
// table computes it from Omega_m, h, w0 and wa of a flat cosmology
const otswap::DistanceTable distances(0.3186, 0.67, -1., 0.);
// ----------------------------------------------------------------
// ------------ Correct the redshift-space distortions ------------
// ----------------------------------------------------------------
// The box's redshift, the tracers' linear bias, and the line of sight:
// axis 0, 1 or 2 for x, y or z
const double redshift = 1.;
const double bias = 1.91255;
const unsigned axis = 1;
// sigma, the width of the gaussian average in Mpc/h: 10 Mpc/h is a
// reasonable starting value, but the best value depends on the sample
// and should be checked in each analysis. With
// weightByRealizations, each neighbour would be weighted by its number
// of valid realizations.
//
// The correction reads the tracers from the result. A tracer with no
// valid neighbour within 3 sigma, which in a box happens only for an
// isolated one, is left uncorrected, with NaN coordinates;
// RealSpaceCatalog::status says which
const double sigma = 10.;
otswap::CorrectionConfig correction;
correction.weightByRealizations = false;
correction.log = &std::cout;
const otswap::RealSpaceCatalog catalogue = otswap::realSpaceBox(result, axis, redshift, distances, bias,
sigma, correction);
// ----------------------------------------------
// ------------ Summarize the result ------------
// ----------------------------------------------
// The single factor f/(b + 3f/5) of the box, and the shift applied
// along the line of sight to each corrected tracer, both kept in the
// catalogue
std::cout << "RSD factor f/(b + 3f/5): " << catalogue.factor[0] << std::endl;
std::size_t corrected = 0;
double sum = 0., squares = 0.;
for (std::size_t i = 0; i < nTracers; ++i) {
if (std::isnan(catalogue.shift[i])) continue;
++corrected;
sum += catalogue.shift[i];
squares += catalogue.shift[i] * catalogue.shift[i];
}
std::cout << "Shift along the line of sight over the " << corrected << " corrected tracers: mean "
<< sum / (double)corrected << ", rms " << std::sqrt(squares / (double)corrected) << " Mpc/h" << std::endl;
// ------------------------------------------
// ------------ Write the output ------------
// ------------------------------------------
// One row per tracer, in the order of the input: the corrected
// position, NaN for an uncorrected tracer, then the number of valid
// realizations, the two diagnostics of the average and the status of
// the tracer
otswap::io::writeRealSpaceCatalog(output + "/reconstructed_catalogue_box.dat", catalogue);
std::cout << "Written: output/reconstructed_catalogue_box.dat" << std::endl;
}
catch (const otswap::Error& e) { std::cerr << e.what() << std::endl; return 1; }
return 0;
}
Redshift-space correction in a lightcone¶
Runs the lightcone reconstruction, without the redshift cut, then moves each tracer along its line of sight to its real-space position, with the growth rate from the distance table and the bias from a b(z) table.
Reads:
Writes output/reconstructed_catalogue_lightcone.dat, with otswap.io.write_real_space_catalog(); its columns
and keywords are on the file formats page.
# SPDX-License-Identifier: GPL-2.0-or-later
# Copyright (C) 2026 by Simone Sartori, simone.sartori@inaf.it
# ===========================================================================
# Example code: how to move a lightcone catalogue from redshift to real space
# ===========================================================================
"""This example explains how to correct a survey catalogue for
redshift-space distortions with the OT reconstruction, in lightcone
geometry. It runs the reconstruction of the lightcone example, then moves
each tracer along its line of sight to the position it would have in real
space.
The correction of tracer i is a shift along its line of sight,
s_i = f(z_i) / (b(z_i) + 3 f(z_i) / 5) * <Psi . r_hat>_i ,
with f the linear growth rate, b the tracers' linear bias, and
<Psi . r_hat>_i the line-of-sight component of the reconstructed mean
displacement, averaged with a gaussian of width sigma over the tracers
around i. The tracer keeps its right ascension and declination, and moves to
the redshift of comoving distance d(z_i) + s_i.
Usage: python full_recons_lightcone.py [data_dir]. The catalogues are read
from data_dir, by default the examples' data folder; the output is written to
the folder output/ next to this file, created if missing.
"""
import sys
from pathlib import Path
import numpy as np
import otswap
def main():
# -------------------------------------------------------------
# ------------ Locate the input and output folders ------------
# -------------------------------------------------------------
# The data folder is the examples' one, next to this example's folder,
# unless a first command-line argument replaces it
here = Path(__file__).resolve().parent
data = Path(sys.argv[1]) if len(sys.argv) > 1 else here.parent / "data"
output = here / "output"
output.mkdir(parents=True, exist_ok=True)
# ---------------------------------------------------------------
# ------------ Read the tracer and random catalogues ------------
# ---------------------------------------------------------------
# Right ascension and declination in degrees, and redshift; passed as
# they are, with angle_unit="deg"
tracers_sky = otswap.io.read(data / "lightcone_tracers.dat", [0, 1, 2]).values
randoms_sky = otswap.io.read(data / "lightcone_randoms.dat", [0, 1, 2]).values
print(f"Tracers read: {len(tracers_sky)}")
print(f"Randoms read: {len(randoms_sky)}, {len(randoms_sky) / len(tracers_sky):g} per tracer")
# ---------------------------------------
# ------------ Read the mask ------------
# ---------------------------------------
mask = otswap.Mask(str(data / "lightcone_mask.fits"))
print(f"Mask: NSIDE {mask.nside}, {mask.sky_area_deg2:g} deg^2")
# ----------------------------------------------------
# ------------ Set the cosmological model ------------
# ----------------------------------------------------
# The table gives the comoving distances and the linear growth rate f(z)
# of the correction, computed from Omega_m, h, w0 and wa. It must cover
# the corrected distances too, which can fall slightly outside the
# catalogue's redshift range: the default range, z = 0 to 10, leaves
# room for that
distances = otswap.DistanceTable.flat(0.286, 0.7, w0=-1.0, wa=0.0)
# ------------------------------------------------
# ------------ Run the reconstruction ------------
# ------------------------------------------------
# As in the lightcone example, without a redshift cut: 8 realizations,
# a convergence threshold of 1e-2, a fixed seed, 30 redshift bins for
# the mean particle separation, and the displacements crossing an
# unobserved pixel of the mask rejected
result = otswap.reconstruct_lightcone(tracers_sky, randoms_sky, mask=mask, n_bins=30,
distances=distances, angle_unit="deg", n_realizations=8,
convergence=1e-2, seed=12345, reject_crossings=True,
max_unobserved_pixels_crossed=0, verbosity="normal")
# ---------------------------------------------
# ------------ Read the bias table ------------
# ---------------------------------------------
# The tracers' linear bias b(z), as nodes in redshift: b is interpolated
# linearly between them, and extrapolated linearly beyond the first and
# last, with an otswap.ExtrapolationWarning saying for how many tracers.
# This table covers the whole catalogue
bias = otswap.io.read_bias_table(data / "lightcone_bias.dat")
# ----------------------------------------------------------------
# ------------ Correct the redshift-space distortions ------------
# ----------------------------------------------------------------
# sigma, the width of the gaussian average in Mpc/h: 10 Mpc/h is a
# reasonable starting value, but the best value depends on the sample
# and should be checked in each analysis. With
# weight_by_realizations=True, each neighbour would be weighted by its
# number of valid realizations. "detailed" prints, besides the line of
# the call, how many tracers moved with their neighbours' average and
# why the others were left uncorrected.
#
# The correction reads the tracers from the result. A tracer without a
# valid realization still receives the average of its neighbours; one
# with no valid neighbour within 3 sigma, or left out of the
# reconstruction by the mask, is left uncorrected, with NaN coordinates;
# catalogue.status says which
catalogue = otswap.real_space_lightcone(result, distances=distances, bias=bias, sigma=10.0,
weight_by_realizations=False, verbosity="detailed")
# ----------------------------------------------
# ------------ Summarize the result ------------
# ----------------------------------------------
# The factor f/(b + 3f/5) at the redshifts of the tracers that took part
# in the reconstruction, and the shift applied along the line of sight
# to each corrected tracer, both kept in the catalogue
factor = catalogue.factor[~np.isnan(catalogue.factor)]
print(f"RSD factor f/(b + 3f/5): from {factor.min():g} to {factor.max():g}")
shift = catalogue.shift[~np.isnan(catalogue.shift)]
print(f"Shift along the line of sight over the {len(shift)} corrected tracers: mean {shift.mean():g}, "
f"rms {np.sqrt(np.mean(shift ** 2)):g} Mpc/h")
# ------------------------------------------
# ------------ Write the output ------------
# ------------------------------------------
# One row per tracer, in the order of the input: the right ascension and
# declination, unchanged, and the corrected redshift; the corrected
# Cartesian position; the number of valid realizations, the two
# diagnostics of the average and the status of the tracer; NaN positions
# for an uncorrected tracer. The writer is the C++ library's, so the file
# is the C++ example's, byte for byte
otswap.io.write_real_space_catalog(output / "reconstructed_catalogue_lightcone.dat", catalogue)
print("Written: output/reconstructed_catalogue_lightcone.dat")
if __name__ == "__main__":
main()
// SPDX-License-Identifier: GPL-2.0-or-later
// Copyright (C) 2026 by Simone Sartori, simone.sartori@inaf.it
// ===========================================================================
// Example code: how to move a lightcone catalogue from redshift to real space
// ===========================================================================
/*
This example explains how to correct a survey catalogue for
redshift-space distortions with the OT reconstruction, in lightcone
geometry. It runs the reconstruction of the lightcone example, then
moves each tracer along its line of sight to the position it would
have in real space.
The correction of tracer i is a shift along its line of sight,
s_i = f(z_i) / (b(z_i) + 3 f(z_i) / 5) * <Psi . r_hat>_i ,
with f the linear growth rate, b the tracers' linear bias, and
<Psi . r_hat>_i the line-of-sight component of the reconstructed mean
displacement, averaged with a gaussian of width sigma over the
tracers around i. The tracer keeps its right ascension and
declination, and moves to the redshift of comoving distance
d(z_i) + s_i.
Usage: full_recons_lightcone [data_dir]. The catalogues are read from
data_dir, by default the examples' data folder; the output is written
to the folder output/ next to this file, created if missing.
*/
#include <algorithm>
#include <cmath>
#include <filesystem>
#include <iostream>
#include <string>
#include <vector>
#include "otswap/OT.h"
// ==========================================================================
int main (int argc, char** argv)
{
try {
// -------------------------------------------------------------
// ------------ Locate the input and output folders ------------
// -------------------------------------------------------------
// OTSWAP_EXAMPLES_DATA and OTSWAP_EXAMPLE_OUTPUT are set by this
// example's CMakeLists.txt to the data folder and to output/ next to
// this file. A first command-line argument replaces the data folder
const std::string data = (argc > 1) ? argv[1] : OTSWAP_EXAMPLES_DATA;
const std::string output = OTSWAP_EXAMPLE_OUTPUT;
std::filesystem::create_directories(output);
// ---------------------------------------------------------------
// ------------ Read the tracer and random catalogues ------------
// ---------------------------------------------------------------
// Right ascension, declination and redshift for each object, one row
// after the other. The angles are read in degrees and converted to
// radians by skyToRadians, with the factor numpy.deg2rad and otswap's
// Python interface use
const std::vector<double> tracersSky =
otswap::skyToRadians(otswap::io::read(data + "/lightcone_tracers.dat", {"0", "1", "2"}).values);
const std::vector<double> randomsSky =
otswap::skyToRadians(otswap::io::read(data + "/lightcone_randoms.dat", {"0", "1", "2"}).values);
const std::size_t nTracers = tracersSky.size() / 3;
const std::size_t nRandoms = randomsSky.size() / 3;
std::cout << "Tracers read: " << nTracers << std::endl;
std::cout << "Randoms read: " << nRandoms << ", " << (double)nRandoms / (double)nTracers << " per tracer" << std::endl;
// ---------------------------------------
// ------------ Read the mask ------------
// ---------------------------------------
const otswap::Mask mask(data + "/lightcone_mask.fits");
std::cout << "Mask: NSIDE " << mask.nside() << ", " << mask.skyAreaDeg2() << " deg^2" << std::endl;
// ----------------------------------------------------
// ------------ Set the cosmological model ------------
// ----------------------------------------------------
// The table gives the comoving distances and the linear growth rate
// f(z) of the correction, computed from Omega_m, h, w0 and wa. It must
// cover the corrected distances too, which can fall slightly outside
// the catalogue's redshift range: the default range, z = 0 to 10,
// leaves room for that
const otswap::DistanceTable distances(0.286, 0.7, -1., 0.);
// ------------------------------------------------
// ------------ Run the reconstruction ------------
// ------------------------------------------------
// As in the lightcone example, without a redshift cut: 8
// realizations, a convergence threshold of 1e-2, a fixed seed, 30
// redshift bins for the mean particle separation, and the
// displacements crossing an unobserved pixel of the mask rejected.
// The report goes to std::cout, in order with this program's lines
otswap::Config config;
config.nRealizations = 8;
config.convergence = 1.e-2;
config.seed = 12345;
config.rejectCrossings = true;
config.maxUnobservedPixelsCrossed = 0;
config.verbosity = otswap::Verbosity::Normal;
config.log = &std::cout;
const otswap::Result result = otswap::reconstructLightcone(tracersSky, randomsSky, mask, 30, distances, config);
// ---------------------------------------------
// ------------ Read the bias table ------------
// ---------------------------------------------
// The tracers' linear bias b(z), as nodes in redshift: b is
// interpolated linearly between them, and extrapolated linearly beyond
// the first and last, which the correction reports and counts in
// RealSpaceCatalog::nExtrapolated. This table covers the whole
// catalogue
const otswap::BiasTable bias = otswap::io::readBiasTable(data + "/lightcone_bias.dat");
// ----------------------------------------------------------------
// ------------ Correct the redshift-space distortions ------------
// ----------------------------------------------------------------
// sigma, the width of the gaussian average in Mpc/h: 10 Mpc/h is a
// reasonable starting value, but the best value depends on the sample
// and should be checked in each analysis. With
// weightByRealizations, each neighbour would be weighted by its number
// of valid realizations. Detailed reports, besides the line of the
// call, how many tracers moved with their neighbours' average and why
// the others were left uncorrected.
//
// The correction reads the tracers from the result. A tracer without a
// valid realization still receives the average of its neighbours; one
// with no valid neighbour within 3 sigma, or left out of the
// reconstruction by the mask, is left uncorrected, with NaN
// coordinates; RealSpaceCatalog::status says which
const double sigma = 10.;
otswap::CorrectionConfig correction;
correction.weightByRealizations = false;
correction.verbosity = otswap::Verbosity::Detailed;
correction.log = &std::cout;
const otswap::RealSpaceCatalog catalogue = otswap::realSpaceLightcone(result, distances, bias, sigma, correction);
// ----------------------------------------------
// ------------ Summarize the result ------------
// ----------------------------------------------
// The factor f/(b + 3f/5) at the redshifts of the tracers that took
// part in the reconstruction, and the shift applied along the line of
// sight to each corrected tracer, both kept in the catalogue
double factorMin = INFINITY, factorMax = -INFINITY, sum = 0., squares = 0.;
std::size_t corrected = 0;
for (std::size_t i = 0; i < nTracers; ++i) {
if (!std::isnan(catalogue.factor[i])) {
factorMin = std::min(factorMin, catalogue.factor[i]);
factorMax = std::max(factorMax, catalogue.factor[i]);
}
if (std::isnan(catalogue.shift[i])) continue;
++corrected;
sum += catalogue.shift[i];
squares += catalogue.shift[i] * catalogue.shift[i];
}
std::cout << "RSD factor f/(b + 3f/5): from " << factorMin << " to " << factorMax << std::endl;
std::cout << "Shift along the line of sight over the " << corrected << " corrected tracers: mean "
<< sum / (double)corrected << ", rms " << std::sqrt(squares / (double)corrected) << " Mpc/h" << std::endl;
// ------------------------------------------
// ------------ Write the output ------------
// ------------------------------------------
// One row per tracer, in the order of the input: the right ascension
// and declination, unchanged, and the corrected redshift; the corrected
// Cartesian position; the number of valid realizations, the two
// diagnostics of the average and the status of the tracer; NaN
// positions for an uncorrected tracer
otswap::io::writeRealSpaceCatalog(output + "/reconstructed_catalogue_lightcone.dat", catalogue);
std::cout << "Written: output/reconstructed_catalogue_lightcone.dat" << std::endl;
}
catch (const otswap::Error& e) { std::cerr << e.what() << std::endl; return 1; }
return 0;
}