Polytope.hpp

Go to the documentation of this file.

#ifndef AI_TOOLBOX_UTILS_POLYTOPE_HEADER_FILE
#define AI_TOOLBOX_UTILS_POLYTOPE_HEADER_FILE
 
#include <array>
#include <optional>
 
#include <Eigen/Dense>
 
#include <AIToolbox/TypeTraits.hpp>
#include <AIToolbox/Utils/Core.hpp>
#include <AIToolbox/Utils/Combinatorics.hpp>
#include <AIToolbox/Utils/IndexMap.hpp>
 
#include <AIToolbox/Utils/LP.hpp>
 
namespace AIToolbox {
    using Hyperplane = Vector;
 
    using Point = ProbabilityVector;
 
    using PointSurface = std::pair<std::vector<Point>, std::vector<double>>;
 
    using CompactHyperplanes = Matrix2D;
 
    inline bool dominates(const Hyperplane & lhs, const Hyperplane & rhs) {
        return (lhs.array() - rhs.array() >= -equalToleranceSmall).minCoeff() ||
               (lhs.array() - rhs.array() >= -lhs.array().min(rhs.array()) * equalToleranceGeneral).minCoeff();
    };
 
    template <typename Iterator, typename P = std::identity>
    Iterator findBestAtPoint(const Point & point, Iterator begin, Iterator end, double * value = nullptr, P p = P{}) {
        auto bestMatch = begin;
        double bestValue = point.dot(std::invoke(p, *bestMatch));
 
        while ( (++begin) < end ) {
            const double currValue = point.dot(std::invoke(p, *begin));
            if ( currValue > bestValue || ( currValue == bestValue && veccmp(std::invoke(p, *begin), std::invoke(p, *bestMatch)) > 0 ) ) {
                bestMatch = begin;
                bestValue = currValue;
            }
        }
        if ( value ) *value = bestValue;
        return bestMatch;
    }
 
    template <typename Iterator, typename P = std::identity>
    Iterator findBestAtSimplexCorner(const size_t corner, Iterator begin, Iterator end, double * value = nullptr, P p = P{}) {
        auto bestMatch = begin;
        double bestValue = std::invoke(p, *bestMatch)[corner];
 
        while ( (++begin) < end ) {
            const double currValue = std::invoke(p, *begin)[corner];
            if ( currValue > bestValue || ( currValue == bestValue && veccmp(std::invoke(p, *begin), std::invoke(p, *bestMatch)) > 0 ) ) {
                bestMatch = begin;
                bestValue = currValue;
            }
        }
        if ( value ) *value = bestValue;
        return bestMatch;
    }
 
    template <typename Iterator, typename P = std::identity>
    Iterator findBestDeltaDominated(const Point & point, const Hyperplane & plane, double delta, Iterator begin, Iterator end, P p = P{}) {
        auto retval = end;
 
        const Hyperplane * maxPlane = &plane;
        double maxVal = point.dot(*maxPlane);
 
        for (auto it = begin; it < end; ++it) {
            const Hyperplane * newPlane = &std::invoke(p, *it);
            const double newVal = point.dot(*newPlane);
            if (newVal > maxVal) {
                const double deltaValue = (newVal - maxVal) / (*newPlane - *maxPlane).norm();
                if (deltaValue > delta) {
                    maxVal = newVal;
                    maxPlane = newPlane;
                    retval = it;
                }
            }
        }
        return retval;
    }
 
    template <typename Iterator, typename P = std::identity>
    Iterator extractBestAtPoint(const Point & point, Iterator begin, Iterator bound, Iterator end, P p = P{}) {
        auto bestMatch = findBestAtPoint(point, begin, end, nullptr, p);
 
        if ( bestMatch >= bound )
            std::iter_swap(bestMatch, bound++);
 
        return bound;
    }
 
    template <typename Iterator, typename P = std::identity>
    Iterator extractBestAtSimplexCorners(const size_t S, Iterator begin, Iterator bound, Iterator end, P p = P{}) {
        if ( end == bound ) return bound;
 
        // For each corner
        for ( size_t s = 0; s < S; ++s ) {
            auto bestMatch = findBestAtSimplexCorner(s, begin, end, nullptr, p);
 
            if ( bestMatch >= bound )
                std::iter_swap(bestMatch, bound++);
        }
        return bound;
    }
 
    template <typename PIterator, typename VIterator, typename P = std::identity>
    PIterator extractBestUsefulPoints(PIterator pbegin, PIterator pend, VIterator begin, VIterator end, P p = P{}) {
        const auto pointsN  = std::distance(pbegin, pend);
        const auto entriesN = std::distance(begin, end);
 
        std::vector<std::pair<PIterator, double>> bestValues(entriesN, {pend, std::numeric_limits<double>::lowest()});
        const auto maxBound = pointsN < entriesN ? pend : pbegin + entriesN;
 
        // So the idea here is that we advance IT only if we found a Point
        // which supports a previously unsupported Hyperplane. This allows us
        // to avoid doing later work for compacting the points before the
        // bound.
        //
        // If instead the found Point takes into consideration an already
        // supported Hyperplane, then it either is better or not. If it's
        // better, we swap it with whatever was before. In both cases, the
        // Point to discard ends up at the end and we decrease the bound.
        auto it = pbegin;
        auto bound = pend;
        while (it < bound && it < maxBound) {
            double value;
            const auto vId = std::distance(begin, findBestAtPoint(*it, begin, end, &value, p));
            if (bestValues[vId].second < value) {
                if (bestValues[vId].first == pend) {
                    bestValues[vId] = {it++, value};
                    continue;
                } else {
                    bestValues[vId].second = value;
                    std::iter_swap(bestValues[vId].first, it);
                }
            }
            std::iter_swap(it, --bound);
        }
        if (it == bound) return it;
 
        // If all Hyperplanes have been supported by at least one Point, then
        // we can finish up the rest with less swaps and checks. Here we only
        // swap with the best if needed, otherwise we don't have to do
        // anything.
        //
        // This is because we can return one Point per Hyperplane at the most,
        // so if we're here the bound is not going to move anyway.
        while (it < bound) {
            double value;
            const auto vId = std::distance(begin, findBestAtPoint(*it, begin, end, &value, p));
            if (bestValues[vId].second < value) {
                bestValues[vId].second = value;
                std::iter_swap(bestValues[vId].first, it);
            }
            ++it;
        }
        return maxBound;
    }
 
    template <typename NewIt, typename OldIt, typename P1 = std::identity, typename P2 = std::identity>
    PointSurface findVerticesNaive(NewIt beginNew, NewIt endNew, OldIt alphasBegin, OldIt alphasEnd, P1 p1 = P1{}, P2 p2 = P2{}) {
        PointSurface vertices;
 
        const size_t alphasSize = std::distance(alphasBegin, alphasEnd);
        if (alphasSize == 0) return vertices;
        const size_t S = std::invoke(p2, *alphasBegin).size();
 
        // This enumerator allows us to compute all possible subsets of S-1
        // elements. We use it on both the alphas, and the boundaries, thus the
        // number of elements we iterate over is alphasSize + S.
        SubsetEnumerator enumerator(S - 1, 0ul, alphasSize + S);
 
        // This is the matrix on the left side of Ax = b (where A is m)
        Matrix2D m(S + 1, S + 1);
        m.row(0)[S] = -1; // First row is always a vector
 
        Vector boundary(S+1);
        boundary[S] = 0.0; // The boundary doesn't care about the value
 
        // This is the vector on the right side of Ax = b
        Vector b(S+1); b.setZero();
 
        Vector result(S+1);
 
        // Common matrix/vector setups
 
        for (auto newVIt = beginNew; newVIt != endNew; ++newVIt) {
            m.row(0).head(S) = std::invoke(p1, *newVIt);
 
            enumerator.reset();
 
            // Get subset of planes, find corner with LU
            size_t last = 0;
            while (enumerator.isValid()) {
                // Reset boundaries to care about all dimensions
                boundary.head(S).fill(1.0);
                size_t counter = last + 1;
                // Note that we start from last to avoid re-copying vectors
                // that are already in the matrix in their correct place.
                for (auto i = last; i < enumerator->size(); ++i) {
                    // For each value in the enumerator, if it is less than
                    // alphasSize it is referring to an alphaVector we need to
                    // take into account.
                    const auto index = (*enumerator)[i];
                    if (index < alphasSize) {
                        // Copy the right vector in the matrix.
                        m.row(counter).head(S) = std::invoke(p2, *std::next(alphasBegin, index));
                        m.row(counter)[S] = -1;
                        ++counter;
                    } else {
                        // We limit the index-th dimension (minus alphasSize to scale in a 0-S range)
                        boundary[index - alphasSize] = 0.0;
                    }
                }
                m.row(counter) = boundary;
                b[counter] = 1.0;
                ++counter;
 
                // Note that we only need to consider the first "counter" rows,
                // as the boundaries get merged in a single one.
                result = m.topRows(counter).colPivHouseholderQr().solve(b.head(counter));
 
                b[counter-1] = 0.0;
 
                // Add to found only if valid, otherwise skip.
                const double max = result.head(S).maxCoeff();
                if ((result.head(S).array() >= 0).all() && (max < 1.0) && checkDifferentSmall(max, 1.0)) {
                    vertices.first.emplace_back(result.head(S));
                    vertices.second.emplace_back(result[S]);
                }
 
                // Advance, and take the id of the first index changed in the
                // next combination.
                last = enumerator.advance();
 
                // If the index went over the alpha list, then we'd only have
                // boundaries, but we don't care about those cases (since we
                // assume we already have the corners of the simplex computed).
                // Thus, terminate.
                if ((*enumerator)[0] >= alphasSize) break;
            }
        }
        return vertices;
    }
 
    template <typename Range, typename P = std::identity>
    PointSurface findVerticesNaive(const Range & range, P p = P{}) {
        PointSurface retval;
 
        std::array<size_t, 1> indexToSkip;
 
        auto goodBegin = range.cbegin();
        for (size_t i = 0; i < range.size(); ++i, ++goodBegin) {
            // For each alpha, we find its vertices against the others.
            indexToSkip[0] = i;
            IndexSkipMap map(&indexToSkip, range);
 
            const auto cbegin = map.cbegin();
            const auto cend   = map.cend();
 
            // Note that the range we pass here is made by a single vector.
            auto newVertices = findVerticesNaive(goodBegin, goodBegin + 1, cbegin, cend, p, p);
 
            retval.first.insert(std::end(retval.first),
                std::make_move_iterator(std::begin(newVertices.first)),
                std::make_move_iterator(std::end(newVertices.first))
            );
            retval.second.insert(std::end(retval.second),
                std::make_move_iterator(std::begin(newVertices.second)),
                std::make_move_iterator(std::end(newVertices.second))
            );
        }
        return retval;
    }
 
    double computeOptimisticValue(const Point & p, const std::vector<Point> & points, const std::vector<double> & values);
 
    std::tuple<double, Vector> LPInterpolation(const Point & p, const CompactHyperplanes & ubQ, const PointSurface & ubV);
 
    std::tuple<double, Vector> sawtoothInterpolation(const Point & p, const CompactHyperplanes & ubQ, const PointSurface & ubV);
 
    class WitnessLP {
        public:
            WitnessLP(size_t S);
 
            void addOptimalRow(const Hyperplane & v);
 
            std::optional<Point> findWitness(const Hyperplane & v);
 
            void reset();
 
            void allocate(size_t rows);
 
        private:
            size_t S;
            LP lp_;
    };
}
 
#endif