DMTN-043: Redesign of afw::math::Statistics

This document captures the redesign work of the lsst::afw::math::Statistics functionality.

Warning

This is an initial draft. A partially working proof of concept implementation,
which sketches out the general design and has been benchmarked to be equivalent
to the current Statistics performance on 4k x 4k noise images, is available at:

Where img, msk, wgt and var are typically a 1 dimensional numpy.ndarray, or any of the other supported overloads (see below).

In C++ the signature is:

template<typenameImageT,typenameMaskT,typenameWeightT,typenameVarianceT>StatisticsResultcomputeStatistics(ImageTconst&img,MaskTconst*msk,// (nullptr if not used)WeightTconst*wgt,// (nullptr if not used)VarianceTconst*var,// (nullptr if not used)intconstflags,StatisticsControlconst&sctrl=StatisticsControl());

Where img, msk, wgt and var are typically ndarray::Array<T,1,1> of the appropriate type T, and where flags is a bitwise combination of:

Statistics to compute are selected either with keyword arguments (Python) or a bitmask (C++).
Additional configuration settings are provided through an (optional) StatisticsControl object
using standard functionality provided by pex_config.

classNewStatisticsControl{public:NewStatisticsControl():numSigmaClip(3),numIter(3),andMask(0x0),noGoodPixelsMask(0x0),isNanSafe(true),calcErrorFromInputVariance(true),baseCaseSize(100),maskPropagationThresholds(16){}LSST_CONTROL_FIELD(numSigmaClip,double,"Number of standard deviations to clip at");LSST_CONTROL_FIELD(numIter,int,"Number of iterations");LSST_CONTROL_FIELD(andMask,typenameimage::MaskPixel,"and-Mask to specify which mask planes to ignore");LSST_CONTROL_FIELD(noGoodPixelsMask,typenameimage::MaskPixel,"mask to set if no values are acceptable");LSST_CONTROL_FIELD(isNanSafe,bool,"Check for NaNs & Infs before running (slower)");LSST_CONTROL_FIELD(calcErrorFromInputVariance,bool,"Calculate errors from the input variances, if available");LSST_CONTROL_FIELD(baseCaseSize,int,"Size of base case in partial sum for numerical stability");LSST_CONTROL_FIELD(maskPropagationThresholds,std::vector<double>,"Thresholds for when to propagate mask bits, treated like a dict (unset bits are set to 1.0)");};

The general design uses a Strategy like structure but with compile-time (instead of run-time) polymorphism.
It has the following components:

Algorithm: calculates (collect) and combines (combineWith) intermediate values and produces the final StatisticsResult (reduce).
Relies on external state (ExternalData) for nested runs.
Only StandardStatistics is implemented, but can be swapped out later.

Runner<Algorithm,Validator>: performs (operator()) a single pass run of the algorithm on the data.
Holds a (unique) instance of ExternalData.
Only SinglePassStatistics is implemented, but can be swapped out later.

Validator: checks if pixels should be included or excluded.
Are composable with various simple building blocks available.

The externally visible computeStatistics function calls an internal detail::translateOptions,
which has a series of overloads that reduce runtime function arguments to compile-time template parameters.
At the end of this series this results in a call to detail::computeStatistics.

The default runner is SinglePassStatistics. This implementation uses pairwise summation (combination)
down to a baseCaseSize after which it performs a simple linear run through the data.
Internally it keeps one instance of Algorithm::ExternalData which is passed to the Algorithm::collect
and Algorithm::reduce functions (by reference) where needed.

It takes the following template parameters:

Validator: A functor that takes image, mask, weight and variance arguments and returns true
if the value is to be clipped, or false otherwise.

Algorithm: The statistics algorithm to execute.

Constructors:

explicitSinglePassStatistics(constValidator&validator=Validator(),size_tbaseCaseSize=100)
Copies of both input parameters are held as private members.

default copy and move constructors (as well as copy and move assignment operators).

Member functions and operators:

operator(): Run Algorithm once on provided img, msk, wgt and var
using pairwise summation / combination.

During a run of Algorithm::collect each set of img, msk, wgt and var values
is passed to the provided Validator functor. If its return value is false the pixel is clipped.

Validators are composable. The makeCombinedChecker variadic function creates a single validator out of its provided arguments.

It does this by nesting instances of CheckBoth which itself takes two template arguments First and Second and has an operator() that returns the logical and of First::operator() and Second::operator().

The following building blocks are provided:

AlwaysTrue: always returns true;

AlwaysFalse: always returns false;

CheckMask: returns bitwise and of msk value with mask argument provided at construction;

CheckFinite: returns true if img is finite (e.g. not NaN or infinity);

CheckRange: returns true if img is within the interval center+/-limit where
center and limit are provided at construction.

Unused parameters (for msk, wgt or var) are given as nullptr (C++) or None (Python) in the API, but for lower layers are given as instances of UnusedParameter.
A fake vector / iterator that does nothing.
This makes the inner loop in Algorithm::collect straightforward while compiling away into non-existence..