Demo of global <-> local bin coordinate conversions

demo.cpp

#include <algorithm>
#include <array>
#include <exception>
#include <iostream>
#include <utility>
#include <vector>


#define ASSERT(x, msg) if (!(x)) throw std::runtime_error(msg)


struct IAxis {
    virtual bool CanGrow() const = 0;
    virtual int GetNBinsNoOver() const = 0;
    virtual int GetNOverflowBins() const = 0;
    int GetNBins() const {
        return GetNBinsNoOver() + GetNOverflowBins();
    }
    virtual int FindBin(double x) const = 0;
    virtual double GetBinFrom(int bin) const = 0;
    virtual double GetBinTo(int bin) const = 0;
};


class EqAxisMock : public IAxis
{
public:
    EqAxisMock(double from, double to, int nBinsNoOver, bool canGrow) :
        m_from(from), m_to(to), m_nBinsNoOver(nBinsNoOver), m_canGrow(canGrow)
    {
        ASSERT(from < to, "Axis limits must be provided in the right order");
        ASSERT(nBinsNoOver > 0, "Invalid number of bins");
    }

    bool CanGrow() const final override {
        return m_canGrow;
    }

    int GetNBinsNoOver() const final override {
        return m_nBinsNoOver;
    }

    int GetNOverflowBins() const final override {
        return m_canGrow ? 0 : 2;
    }

    int FindBin(double x) const final override {
        if (x < m_from) {
            return m_canGrow ? 0 : -1;
        } else if (x >= m_to) {
            return m_canGrow ? 0 : -2;
        } else {
            return (x - m_from) / GetBinWidth() + 1;
        }
    }

    double GetBinFrom(int bin) const final override {
        switch (bin) {
            case -1:
                ASSERT(!m_canGrow,
                       "Growable axes don't have an underflow bin");
                return -std::numeric_limits<double>::infinity();
            case -2:
                ASSERT(!m_canGrow,
                       "Growable axes don't have an overflow bin");
                return m_to;
            default:
                ASSERT(bin != 0, "0 is not a valid local bin index");
                ASSERT(bin <= GetNBinsNoOver(),
                       "Local bin index is out of axis range");
                return m_from + (bin - 1) * GetBinWidth();
        }
    }

    double GetBinTo(int bin) const final override {
        switch (bin) {
            case -1:
                ASSERT(!m_canGrow,
                       "Growable axes don't have an underflow bin");
                return m_from;
            case -2:
                ASSERT(!m_canGrow,
                       "Growable axes don't have an overflow bin");
                return std::numeric_limits<double>::infinity();
            default:
                ASSERT(bin != 0, "0 is not a valid local bin index");
                ASSERT(bin <= GetNBinsNoOver(),
                       "Local bin index is out of axis range");
                return m_from + bin * GetBinWidth();
        }
    }

private:
    double m_from, m_to;
    int m_nBinsNoOver;
    bool m_canGrow;

    // NOTE: Specific to equidistant axis, so general algorithm shouldn't use it
    double GetBinWidth() const {
        return (m_to - m_from) / m_nBinsNoOver;
    }
};


// ---


// Apply a function to each element of an array and return the array of results
//
// You'll need to specify the "Out" type because C++ is bad at type inference.
//
template <typename Out, typename In, std::size_t NDIMS, typename UnaryFunction>
std::array<Out, NDIMS> Map(const std::array<In, NDIMS>& inputs,
                           UnaryFunction&& f) {
    std::array<Out, NDIMS> outputs;
    for (std::size_t dim = 0; dim < NDIMS; ++dim) {
        outputs[dim] = f(inputs[dim]);
    }
    return outputs;
}


// Shorthand for calling GetNBins on an array of axes
template <std::size_t NDIMS>
std::array<int, NDIMS> GetNBins(const std::array<const IAxis*, NDIMS>& axes) {
    return Map<int>(axes,
                    [](const IAxis* axis) { return axis->GetNBins(); });
}


// Shorthand for calling GetNBinsNoOver on an array of axes
template <std::size_t NDIMS>
std::array<int, NDIMS> GetNBinsNoOver(const std::array<const IAxis*, NDIMS>& axes) {
    return Map<int>(axes,
                    [](const IAxis* axis) { return axis->GetNBinsNoOver(); });
}


// Given an array of inputs and an array of axis pointers of the same size,
// apply a binary function to each (input, axis) pair and return the results.
//
// A reader experienced with functional programming may recognize a Zip->Map
// graph and wonder why I'm not just implementing Zip. The answer is that C++'s
// poor type inference, limited type system vocabulary and super-verbose lambda
// syntax make any nontrivial combination of functional programming constructs
// look like an unreadable mess.
//
// Specialized combinators ask less from the broken parts of C++, and are thus
// preferrable over their general counterparts when you don't need generality.
//
template <typename Out, typename In, std::size_t NDIMS, typename BinaryFunction>
std::array<Out, NDIMS> ZipAxisMap(const std::array<In, NDIMS>& inputs,
                                  const std::array<const IAxis*, NDIMS>& axes,
                                  BinaryFunction&& f) {
    std::array<Out, NDIMS> outputs;
    for (std::size_t dim = 0; dim < NDIMS; ++dim) {
        outputs[dim] = f(inputs[dim], *axes[dim]);
    }
    return outputs;
}


// Shorthand for calling GetBinFrom from an array of axes to an array of bins
template <std::size_t NDIMS>
std::array<double, NDIMS> GetBinFrom(const std::array<int, NDIMS>& bins,
                                     const std::array<const IAxis*, NDIMS>& axes) {
    return ZipAxisMap<double>(
        bins,
        axes,
        [](int bin, const IAxis& axis) { return axis.GetBinFrom(bin); }
    );
}


// Shorthand for calling GetBinTo from an array of axes to an array of bins
template <std::size_t NDIMS>
std::array<double, NDIMS> GetBinTo(const std::array<int, NDIMS>& bins,
                                   const std::array<const IAxis*, NDIMS>& axes) {
    return ZipAxisMap<double>(
        bins,
        axes,
        [](int bin, const IAxis& axis) { return axis.GetBinTo(bin); }
    );
}


// Compute a zero-based global bin index given...
//
// - A set of zero-based per-axis bin indices
// - The number of considered bins on each axis (can be either GetNBinsNoOver
//   or GetNBins depending on what you are trying to do)
// - A policy of treating all bins are regular (i.e. no negative indices)
//
template <std::size_t NDIMS>
int ComputeGlobalBinRaw(const std::array<int, NDIMS>& zero_based_bins,
                        const std::array<int, NDIMS>& axis_bin_counts) {
    int result = 0;
    int bin_size = 1;
    for (std::size_t dim = 0; dim < NDIMS; ++dim) {
        ASSERT(zero_based_bins[dim] >= 0,
               "Expected a zero-based local bin index, got a negative one");
        ASSERT(axis_bin_counts[dim] >= 1,
               "Expected a positive bin count, got a negative or null one");
        ASSERT(zero_based_bins[dim] < axis_bin_counts[dim],
               "Input local bin index is outside of input axis range");
        result += zero_based_bins[dim] * bin_size;
        bin_size *= axis_bin_counts[dim];
    }
    return result;
}


// Compute zero-based local bin indices given...
//
// - A zero-based global bin index
// - The number of considered bins on each axis (can be either GetNBinsNoOver
//   or GetNBins depending on what you are trying to do)
// - A policy of treating all bins are regular (i.e. no negative indices)
//
template <std::size_t NDIMS>
std::array<int, NDIMS> ComputeLocalBinsRaw(
    int global_bin,
    const std::array<int, NDIMS>& axis_bin_counts
) {
    ASSERT(global_bin >= 0,
           "Expected a zero-based global bin index, got a negative one");
    std::array<int, NDIMS> result;
    for (std::size_t dim = 0; dim < NDIMS; ++dim) {
        ASSERT(axis_bin_counts[dim] >= 1,
               "Expected a positive bin count, got a negative or null one");
        result[dim] = global_bin % axis_bin_counts[dim];
        global_bin /= axis_bin_counts[dim];
    }
    ASSERT(global_bin == 0, "Input global bin index is out of range");
    return result;
}


// Convert a local axis bins from the standard -1/-2 under/overflow bin indexing
// convention to a "virtual bin" convention where the underflow bin has index 0
// and the overflow bin has index N+1 where N is the axis' regular bin count.
//
// For growable axes, subtract 1 from regular indices so that the indexing
// convention remains zero-based (this means that there will be no "holes" in
// global binning, which matters more than the choice of regular index base)
//
int LocalBinToVirtualBin(int bin, const IAxis& axis) {
    switch (bin) {
        case -1:
            ASSERT(!axis.CanGrow(), "Growable axes shouldn't underflow");
            return 0;
        case -2:
            ASSERT(!axis.CanGrow(), "Growable axes shouldn't overflow");
            return axis.GetNBins() - 1;
        default:
            ASSERT(bin != 0, "0 is not a valid local bin index");
            ASSERT(bin <= axis.GetNBinsNoOver(),
                   "Input local bin index is out of axis range");
            return bin - static_cast<int>(axis.CanGrow());
    }
}


// Convert back from zero-based virtual bins to the standard under/overflow bin
// indexing convention.
//
// For growable axes, add 1 in order to go back to the usual 1-based regular bin
// indexing convention, thus reversing the effect of LocalBinToVirtualBin.
//
int VirtualBinToLocalBin(int virtual_bin, const IAxis& axis) {
    if ((!axis.CanGrow()) && (virtual_bin == 0)) {
        return -1;
    } else if ((!axis.CanGrow()) && (virtual_bin == (axis.GetNBins() - 1))) {
        return -2;
    } else {
        const int regular_bin_offset = -static_cast<int>(axis.CanGrow());
        ASSERT(virtual_bin >= 1 + regular_bin_offset,
               "Expected a zero-based virtual bin index, got a negative one");
        ASSERT(virtual_bin <= axis.GetNBinsNoOver() + regular_bin_offset,
               "Input virtual bin index is out of axis range");
        return virtual_bin - regular_bin_offset;
    }
}


// Compute the global index of a certain bin on an N-dimensional histogram,
// knowing the local bin indices as returned by calling FindBin on each axis.
template <std::size_t NDIMS>
int ComputeGlobalBin(const std::array<int, NDIMS>&          bins,
                     const std::array<const IAxis*, NDIMS>& axes) {
    // Get regular bins out of the way
    if (std::all_of(bins.cbegin(), bins.cend(),
                    [](int bin) { return bin >= 1; }))
    {
        return ComputeGlobalBinRaw(
            Map<int>(bins, [](int bin) { return bin - 1; }),
            GetNBinsNoOver(axes)
        ) + 1;
    }

    // Convert bin indices to a zero-based coordinate system where the underflow
    // bin (if any) has coordinate 0 and the overflow bin (if any) has
    // coordinate N-1, where N is the axis' total number of bins.
    std::array<int, NDIMS> virtual_bins =
        ZipAxisMap<int>(bins, axes, LocalBinToVirtualBin);

    // Deduce what the global bin index would be in this coordinate system that
    // unifies regular and overflow bins.
    const int global_virtual_bin = ComputeGlobalBinRaw(virtual_bins,
                                                       GetNBins(axes));

    // Move to 1-based and negative indexing
    const int neg_1based_virtual_bin = -global_virtual_bin - 1;

    // At this point, we have an index that represents a count of all bins, both
    // regular and overflow, that are located before the current bin when
    // enumerating histogram bins in row-major order.
    //
    // We will next count the number of _regular_ bins which are located before
    // the current bin, and by removing this offset from the above index, we
    // will get a count of overflow bins that are located before the current bin
    // in row-major order. Which is what we want as our overflow bin index.
    //
    int total_regular_bins_before = 0;

    // First, we need to know how many regular bins we leave behind us for each
    // step on each axis, assuming that the bin from which we come was regular.
    //
    // If mathematically inclined, you can also think of this as the size of an
    // hyperplane of regular bins when projecting on lower-numbered dimensions.
    // See the docs of ComputeLocalBins for more on this worldview.
    //
    std::array<int, NDIMS> bin_sizes;
    bin_sizes[0] = 1;
    for (std::size_t dim = 1; dim < NDIMS; ++dim) {
        bin_sizes[dim] = bin_sizes[dim-1] * axes[dim-1]->GetNBinsNoOver();
    }

    // Then, starting from the _last_ histogram dimension...
    for (int dim = NDIMS-1; dim >= 0; --dim) {
        // We can tell how many regular bins lie before us on this axis,
        // accounting for the underflow bin of this axis if it has one.
        const int num_underflow_bins = static_cast<int>(!axes[dim]->CanGrow());
        const int num_regular_bins_before =
            std::max(virtual_bins[dim] - num_underflow_bins, 0);
        total_regular_bins_before += num_regular_bins_before * bin_sizes[dim];

        // If we are on an overflow or underflow bin on this axis, we know that
        // we don't need to look at the remaining axes. Projecting on those
        // dimensions would only take us into an hyperplane of over/underflow
        // bins for the current axis, and we know that by construction there
        // will be no regular bins in there.
        if (bins[dim] < 1) break;
    }

    // Now that we know how many bins lie before us, and how many of those are
    // regular bins, we can trivially deduce how many overflow bins lie before
    // us, and emit that as our global overflow bin index.
    return neg_1based_virtual_bin + total_regular_bins_before;
}


// Given a global histogram bin index as generated by ComputeGlobalBin above,
// go back to per-axis "local" bin coordinates.
template<std::size_t NDIMS>
std::array<int, NDIMS> ComputeLocalBins(
    int global_bin,
    const std::array<const IAxis*, NDIMS> axes
) {
    // Get regular bins out of the way
    if (global_bin >= 1) {
        return Map<int>(
            ComputeLocalBinsRaw(global_bin - 1, GetNBinsNoOver(axes)),
            [](int bin) { return bin + 1; }
        );
    }

    // Convert our negative index to something positive and 0-based, as that is
    // more convenient to work with. Note, however, that this is _not_
    // equivalent to the virtual_bin that we had before, because what we have
    // here is a count of overflow bins, not of all bins...
    ASSERT(global_bin != 0, "Received an invalid global bin index as input");
    const int corrected_virtual_overflow_bin = -global_bin - 1;

    // ...so we need to retrieve and bring back the regular bin count, and this
    // is where the fun begins.
    //
    // The main difficulty is that the number of regular bins is not fixed as
    // one slides along a histogram axis. Using a 2D binning case as a simple
    // motivating example...
    //
    //    -1   -2   -3   -4   <- No regular bins on the underflow line of axis 1
    //    -5    1    2   -6   <- Some of them on middle lines of axis 1
    //    -7    3    4   -8
    //    -9   -10  -11  -12  <- No regular bins on the overflow line of axis 1
    //
    // As we go to higher dimensions, the geometry becomes more complex, but
    // if we replace "line" with "plane", we get a similar picture in 3D when we
    // slide along axis 2:
    //
    //  No regular bins on the    Some of them on the     No regular bins again
    //    UF plane of axis 2    regular planes of ax.2   on the OF plane of ax.2
    //
    //    -1   -2   -3   -4       -17  -18  -19  -20      -29  -30  -31  -32
    //    -5   -6   -7   -8       -21   1    2   -22      -33  -34  -35  -36
    //    -9   -10  -11  -12      -23   3    4   -24      -37  -37  -39  -40
    //    -13  -14  -15  -16      -25  -26  -27  -28      -41  -42  -43  -44
    //
    // We can generalize this to N dimensions by saying that as we slide along
    // the last axis of an N-d histogram, we see an hyperplane full of overflow
    // bins, then some hyperplanes with regular bins in the "middle" surrounded
    // by overflow bins, then a last hyperplane full of overflow bins.
    //
    // From this, we can devise a recursive algorithm to recover the number of
    // regular bins before the overflow bin we're currently looking at:
    //
    // - Start by processing the last histogram axis.
    // - Ignore the first and last hyperplane on this axis, which only contain
    //   underflow and overflow bins respectively.
    // - Count how many complete hyperplanes of regular bins lie before us on
    //   this axis, which we can do indirectly in our overflow bin based
    //   reasoning by computing the perimeter of the regular region and dividing
    //   our "regular" overflow bin count by that amount.
    // - Now we counted previous hyperplanes on this last histogram axis, but
    //   we need to process the hyperplane that our bin is located in, if any.
    //      * For this, we reduce our overflow bin count to a count of
    //        _unaccounted_ overflow bins in the current hyperplane...
    //      * ...which allows us to recursively continue the computation by
    //        processing the next (well, previous) histogram axis in the context
    //        of this hyperplane, in the same manner as above.
    //
    // Alright, now that the general plan is sorted out, let's compute some
    // quantities that we are going to need, namely the total number of bins per
    // hyperplane (overflow and regular) and the number of regular bins per
    // hyperplane on the hyperplanes that have them.
    //
    std::array<int, NDIMS-1> bins_per_hyperplane;
    std::array<int, NDIMS-1> regular_bins_per_hyperplane;
    int curr_bins_per_hyperplane = 1;
    int curr_regular_bins_per_hyperplane = 1;
    for (std::size_t dim = 0; dim < NDIMS-1; ++dim) {
        curr_regular_bins_per_hyperplane *= axes[dim]->GetNBinsNoOver();
        regular_bins_per_hyperplane[dim] = curr_regular_bins_per_hyperplane;
        curr_bins_per_hyperplane *= axes[dim]->GetNBins();
        bins_per_hyperplane[dim] = curr_bins_per_hyperplane;
    }

    // Given that, and starting from the last axis...
    int unprocessed_previous_overflow_bins = corrected_virtual_overflow_bin;
    int num_regular_bins_before = 0;
    for (std::size_t dim = NDIMS-1; dim > 0; --dim) {
        // Let's start by computing the contribution of the underflow
        // hyperplane (if any), in which we know there will be no regular bins
        const int num_underflow_hyperplanes =
            static_cast<int>(!axes[dim]->CanGrow());
        const int bins_in_underflow_hyperplane =
            num_underflow_hyperplanes * bins_per_hyperplane[dim-1];

        // Next, from the total number of bins per hyperplane and the number of
        // regular bins per hyperplane that has them, we deduce the number of
        // overflow bins per hyperplane that has regular bins.
        const int overflow_bins_per_regular_hyperplane =
            bins_per_hyperplane[dim-1] - regular_bins_per_hyperplane[dim-1];

        // This allows us to answer a key question: are there any under/overflow
        // bins on the hyperplanes that have regular bins? It may not be the
        // case if some axes are growable, and thus don't have overflow bins.
        if (overflow_bins_per_regular_hyperplane != 0) {
            // If so, we start by cutting off the contribution of the underflow
            // and overflow hyperplanes, to focus specifically on regular bins.
            const int overflow_bins_in_regular_hyperplanes =
                std::clamp(
                    unprocessed_previous_overflow_bins
                        - bins_in_underflow_hyperplane,
                    0,
                    overflow_bins_per_regular_hyperplane
                        * axes[dim]->GetNBinsNoOver()
                );

            // We count how many _complete_ "regular" hyperplanes that leaves
            // before us, and account for those in our regular bin count.
            const int num_regular_hyperplanes_before =
                overflow_bins_in_regular_hyperplanes
                    / overflow_bins_per_regular_hyperplane;
            num_regular_bins_before +=
                num_regular_hyperplanes_before
                    * regular_bins_per_hyperplane[dim-1];

            // This only leaves the _current_ hyperplane as a possible source of
            // more regular bins that we haven't accounted for yet. We'll take
            // those into account while processing previous dimensions.
            unprocessed_previous_overflow_bins =
                overflow_bins_in_regular_hyperplanes
                    % overflow_bins_per_regular_hyperplane;
        } else {
            // If there are no overflow bins in regular hyperplane, then the
            // rule changes: observing _one_ overflow bin after the underflow
            // hyperplane means that _all_ regular hyperplanes on this axis are
            // already before us.
            if (unprocessed_previous_overflow_bins >= bins_in_underflow_hyperplane) {
                num_regular_bins_before +=
                    axes[dim]->GetNBinsNoOver()
                        * regular_bins_per_hyperplane[dim-1];
            }

            // In this case, we're done, because the current bin may only lie
            // in the underflow or underflow hyperplane of this axis. Which
            // means that there are no further regular bins to be accounted for
            // in the current hyperplane.
            unprocessed_previous_overflow_bins = 0;
        }

        // No need to continue this loop if we've taken into account all
        // overflow bins that were associated with regular bins.
        if (unprocessed_previous_overflow_bins == 0) break;
    }

    // By the time we reach the first axis, there should only be at most one
    // full row of regular bins before us:
    //
    //    -1  1  2  3  -2
    //     ^            ^
    //     |            |
    //     |        Option 2: one overflow bin before us
    //     |
    // Option 1: no overflow bin before us
    //
    ASSERT(unprocessed_previous_overflow_bins <= 1,
           "Logic error detected in the global -> local bin algorithm");
    num_regular_bins_before +=
        unprocessed_previous_overflow_bins * axes[0]->GetNBinsNoOver();

    // Now that we know the number of regular bins before us, we can add this to
    // to the zero-based overflow bin index that we started with to get a global
    // zero-based bin index accounting for both under/overflow bins and regular
    // bins, just like what we had in the ComputeGlobalBin() implementation.
    const int global_virtual_bin =
        corrected_virtual_overflow_bin + num_regular_bins_before;

    // We can then easily go back to zero-based "virtual" bin indices...
    const std::array<int, NDIMS> virtual_bins =
        ComputeLocalBinsRaw(global_virtual_bin, GetNBins(axes));

    // ...and from that go back to the -1/-2 overflow bin indexing convention.
    return ZipAxisMap<int>(virtual_bins, axes, VirtualBinToLocalBin);
}


// ---


int main() {
    // Set up some test axes
    const EqAxisMock eq_axis_0(3.0, 6.0, 3, true);
    const EqAxisMock eq_axis_1(7.0, 9.5, 5, false);
    const EqAxisMock eq_axis_2(-2., -1., 2, true);

    // Should I print out GetBinFrom and GetBinTo values?
    const bool verbose = true;

    // For the purpose of testing the generality of the code, only access these
    // axes through a minimal virtual interface
    const IAxis& axis_0 = eq_axis_0;
    const IAxis& axis_1 = eq_axis_1;
    const IAxis& axis_2 = eq_axis_2;

    // Enumerate all the bins
    const auto enumerate_bins = [](const IAxis& axis) -> std::vector<int> {
        std::vector<int> all_bins;
        if (!axis.CanGrow()) all_bins.push_back(-1);
        for (int bin = 1; bin <= axis.GetNBinsNoOver(); ++bin) {
            all_bins.push_back(bin);
        }
        if (!axis.CanGrow()) all_bins.push_back(-2);
        return all_bins;
    };
    const std::vector<int> all_bins_0 = enumerate_bins(axis_0);
    const std::vector<int> all_bins_1 = enumerate_bins(axis_1);
    const std::vector<int> all_bins_2 = enumerate_bins(axis_2);

    // We'll test for the 1D, 2D and 3D cases, with these common parts:
    int last_regular_bin;
    int last_overflow_bin;
    auto test_setup = [&] {
        last_regular_bin = 0;
        last_overflow_bin = 0;
    };
    const auto print_array = [](const auto& array) {
        std::cout << '(';
        for (std::size_t dim = 0; dim < array.size() - 1; ++dim) {
            std::cout << array[dim] << ", ";
        }
        std::cout << array.back() << ')';
    };
    auto test_iteration = [&](const auto& bins, const auto& axes) {
        // We start from a set of local bin indices
        std::cout << "On true bin ";
        print_array(bins);
        std::cout << ": ";

        // We are able to convert them into a global histogram bin index, and
        // to convert such an index back into local bin indices.
        int global_bin = ComputeGlobalBin(bins, axes);
        auto computed_bins = ComputeLocalBins(global_bin, axes);
        std::cout << "Global bin " << global_bin << " aka ";
        print_array(computed_bins);

        // Once you have those primitives, you can do everything including
        // global FindBin (= local FindBin + local-to-global idx conversion) or
        // global bin boundary computations (= global-to-local idx
        // conversion + local bin boundary computations)
        if (verbose) {
            std::cout << " goes from ";
            print_array(GetBinFrom(bins, axes));
            std::cout << " to ";
            print_array(GetBinTo(bins, axes));
        }
        std::cout << std::endl;

        // Check correctness of the results
        if (std::all_of(bins.cbegin(), bins.cend(),
                        [](int bin) { return bin >= 1; }))
        {
            ASSERT(global_bin == last_regular_bin + 1,
                   "Computed global regular bin index is incorrect");
            last_regular_bin = global_bin;
        } else {
            ASSERT(global_bin == last_overflow_bin - 1,
                   "Computed global overflow bin index is incorrect");
            last_overflow_bin = global_bin;
        }
        ASSERT(computed_bins == bins,
               "Computed local bin indices are incorrect");
    };

    // Test the 1D case
    std::cout << "1D case:" << std::endl;
    test_setup();
    for (const int bin_0: all_bins_0) {
        test_iteration(std::array{bin_0}, std::array{&axis_0});
    }

    // Test the 2D case
    std::cout << std::endl << "2D case:" << std::endl << "---" << std::endl;
    test_setup();
    for (const int bin_1: all_bins_1) {
        for (const int bin_0: all_bins_0) {
            test_iteration(std::array{bin_0,   bin_1},
                           std::array{&axis_0, &axis_1});
        }
        std::cout << "---" << std::endl;
    }

    // Test the 3D case
    std::cout << std::endl << "3D case:" << std::endl << "===" << std::endl;
    test_setup();
    for (const int bin_2: all_bins_2) {
        std::cout << "---" << std::endl;
        for (const int bin_1: all_bins_1) {
            for (const int bin_0: all_bins_0) {
                test_iteration(std::array{bin_0,   bin_1,   bin_2},
                               std::array{&axis_0, &axis_1, &axis_2});
            }
            std::cout << "---" << std::endl;
        }
        std::cout << "===" << std::endl;
    }
}