C++ reference

Generated from the public headers otswap/OT.h, otswap/RSD.h and otswap/io.h by Doxygen, with only what carries a comment included.

otswap/OT.h

Errors

class Error : public std::runtime_error

Base class of every exception otswap throws on invalid input or failure. std::bad_alloc is not converted, and propagates unchanged, from inside a parallel region too. The neighbour search of the redshift-space correction raises meshsearch::Error, a std::runtime_error, for positions that span no volume.

An Error records where it was constructed: the source file, the function and the line. The location is not part of what(), which holds the message alone. An error raised by a shared helper, or rethrown with context after a parallel region or a lookup, carries the location of that helper or of the rethrow.

Public Functions

inline explicit Error(const std::string &what, const char *file = __builtin_FILE(), const char *function = __builtin_FUNCTION(), int line = __builtin_LINE())
Parameters:

  • what – the message, returned by what()

  • file, function, line – the location; left to their defaults, they are those of the expression that constructs the Error

inline const char *file() const noexcept

Source file in which the Error was constructed, as the compiler named it.

inline const char *function() const noexcept

Name of the function in which the Error was constructed.

inline int line() const noexcept

Line at which the Error was constructed.

Configuration

struct Config

Free parameters of the reconstruction.

The multipliers that govern the seeding radius, the neighbourhood size and the initial pairing are absent: they are internal constants of the algorithm, not parameters.

Public Members

unsigned nRealizations = 1

Independent reconstructions to run. Each consumes nObjects randoms.

double convergence = 1.e-3

The sweeps stop once the fraction of successful swaps in a sweep, swaps per tracer visited, no longer exceeds this threshold. Every tracer is examined in every sweep.

unsigned seed = 0

Seed of the generator. 0 draws one from std::random_device.

double cellSize = 4.

Grid cell size, in units of the representative mean particle separation. Affects speed only, never the result.

bool rejectCrossings = true

With a Mask given to reconstructLightcone, apply rejectMaskCrossings to the result, with maxUnobservedPixelsCrossed, before returning it. Read by the mask overloads of reconstructLightcone only.

unsigned maxUnobservedPixelsCrossed = 0

The threshold of that call; see rejectMaskCrossings. Read by the mask overloads of reconstructLightcone only.

Verbosity verbosity = Verbosity::Normal

What the call reports; see Verbosity.

std::ostream *log = &std::clog

Where the report is written; null writes nothing, as Silent.

enum class otswap::Verbosity

What a reconstruction or a redshift-space correction reports.

Each line starts with “otswap: “ and ends in ‘

’; the counts behind every line are in the object returned, so nothing printed is unavailable to a program.

Values:

enumerator Silent

Nothing.

enumerator Normal

The selection report of a lightcone, when a selection applied (SelectionCounts::message), and the extrapolation of b(z), when it happened; then one line per call with its time and what was lost: 

otswap: reconstructLightcone: <R> realizations of <N> tracers in <T> s; <K> without a valid displacement
otswap: realSpaceLightcone: corrected <C> of <N> tracers in <T> s; <U> left uncorrected

enumerator Detailed

As Normal, with the lines that explain it: the mean particle separation of a box and its source, the mps(z) profile of a lightcone, and the breakdown of the corrected (moved by their neighbours) and of the uncorrected (left out of the reconstruction, or with no valid neighbour).

enum class otswap::Geometry

The geometry a Result was reconstructed in.

Values:

enumerator Box

reconstructBox

enumerator Lightcone

reconstructLightcone

struct RedshiftCut

A closed redshift range for reconstructLightcone. The default, (-inf, +inf), cuts nothing.

Distances

class DistanceTable

Sampled redshift to comoving distance relation, interpolated in both directions, and the linear growth rate on the same grid.

The object is immutable once built and holds no external state, so one instance can serve any number of reconstructions concurrently.

Public Functions

DistanceTable(double OmegaM, double h, double w0, double wa, double zMin = 0., double zMax = 10., unsigned nSamples = 50000)

Build from the parameters of a flat cosmology.

Samples nSamples points over [zMin, zMax] and integrates 1/E(z). Non-flat cosmologies are served by the table constructor below.

With the default range and sampling, for flat LCDM and w0-wa cosmologies near it, the relative error of distanceAt, redshiftAt and growthRateAt is below 1e-6 for z >= 0.01, and that of growthRateAt over the whole range; the absolute error of distanceAt, and that of redshiftAt expressed as a distance along the line of sight, is below 2e-5 Mpc/h over [0, 10]. Below z = 0.01 the relative error of the distances grows, because the distance itself vanishes at z = 0 while the absolute error does not. The table holds three arrays of nSamples doubles.

The growth rate f = dlnD/dlna is sampled on the same grid, by integrating the linear growth equation from deep in matter domination and taking f from the solution. It is not approximated by Omega_m(z) raised to a power: that form is fitted to LCDM and these parameters describe a w0-wa cosmology.

Parameters:

  • OmegaM – matter density parameter; must lie in (0, 1]

  • h – Hubble constant in units of 100 km/s/Mpc; must be positive

  • w0 – dark energy equation of state, constant term

  • wa – dark energy equation of state, evolution term

  • zMin – lower end of the sampled range; must be non-negative

  • zMax – upper end; must exceed zMin

  • nSamples – sampling points, uniform in redshift; must be at least two

Throws:

Error – if any parameter is out of range.

DistanceTable(std::vector<double> redshift, std::vector<double> distance, std::vector<double> growthRate = {})

Build from a table the caller has computed.

Parameters:

  • redshift – strictly increasing, at least two entries

  • distance – strictly increasing, same length, comoving Mpc/h

  • growthRate – f = dlnD/dlna on the same grid, or empty. When empty the table carries no growth rate and growthRateAt throws. Need not be monotonic; every entry must be finite and positive.

Throws:

Error – if either of the first two arrays is not strictly increasing, if the lengths differ, or if fewer than two points are given.

double distanceAt(double z) const

Comoving distance at z. Throws outside the sampled range. A value outside it by at most 8 rounding errors of the range’s largest magnitude is taken to be at the nearer end, so that values computed from the table’s own ends, such as the radius of a position converted at the largest redshift, are accepted. The same holds for redshiftAt and growthRateAt.

double redshiftAt(double distance) const

Redshift at a comoving distance. Throws outside the sampled range, with the end tolerance of distanceAt.

double growthRateAt(double z) const

Growth rate f = dlnD/dlna at z. Throws outside the sampled range, with the end tolerance of distanceAt, and when the table carries no growth rate.

bool hasGrowthRate() const

True when the table carries a growth rate.

std::vector<double> otswap::toCartesian(const std::vector<double> &sky, const DistanceTable &distances)

Convert sky coordinates to Cartesian comoving coordinates.

Object i at right ascension ra, declination dec and redshift z maps to x = d cos(dec) cos(ra), y = d cos(dec) sin(ra), z = d sin(dec), with d = distances.distanceAt(z). reconstructLightcone converts its sky arrays with this function, so Cartesian arrays made with it are exactly the ones that overload uses.

Parameters:

  • sky – sky coordinates, 3 * nObjects entries, ordered right ascension, declination, redshift; angles in radians. May be empty.

  • distances – the redshift to distance relation.

Throws:

Error – if the array’s size is not a multiple of three, if an entry is not finite, if a declination lies outside [-pi/2, pi/2], or if a redshift falls outside the distance table; the message names the object.

Returns:

Cartesian coordinates in the distance table’s unit, Mpc/h, 3 * nObjects entries ordered x, y, z.

std::vector<double> otswap::toSky(const std::vector<double> &cartesian, const DistanceTable &distances)

Convert Cartesian comoving coordinates to sky coordinates: the inverse of toCartesian, with the same observer at the origin.

Object i at (x, y, z) maps to the right ascension atan2(y, x) folded into [0, 2 pi), the declination asin(z / d), and the redshift at which the table’s comoving distance is d = |(x, y, z)|. A position at the origin gets right ascension and declination 0. toSky(toCartesian(sky)) returns sky within a few rounding errors, with the right ascension folded into [0, 2 pi).

Parameters:

  • cartesian – 3 * nObjects entries, ordered x, y, z, in Mpc/h. May be empty. A row holding a NaN gives a NaN row.

  • distances – the redshift to distance relation; its range bounds the distances accepted.

Throws:

Error – if the array’s size is not a multiple of three, if an entry is infinite, or if a distance falls outside the table (beyond its end tolerance); the message names the object and, for the distance, suggests a table covering a wider redshift range.

Returns:

sky coordinates, 3 * nObjects entries, ordered right ascension, declination, redshift; angles in radians.

std::vector<double> otswap::skyToRadians(std::vector<double> skyDegrees)

Right ascension and declination of a sky array from degrees to radians, multiplied by pi/180; redshifts unchanged.

The factor is the double numpy.deg2rad multiplies by, so an array converted here and one converted by the Python package are the same.

Parameters:

skyDegrees – 3 * nObjects entries, ordered right ascension, declination, redshift; angles in degrees. May be empty.

Throws:

Error – if the size is not a multiple of three, or a finite declination lies outside [-90, 90] degrees; the message names the object.

std::vector<double> otswap::skyToDegrees(std::vector<double> skyRadians)

Right ascension and declination of a sky array from radians to degrees, multiplied by 180/pi; redshifts unchanged. A right ascension that rounds to 360 is folded to 0, so that one in [0, 2 pi) stays in [0, 360). NaN stays NaN.

Throws:

Error – if the size is not a multiple of three.

Mask

class Mask

Healpix veto mask, read from a FITS file or built from a full-sky map in memory.

A pixel is allowed (observed) when its value is greater than 0, whatever its magnitude: fractional values, the smallest positive ones and +inf are allowed. A pixel holding 0, a negative value, -inf, NaN or Healpix’s UNSEEN (-1.6375e30) is unobserved. Every value is accepted, and each is tested once, in double, when the mask is built; none is kept or used as a weight, so skyAreaDeg2 is a count of allowed pixels, not a sum of their values. The rule is the same for both sources.

NSIDE may be any integer from 1 to 2^29, Healpix’s limit, and must be a power of 2 in the NESTED scheme. The mask holds one byte per pixel, 12 NSIDE^2 bytes: 12.6 MB at NSIDE 1024, 805 MB at NSIDE 8192, 3.2 GB at NSIDE 16384. A map too large for memory raises std::bad_alloc.

The object is immutable once built.

Public Functions

explicit Mask(const std::string &fitsFile)
Parameters:

fitsFile – Healpix map: a binary table in HDU 2 with 1024 values to a row, in any numeric column type; NSIDE and ORDERING are taken from the header, in either RING or NESTED scheme. The values are read a block of rows at a time, so reading needs memory for the mask’s bytes only. Each value is tested as cfitsio converts it to double, after any TSCAL and TZERO; a null value of an integer column (TNULL) is unobserved.

Throws:

Error – if the file cannot be opened or read, if it is not a full-sky Healpix map in that layout, if NSIDE exceeds 2^29, or if the ordering is NESTED and NSIDE is not a power of 2.

Mask(const double *values, std::size_t count, PixelOrdering ordering)

Build from the pixel values of a full-sky Healpix map.

NSIDE follows from the number of values, which must be 12 NSIDE^2. The values are read once, during construction, and not kept.

Parameters:

  • values – one value per pixel, in the order of ordering

  • count – the number of values

  • ordering – the pixel ordering of values; no default

Throws:

Error – if count is 0 or not 12 NSIDE^2 for an integer NSIDE, if NSIDE exceeds 2^29, if the ordering is NESTED and NSIDE is not a power of 2, or if values is null; the message gives the number of values and the NSIDE it implies.

Mask(const std::vector<double> &values, PixelOrdering ordering)

As the constructor from a pointer and a count, over all of values.

bool allows(double rightAscension, double declination) const

True when the direction falls in an allowed pixel. Radians.

Throws:

Error – if the right ascension is not finite, or the declination is not in [-pi/2, pi/2].

double skyAreaDeg2() const

Sky area covered by the allowed pixels, in square degrees.

int nside() const

NSIDE of the map, from 1 to 2^29.

enum class otswap::PixelOrdering

Pixel ordering of a full-sky Healpix map given in memory.

Values:

enumerator Ring

the RING scheme

enumerator Nested

the NESTED scheme

Reconstruction

Result otswap::reconstructBox(const std::vector<double> &tracers, const std::vector<double> &randoms, double mps, const Config &config)

Reconstruct in box geometry, with a constant mean particle separation.

Note

The box is not periodic: nothing flows through its faces, so modes on the scale of the box itself cannot be reconstructed.

Note

Given a seed, the result does not depend on the number of threads. The assignment of given randoms to realizations, and in each realization the random draw, the seeding pass and the swap loop, draw from separate streams, each seeded from the run seed, the realization index and the stream alone.

Parameters:

  • tracers – Cartesian coordinates, 3 * nObjects entries.

  • randoms – Cartesian coordinates of the randoms, or empty. When empty, config.nRealizations * nObjects randoms are drawn uniformly in the bounding box of the tracers. When given, the array must hold at least config.nRealizations * nObjects of them, since each realization consumes a disjoint subset; they may lie outside the tracers’ bounding box.

  • mps – mean particle separation; must be positive.

  • config – the parameters of the run; rejectCrossings and maxUnobservedPixelsCrossed are not read.

Throws:

Error – if any array has a size that is not a multiple of three, if the tracer array is empty, if the randoms are too few, if any coordinate is not finite, or if mps is not positive.

Returns:

the displacement field; every entry of valid is 1, validRealizations is uniformly config.nRealizations, and meanDisplacement is the mean over all realizations. Result::mps is mps.

Result otswap::reconstructBox(const std::vector<double> &tracers, const std::vector<double> &randoms, const Config &config)

Reconstruct in box geometry, with the mean particle separation of the tracers: mps = (V / N)^(1/3), V the volume of their bounding box, the box that drawn randoms fill, and N their number. The value used is Result::mps.

Throws:

Error – as above, and if the tracers span no volume (all of them in one plane, on one line or at one point).

Result otswap::reconstructLightcone(const std::vector<double> &tracersSky, const std::vector<double> &randomsSky, double skyAreaDeg2, unsigned nBins, const DistanceTable &distances, const Config &config, const RedshiftCut &cut = {})

Reconstruct in lightcone geometry, from sky coordinates.

The mean particle separation is measured from the tracers themselves: nBins uniform bins in redshift over the observed range, mps = (N / V)^(-1/3) per bin, with V the shell volume implied by skyAreaDeg2 and the distance table. A bin holding fewer than 1/(9 eps^2) tracers, eps = 0.02, is refused: below that count the Poisson error on its mps exceeds eps. The bin values are the nodes of a piecewise linear profile, extrapolated linearly beyond the terminal ones.

Randoms are required and are not generated: they carry the survey geometry, the selection function and the completeness, none of which this package models.

The counts are recorded in Result::selection and, with a cut that has a finite bound, reported as config.verbosity says.

Parameters:

  • tracersSky – sky coordinates, 3 * nObjects entries, radians.

  • randomsSky – sky coordinates of the randoms; at least config.nRealizations * nObjects of them. They may lie outside the Cartesian bounding box of the tracers.

  • skyAreaDeg2 – effective area of the survey. When a mask is available, pass the mask instead, to the overload below, which takes the area from it; otherwise supply the published value: a bounding box in right ascension and declination overestimates it for any footprint that is not rectangular, and the error propagates into the mean particle separation.

  • nBins – redshift bins used to measure mps(z).

  • distances – the redshift to distance relation, with which the sky coordinates are converted and the shell volumes computed; it must cover the redshifts of the objects kept.

  • config – the parameters of the run; rejectCrossings and maxUnobservedPixelsCrossed are read by the mask overloads only.

  • cut – tracers and randoms whose redshift lies outside [cut.min, cut.max] are left out of the reconstruction, before anything else: the randoms are dropped, and the tracers keep their row, flagged in Result::outsideRedshiftCut. mps(z) is measured on the tracers kept, and only their redshifts need lie in the distance table. The default cuts nothing.

Throws:

Error – if a bin holds too few tracers for its density to be meaningful; the message reports how many bins the catalog supports. Also if the profile extrapolates to a mps that is not positive at a tracer’s redshift, if any array is malformed, if a declination lies outside [-pi/2, pi/2], if the randoms are too few, or if a redshift falls outside the distance table. Also if a bound of the cut is NaN or min > max, or if too few tracers or randoms remain inside it; the message gives the counts kept and dropped.

Result otswap::reconstructLightcone(const std::vector<double> &tracersSky, const std::vector<double> &randomsSky, const Mask &mask, unsigned nBins, const DistanceTable &distances, const Config &config, const RedshiftCut &cut = {})

Reconstruct in lightcone geometry, from sky coordinates, within the observed pixels of a mask.

As the overload with the sky area, with the area taken from the mask, mask.skyAreaDeg2(), and the mask applied twice:

  • before the reconstruction, tracers and randoms on unobserved pixels are left out, decided on their right ascension and declination as Mask::allows decides: the randoms are dropped, and the tracers keep their row, flagged in Result::outsideMask;

  • after it, when config.rejectCrossings is set, rejectMaskCrossings is applied with config.maxUnobservedPixelsCrossed, so the result carries filteredNside = mask.nside().

The mask and the redshift cut are both evaluated on every object, so a tracer may be flagged by both; the objects kept are those that pass both. The counts are recorded in Result::selection and reported as config.verbosity says.

To filter with another threshold afterwards, set config.maxUnobservedPixelsCrossed instead: a second rejectMaskCrossings with the same mask and a higher threshold changes nothing, since the filter only clears entries. Or set config.rejectCrossings to false and call rejectMaskCrossings on the result.

Parameters:

  • tracersSky – as in the overload with the sky area.

  • randomsSky – as in the overload with the sky area.

  • nBins – as in the overload with the sky area.

  • distances – as in the overload with the sky area.

  • mask – the survey’s veto mask; borrowed for the call only.

  • config – as in the overload with the sky area; rejectCrossings and maxUnobservedPixelsCrossed are read here.

  • cut – as in the overload with the sky area, evaluated with the mask on every object.

Throws:

Error – as the overload with the sky area, the counts in the message naming the mask as well as the cut.

Result otswap::reconstructLightcone(const std::vector<double> &tracers, const std::vector<double> &randoms, const std::vector<double> &tracersSky, const std::vector<double> &randomsSky, double skyAreaDeg2, unsigned nBins, const DistanceTable &distances, const Config &config, const RedshiftCut &cut = {})

Reconstruct in lightcone geometry, with the Cartesian coordinates already computed.

Identical to the overload above, except that the conversion is skipped. The sky coordinates are still required: the tracers’ redshifts drive mps(z), and the sky random array must describe the same objects as the Cartesian one.

Warning

No check verifies that the Cartesian coordinates agree with the sky ones under the given cosmology; disagreement is silent. Arrays made with toCartesian and the same distance table agree exactly, and give the same result as the sky overload.

Parameters:

  • tracers – Cartesian coordinates, 3 * nObjects entries.

  • randoms – Cartesian coordinates of the randoms.

  • tracersSky – as above.

  • randomsSky – as above.

  • skyAreaDeg2 – as above.

  • nBins – as above.

  • distances – as above.

  • config – as above.

  • cut – as above, decided on the sky redshifts; the same rows are dropped from the Cartesian arrays.

Throws:

Error – as above, and if the Cartesian and sky arrays describe a different number of objects.

Result otswap::reconstructLightcone(const std::vector<double> &tracers, const std::vector<double> &randoms, const std::vector<double> &tracersSky, const std::vector<double> &randomsSky, const Mask &mask, unsigned nBins, const DistanceTable &distances, const Config &config, const RedshiftCut &cut = {})

Reconstruct in lightcone geometry, with the Cartesian coordinates already computed, within the observed pixels of a mask.

The overload above with the conversion skipped, as for the overloads with the sky area. The mask is applied on the sky coordinates; the same rows are dropped from the Cartesian arrays.

Warning

No check verifies that the Cartesian coordinates agree with the sky ones under the given cosmology; disagreement is silent.

Throws:

Error – as above, and if the Cartesian and sky arrays describe a different number of objects.

void otswap::rejectMaskCrossings(Result &result, const Mask &mask, unsigned maxUnobservedPixelsCrossed = 0)

Reject the displacements whose path crosses masked sky.

The segment from the Eulerian to the Lagrangian position lies in the plane through the origin and both endpoints, so its projection on the sky is exactly the great circle arc between them. The arc is bisected adaptively: a sub-arc whose ends fall in the same pixel is not split further, nor one shorter than a millionth of the pixel size, nor one shorter than a pixel whose end pixels share an edge and have no unobserved pixel as a common neighbour. The pixel of every midpoint is examined, and the distinct unobserved pixels met are counted. The pixels holding the two endpoints are exempt.

A displacement is rejected when its arc crosses more than maxUnobservedPixelsCrossed distinct unobserved pixels; the search stops as soon as the count exceeds it. With 0, the first one rejects. A displacement with an endpoint at the origin is kept, and one whose endpoints lie in exactly opposite directions, which defines no arc, is rejected.

A rejected displacement has its entry of valid cleared; entries already cleared stay cleared, so the filter composes with any other and applying it twice with the same mask and threshold changes nothing. The fields that follow from valid are then recomputed, as recomputeMeans does: validRealizations, meanDisplacement, lagrangian and lagrangianSky.

Result::selection records the filter: maxUnobservedPixelsCrossed the threshold of this call, displacements the valid displacements before the first filter applied to the result, and displacementsCrossingMask those rejected by every filter so far.

A tracer flagged in Result::outsideRedshiftCut or Result::outsideMask is left as it is: its displacement and matchedRandom rows may be NaN, and its entries of valid must be 0. Every other displacement must be finite.

Warning

The count is of pixels, not an angular length, so it depends on NSIDE: the same arc crosses about twice as many pixels when NSIDE doubles, and the same threshold is a different criterion.

Parameters:

  • result – a lightcone result, updated in place.

  • mask – the mask whose unobserved pixels are counted.

  • maxUnobservedPixelsCrossed – distinct unobserved pixels tolerated along the arc, the endpoint pixels excluded.

Throws:

Error – if the result is malformed (including a NaN outside the flagged tracers’ rows, a flagged tracer with a valid entry, or a flag, tracers, tracersSky, lagrangian or lagrangianSky array whose size is neither 0 nor the one of its layout), if an entry of valid is neither 0 nor 1, or if the result was already filtered against a mask of a different NSIDE.

void otswap::recomputeMeans(Result &result)

Recompute every field of a result that follows from valid: validRealizations and meanDisplacement, then lagrangian and lagrangianSky.

meanDisplacement is summed over the valid realizations in realization order and divided by their count, as the reconstruction does, so a result whose displacement and valid were read back from a file gives the reconstruction’s bits. lagrangian is tracers + meanDisplacement, left empty when tracers is; lagrangianSky is computed with result.distances, and left empty without a distance table.

Throws:

Error – if displacement or valid has the wrong size, an entry of valid is neither 0 nor 1, or a flag, tracers or tracersSky array has a size that is neither 0 nor the one of its layout.

Results

struct Result

Displacement field produced by a reconstruction.

Displacements point from the Eulerian (observed) position to the Lagrangian one: the Lagrangian position of object i is its coordinate plus its displacement.

Public Members

Geometry geometry = Geometry::Box

The geometry of the reconstruction, which the writers and the redshift-space correction follow.

Config config

The configuration the reconstruction ran with. seed is the one used: the drawn one when 0 was given, so that the result can be reproduced.

std::vector<double> tracers

Cartesian position of each tracer, the start of its displacements: the array given to reconstructBox or to a Cartesian overload of reconstructLightcone, or toCartesian of the sky array. A tracer left out by the redshift cut or the mask has its conversion when its redshift lies in the distance table, NaN otherwise. Flat, [object][xyz]; size 3 * nObjects.

std::vector<double> tracersSky

Sky coordinates of each tracer, as given to reconstructLightcone: right ascension, declination, redshift; radians. Flat, [object][3]; size 3 * nObjects in a lightcone result, empty in a box result.

std::vector<double> displacement

Displacements. Flat, [realization][object][xyz]. Size 3 * nRealizations * nObjects.

std::vector<double> matchedRandom

Cartesian coordinates of the random matched to each object. Same layout and size as displacement.

std::vector<std::uint8_t> valid

Validity of each displacement: 1 valid, 0 rejected by a filter. Flat, [realization][object]. Size nRealizations * nObjects. Every entry is 1 until a filter is applied; a filter only ever clears entries.

std::vector<double> meanDisplacement

Mean displacement per object, over its valid realizations only, summed in realization order, then divided by their count. NaN in all three components for an object with no valid realization. Flat, [object][xyz]. Size 3 * nObjects.

std::vector<unsigned> validRealizations

Valid realizations per object: the count of 1 entries of valid for that object, and the number of terms in its meanDisplacement. Size nObjects.

int filteredNside = 0

NSIDE of the mask this result was last filtered against, or 0 when it has not been filtered. Set by rejectMaskCrossings, which refuses to filter the same result against a mask of a different NSIDE.

std::vector<std::uint8_t> outsideRedshiftCut

Whether each tracer lay outside the redshift cut of reconstructLightcone and so took no part in the reconstruction: 1 for such a tracer, 0 otherwise. Flat, [object]. Size nObjects, all 0 when no cut was given and in a box result.

A cut tracer has NaN in every component of its displacement and matchedRandom rows, valid 0 in every realization, validRealizations 0 and meanDisplacement NaN. No filter changes the flag. A Result built by hand may leave the field empty, which every function reads as all 0.

std::vector<std::uint8_t> outsideMask

Whether each tracer fell on an unobserved pixel of the Mask given to reconstructLightcone, and so took no part in the reconstruction: 1 for such a tracer, 0 otherwise. Flat, [object]. Size nObjects, all 0 when no mask was given and in a box result. Independent of outsideRedshiftCut: a tracer may carry both flags. A flagged tracer’s rows are as for outsideRedshiftCut, and no filter changes the flag. An empty field is read as all 0.

SelectionCounts selection

What reconstructLightcone left out, and why, and what the mask filters rejected; updated by every rejectMaskCrossings.

std::vector<double> lagrangian

Mean Lagrangian position of each tracer, tracers + meanDisplacement. NaN rows where meanDisplacement is NaN. Flat, [object][xyz]; size 3 * nObjects. Recomputed with meanDisplacement by rejectMaskCrossings and recomputeMeans.

std::vector<double> lagrangianSky

Sky coordinates of each tracer’s mean Lagrangian position, as toSky computes them from lagrangian with the reconstruction’s distance table: right ascension in [0, 2 pi), declination, redshift; radians. Flat, [object][3]; size 3 * nObjects in a lightcone result, empty in a box result. Recomputed with lagrangian.

NaN rows where lagrangian is NaN: the tracers flagged in outsideRedshiftCut or outsideMask and those with no valid realization. A position whose distance falls outside the table, which a table that does not start at distance 0 allows, has a NaN redshift and keeps its right ascension and declination.

double mps = std::numeric_limits<double>::quiet_NaN()

Box: the mean particle separation used, given or computed, Mpc/h. NaN in a lightcone result, whose separation is mpsProfile.

MpsProfile mpsProfile

Lightcone: the mean particle separation profile measured on the tracers kept. Empty in a box result.

std::shared_ptr<const DistanceTable> distances

Lightcone: the distance table of the reconstruction, shared by the copies of this result, with which rejectMaskCrossings and recomputeMeans recompute lagrangianSky. Null in a box result.

double elapsedSeconds = 0.

Wall time of the call that made the result, in seconds. Two runs with the same seed give the same result but for this field.

struct MpsProfile

The mean particle separation of a lightcone as a function of redshift, measured by reconstructLightcone.

The tracers kept are binned in nBins uniform redshift bins over [redshiftMin, redshiftMax], and mps = (N / V)^(-1/3) in each, V the shell volume implied by the sky area. The bin values are the nodes of a piecewise linear profile, extrapolated linearly beyond the terminal ones; at() evaluates it. Empty in a box result.

Public Functions

double at(double z) const

The mps at redshift z, in Mpc/h: linear between the nodes, extrapolated linearly beyond the terminal ones along the terminal segment, constant with a single node. This is the value the reconstruction used at each tracer.

Throws:

Error – if the profile is empty, or if the value is not positive, which only an extrapolation can give.

Public Members

std::vector<double> redshift

bin centres

std::vector<double> mps

Mpc/h at each centre.

std::vector<std::size_t> count

tracers in each bin

double redshiftMin = std::numeric_limits<double>::quiet_NaN()

The binned range: the smallest and largest redshift of the tracers kept.

double representative = std::numeric_limits<double>::quiet_NaN()

The count-weighted median of the nodes, which sizes the grids of the reconstruction.

struct SelectionCounts

The objects a lightcone reconstruction left out, and why, and the displacements its mask filter rejected.

The tracers removed are tracersOutsideRedshiftCut + tracersOutsideMask

  • tracersOutsideBoth, and the same for the randoms. All counts are 0, and nothing is applied, in a box result.

Public Functions

std::string message() const

The selection report of a reconstruction.

One line for the tracers and one for the randoms, each listing the selections applied, then one for the rejected crossings when a filter was applied; each line ends in ‘

’. A line is written even when nothing was removed. In the form:

otswap: kept <n> of <N> tracers: <a> outside the redshift cut [<zmin>, <zmax>], <b> outside the mask (<c> outside both)
otswap: kept <m> of <M> randoms: <a> outside the redshift cut [<zmin>, <zmax>], <b> outside the mask (<c> outside both)
otswap: rejected <r> of <D> displacements crossing more than <L> unobserved pixels

Returns:

the report, or an empty string when no selection was applied.

Public Members

std::optional<RedshiftCut> redshiftCut

the cut, when a bound was finite

bool maskApplied = false

a Mask was given

std::size_t tracers = 0

tracers given

std::size_t tracersOutsideRedshiftCut = 0

of them, outside the cut

std::size_t tracersOutsideMask = 0

of them, on an unobserved pixel

std::size_t tracersOutsideBoth = 0

counted in both of the above

std::size_t randoms = 0

randoms given

std::size_t randomsOutsideRedshiftCut = 0

of them, outside the cut

std::size_t randomsOutsideMask = 0

of them, on an unobserved pixel

std::size_t randomsOutsideBoth = 0

counted in both of the above

std::optional<unsigned> maxUnobservedPixelsCrossed

The threshold of the latest rejectMaskCrossings applied to the result, by reconstructLightcone or afterwards; empty when none was.

std::size_t displacements = 0

valid displacements before the first filter

std::size_t displacementsCrossingMask = 0

of them, rejected by every filter so far

otswap/RSD.h

The chain

RealSpaceCatalog otswap::realSpaceLightcone(const Result &result, const DistanceTable &distances, const BiasTable &bias, double sigma, const CorrectionConfig &config = {})

Move a lightcone catalogue from redshift space to real space.

The chain, for the tracers of the result:

  1. the projection of result.meanDisplacement on the line of sight of each tracer, at result.tracers (radialProjection);

  2. its gaussian average of width sigma over the tracers with a valid realization (neighbourAverage); a tracer without one still receives the average around its own position;

  3. the factor f / (b + 3 f / 5) at each observed redshift, as rsdFactor computes it; the tracers at which b(z) is extrapolated are counted in RealSpaceCatalog::nExtrapolated, and reported;

  4. each tracer keeps its right ascension and declination and moves to the redshift of comoving distance d(z_i) + s_i, where s_i is the factor times the average; its Cartesian position moves by s_i along its line of sight.

Tracers flagged in result.outsideRedshiftCut or result.outsideMask are left out: they take no part in any average, are not moved, have diagnostics 0, and the status LeftOut.

Parameters:

  • result – a lightcone reconstruction, with its tracers and tracersSky, as reconstructLightcone returns it.

  • distances – a table carrying the growth rate, covering the corrected distances as well as the observed ones.

  • bias – the b(z) table.

  • sigma – width of the average, in Mpc/h; 0 for none. 10 Mpc/h is a reasonable starting value; the best value depends on the sample, and should be checked in each analysis.

  • config – the weighting, and what the call reports.

Throws:

  • Error – if result is not a lightcone result, is malformed or lacks its tracers or tracersSky, if a tracer kept is not finite or is at the origin, if the table carries no growth rate, if the bias table is malformed or extrapolates to a bias that is not positive, or if a corrected comoving distance is not positive or lies outside the distance table; the message names the object, and in the last case suggests a table covering a wider redshift range.

  • meshsearch::Error – as neighbourAverage.

RealSpaceCatalog otswap::realSpaceBox(const Result &result, unsigned axis, double redshift, const DistanceTable &distances, double bias, double sigma, const CorrectionConfig &config = {})

Move a box catalogue from redshift space to real space, with the line of sight along a Cartesian axis.

The chain is that of realSpaceLightcone, on result.tracers, with the projection on the axis, a single factor f / (b + 3 f / 5) at the box’s redshift, and the shift applied along the axis. Positions are not wrapped.

The box is not periodic, here or in reconstructBox: nothing flows through its faces, so modes on the scale of the box itself cannot be reconstructed.

Parameters:

  • result – a box reconstruction, with its tracers.

  • axis – the line of sight: 0, 1 or 2.

  • redshift – the box’s redshift, at which f is taken.

  • distances – a table carrying the growth rate at that redshift.

  • bias – the tracers’ linear bias; finite and positive.

  • sigma – as in realSpaceLightcone.

  • config – the weighting, and what the call reports.

Throws:

  • Error – if result is not a box result, is malformed or lacks its tracers, if the axis is not 0, 1 or 2, or as rsdFactorBox.

  • meshsearch::Error – as neighbourAverage.

struct CorrectionConfig

Options of realSpaceLightcone and realSpaceBox.

Public Members

bool weightByRealizations = false

Weight each neighbour by its count of valid realizations as well, as in neighbourAverage.

Verbosity verbosity = Verbosity::Normal

What the call reports; see Verbosity.

std::ostream *log = &std::clog

Where the report is written; null writes nothing, as Silent.

struct RealSpaceCatalog

A catalogue moved to real space. Row i describes input object i.

Public Members

Geometry geometry = Geometry::Box

The geometry of the reconstruction corrected.

std::vector<double> sky

Lightcone: right ascension, folded into [0, 2 pi), and declination of the tracer, and its corrected redshift; radians. Flat, [object][3]. NaN rows for the tracers not moved. Empty in a box.

std::vector<double> cartesian

The corrected Cartesian position, Mpc/h. Box: the tracer with its axis coordinate moved. Lightcone: r (1 + shift / |r|), the tracer r moved along its line of sight. Flat, [object][xyz]. NaN rows for the tracers not moved.

std::vector<double> shift

The shift applied along the line of sight, in Mpc/h: the factor times the averaged projection. NaN for a tracer not moved. nObjects entries.

std::vector<double> factor

The factor f / (b + 3 f / 5) at each tracer’s observed redshift; in a box, the single factor repeated. NaN for a tracer left out of the reconstruction. nObjects entries.

std::vector<CorrectionStatus> status

What the correction did with each tracer. nObjects entries.

std::vector<unsigned> validRealizations

The valid realizations of each tracer, copied from the result corrected. nObjects entries.

std::vector<unsigned> nNeighbours

The diagnostics of neighbourAverage: valid tracers averaged, and the sum of their valid realizations; 0 for a tracer left out. nObjects entries each.

std::size_t nExtrapolated = 0

Tracers at whose redshift b(z) was extrapolated beyond its table; 0 in a box.

double sigma = std::numeric_limits<double>::quiet_NaN()

The width of the average and its weighting, as given.

unsigned axis = 0

Box: the line-of-sight axis, the redshift and the bias, as given. 0 and NaN in a lightcone.

double elapsedSeconds = 0.

Wall time of the call that made the catalogue, in seconds.

enum class otswap::CorrectionStatus : std::uint8_t

What the correction did with each tracer.

Values:

enumerator Corrected

had a valid realization; moved

enumerator MovedByNeighbours

had no valid realization; moved with its neighbours’ average

enumerator NoValidNeighbour

no valid tracer within 3 sigma (with sigma = 0: no valid realization); not moved

enumerator LeftOut

outside the reconstruction’s redshift cut or mask; not moved

The four steps

std::vector<double> otswap::radialProjection(const std::vector<double> &displacement, const std::vector<double> &positions)

Component of each displacement along the line of sight of its object, r_i . d_i / |r_i|: positive when the displacement points away from the observer.

Parameters:

  • displacement – 3 * N, for example Result::meanDisplacement. NaN is allowed, and gives NaN for that object; infinities are not.

  • positions – Cartesian, 3 * N, for example Result::tracers; every position finite and away from the origin.

Throws:

Error – if the sizes differ or are not multiples of three, if a position is not finite or is at the origin, or if a displacement is infinite; the message names the object.

std::vector<double> otswap::axisProjection(const std::vector<double> &displacement, unsigned axis)

Component of each displacement along a Cartesian axis, the line of sight of a box.

Parameters:

  • displacement – 3 * N; NaN allowed, infinities not.

  • axis – 0, 1 or 2.

Throws:

Error – if the size is not a multiple of three, if the axis is not 0, 1 or 2, or if a displacement is infinite.

NeighbourAverage otswap::neighbourAverage(const std::vector<double> &positions, const std::vector<double> &values, const std::vector<unsigned> &validRealizations, double sigma, bool weightByRealizations = false)

Gaussian average of a value over the neighbours of each object.

An object is valid when its validRealizations is greater than zero. Each object, valid or not, receives the weighted mean of the values of the valid objects within 3 sigma of its own position, its own value included when it is valid. The weight of a neighbour at distance d is exp(-d^2 / (2 sigma^2)), multiplied by its validRealizations when weightByRealizations is set, and the weights are normalised on their sum, so a constant value is returned unchanged, at the edges of the catalogue too. An object with no valid object within 3 sigma receives NaN.

With sigma = 0 no average is taken and values is returned unchanged.

The neighbours are found on a grid built inside, over the bounding box of the positions given; nothing passed in is modified. The result does not depend on that grid within rounding, and is the same bit for bit for the same inputs, whatever the number of threads.

Parameters:

  • positions – Cartesian, 3 * N, finite.

  • values – N entries; finite wherever validRealizations > 0, and not read elsewhere.

  • validRealizations – N entries.

  • sigma – width of the gaussian, in the unit of the positions; finite and non-negative. 10 Mpc/h is a reasonable starting value; the best value depends on the sample, and should be checked in each analysis.

  • weightByRealizations – weight each neighbour by its count of valid realizations as well.

Throws:

  • Error – if the sizes differ, if a position is not finite, if the value of a valid object is not finite, or if sigma is negative or not finite; the message names the object.

  • meshsearch::Error – if sigma > 0 and the positions span no volume (a single object, or all of them in one plane or on one line): the neighbour grid cannot be built.

Returns:

the average and two diagnostics per object.

struct NeighbourAverage

The result of neighbourAverage: N entries each.

Public Members

std::vector<double> values

The average at each object; NaN with no valid object within 3 sigma.

std::vector<unsigned> nNeighbours

The number of valid objects averaged, the object itself included when it is valid. With sigma = 0, 1 for a valid object and 0 otherwise.

std::vector<unsigned> nRealizationsAveraged

The sum of validRealizations over those objects.

std::vector<double> otswap::rsdFactor(const std::vector<double> &redshift, const DistanceTable &distances, const BiasTable &bias, std::size_t *nExtrapolated = nullptr)

The factor f / (b + 3 f / 5) that turns the averaged line-of-sight displacement into the shift to real space, at each redshift.

f is distances.growthRateAt(z), b(z) the bias table at z. Nothing is written: the redshifts at which b(z) is extrapolated beyond the table are counted in *nExtrapolated, when given.

Parameters:

  • redshift – observed redshifts, inside the distance table’s range.

  • distances – a table carrying the growth rate.

  • bias – the b(z) table.

  • nExtrapolated – if not null, receives the number of redshifts outside the bias table.

Throws:

Error – if the distance table carries no growth rate, if a redshift is not finite or lies outside the distance table, if the bias table is malformed, or if the extrapolated bias at a redshift is not positive; the message names the object.

double otswap::rsdFactorBox(double redshift, const DistanceTable &distances, double bias)

The factor f / (b + 3 f / 5) of a box at a single redshift, f = distances.growthRateAt(redshift), with a constant bias.

Parameters:

  • redshift – the box’s redshift, inside the distance table’s range.

  • distances – a table carrying the growth rate.

  • bias – finite and positive.

Throws:

Error – if the redshift is not finite or lies outside the distance table, if the table carries no growth rate, or if the bias is not positive.

struct BiasTable

A b(z) table: the linear bias of the tracers at redshift nodes, interpolated linearly between them and extrapolated linearly beyond them from the first or last segment. Read one from a file with io::readBiasTable.

Public Members

std::vector<double> redshift

at least two, finite, strictly increasing

std::vector<double> bias

same length, finite and positive

std::vector<double> otswap::shiftRadially(const std::vector<double> &positions, const std::vector<double> &shift)

Move each position by shift_i along its own line of sight: r_i (1 + shift_i / |r_i|). A positive shift moves the object away from the observer.

Parameters:

  • positions – Cartesian, 3 * N; every position finite and away from the origin.

  • shift – N entries; NaN gives NaN in the whole row, infinities are refused.

Throws:

Error – if the sizes differ, if a position is not finite or is at the origin, or if a shift is infinite; the message names the object.

std::vector<double> otswap::shiftAlongAxis(const std::vector<double> &positions, const std::vector<double> &shift, unsigned axis)

Move each position by shift_i along a Cartesian axis. Positions are not wrapped: a box is not treated as periodic.

Parameters:

  • positions – Cartesian, 3 * N, finite.

  • shift – N entries; NaN gives NaN along the axis, infinities are refused.

  • axis – 0, 1 or 2.

otswap/io.h

Tables

struct Column

One column of a table to write.

‘D’ (double) and ‘J’ (32-bit integer) columns hold their values in data; ‘K’ (64-bit integer) columns in integers, exactly. The other vector is ignored.

Public Members

char type = 'D'

‘D’, ‘J’ or ‘K’

std::string description

ASCII header comment; FITS TCOMMn.

std::string unit

ASCII “[unit]” after the description; FITS TUNITn.

std::vector<double> data

the values of a ‘D’ or ‘J’ column

std::vector<std::int64_t> integers

the values of a ‘K’ column

struct Table

A table read from file.

Public Members

std::size_t nColumns = 0

the columns read as double

std::vector<double> values

flat, row-major [row][column]

std::size_t nIntegerColumns = 0

the columns read as 64-bit integers

std::vector<std::int64_t> integers

flat, row-major [row][integer column]

struct Keyword

A header keyword: a FITS keyword in the header of HDU 2, an ASCII “## NAME = value / comment” line. The name is at most 8 characters of A-Z, 0-9, ‘-’ and ‘_’. In FITS a value that reads entirely as an integer or a finite number is written as one, anything else as a string.

struct WriteOptions

Options of io::write.

Public Members

int precision = 9

Significant digits of a ‘D’ value in ASCII, from 1 to 17; 17 reads back as the same double. FITS stores the doubles themselves.

std::vector<Keyword> keywords

Keywords written in the header, in this order.

Table otswap::io::read(const std::string &file, const std::vector<std::string> &columns, char delimiter = ' ', char comment = '#')

Read selected columns of a table, as double.

In FITS the columns are named, case insensitively, and must be scalar. In ASCII a column given by digits alone is a 0-based index; any other string is a name, looked up in the last line starting with “###” before the first data row, which is where io::write puts the names. Lines that are empty or start with comment are skipped, and fields are separated by delimiter, blanks or tabs.

NaN and infinities are returned as stored, in both formats: a file io::write wrote with NaN rows reads back with them. A null value of a FITS integer column (TNULL) is an error.

Throws:

Error – if the file is missing or unreadable, a column is absent (the message lists the names found), a column is requested twice, a requested index is too large, a row is too short, or a null integer is met.

Table otswap::io::read(const std::string &file, const std::vector<std::string> &columns, const std::vector<std::string> &integerColumns, char delimiter = ' ', char comment = '#')

Read selected columns as double, and others exactly as 64-bit integers.

columns are returned in Table::values, integerColumns in Table::integers, each in the order requested. In ASCII an integer field is an optional ‘-’ followed by digits; in FITS an integer column has an integer type (B, I, J or K).

Throws:

Error – as above, and if an integer field holds anything else (a decimal point, an exponent, a ‘+’, nan), lies outside the 64-bit range, is null, or comes from a FITS column of floating type; the message names the row and the column.

void otswap::io::write(const std::string &file, const std::vector<Column> &columns, const WriteOptions &options = {})

Write a table. An existing file is overwritten.

Throws:

Error – if no column is given, a column type is not ‘D’, ‘J’ or ‘K’, the columns differ in length (data for ‘D’ and ‘J’, integers for ‘K’), a ‘J’ value is NaN, not an integer, or outside [-2^31, 2^31) (the message names the column and the row), a keyword name is invalid, the precision is outside [1, 17], or the file cannot be written.

void otswap::io::write(const std::string &file, const std::vector<Column> &columns, std::size_t nRows, const std::function<void(std::size_t, std::vector<double>&)> &fillRow, const WriteOptions &options = {})

Write a table whose rows are produced on demand.

The data and integers of each column are ignored; fillRow is called once per row, in order, with a buffer of one entry per column. Nothing is held in memory, so this is the variant for catalogue-sized output. ‘K’ columns need the overload below.

Throws:

Error – as above, and if a column is of type ‘K’.

void otswap::io::write(const std::string &file, const std::vector<Column> &columns, std::size_t nRows, const std::function<void(std::size_t, std::vector<double>&, std::vector<std::int64_t>&)> &fillRow, const WriteOptions &options = {})

As above, with ‘K’ columns: fillRow receives two buffers of one entry per column, and each column reads its own entry of the buffer of its type, values for ‘D’ and ‘J’, integers for ‘K’; the entry of the other buffer is ignored.

BiasTable otswap::io::readBiasTable(const std::string &file, char delimiter = ' ', char comment = '#')

Read a b(z) table.

The format follows the extension: FITS from the columns named REDSHIFT and BIAS (case insensitively); ASCII from columns 0 and 1, with the given delimiter and comment character.

The table is checked as rsdFactor checks a BiasTable, so a table read without error is one rsdFactor accepts.

Throws:

Error – if the file cannot be read (as io::read), if it holds fewer than two rows, if a redshift is not finite or not greater than the one before, or if a bias is not finite or not positive; the message names the file and the row.

Writers

io::writeDisplacements, io::writeDisplacementField, io::writeRealSpaceCatalog, io::writeMpsProfile and their group enums are documented on the File formats page, with the layouts they write.