QPot.h

Go to the documentation of this file.

 
#ifndef QPOT_H
#define QPOT_H
 
#include <algorithm>
#include <stdexcept>
#include <string>
 
#define Q(x) #x
#define QUOTE(x) Q(x)
 
#define QPOT_WARNING_IMPL(message, file, line, function) \
    std::cout << std::string(file) + ":" + std::to_string(line) + " (" + std::string(function) + \
                     ")" + ": " message ") \n\t";
 
#define QPOT_ASSERT_IMPL(expr, file, line, function) \
    if (!(expr)) { \
        throw std::runtime_error( \
            std::string(file) + ":" + std::to_string(line) + " (" + std::string(function) + ")" + \
            ": assertion failed (" #expr ") \n\t"); \
    }
#ifdef QPOT_ENABLE_ASSERT
#define QPOT_ASSERT(expr) QPOT_ASSERT_IMPL(expr, __FILE__, __LINE__, __FUNCTION__)
#else
#define QPOT_ASSERT(expr)
#endif
 
#ifdef QPOT_ENABLE_DEBUG
#define QPOT_DEBUG(expr) QPOT_ASSERT_IMPL(expr, __FILE__, __LINE__, __FUNCTION__)
#else
#define QPOT_DEBUG(expr)
#endif
 
#define QPOT_REQUIRE(expr) QPOT_ASSERT_IMPL(expr, __FILE__, __LINE__, __FUNCTION__)
 
#ifdef QPOT_ENABLE_WARNING_PYTHON
#define QPOT_WARNING_PYTHON(message) QPOT_WARNING_IMPL(message, __FILE__, __LINE__, __FUNCTION__)
#else
#define QPOT_WARNING_PYTHON(message)
#endif
 
#ifndef QPOT_VERSION
#define QPOT_VERSION "@PROJECT_VERSION@"
#endif
 
namespace QPot {
 
namespace detail {
 
inline std::string unquote(const std::string& arg)
{
    std::string ret = arg;
    ret.erase(std::remove(ret.begin(), ret.end(), '\"'), ret.end());
    return ret;
}
 
} // namespace detail
 
inline std::string version()
{
    return detail::unquote(std::string(QUOTE(QPOT_VERSION)));
}
 
inline std::vector<std::string> version_dependencies()
{
    std::vector<std::string> ret;
 
    ret.push_back("qpot=" + version());
 
    ret.push_back(
        "xtensor=" + detail::unquote(std::string(QUOTE(XTENSOR_VERSION_MAJOR))) + "." +
        detail::unquote(std::string(QUOTE(XTENSOR_VERSION_MINOR))) + "." +
        detail::unquote(std::string(QUOTE(XTENSOR_VERSION_PATCH))));
 
#ifdef XSIMD_VERSION_MAJOR
    ret.push_back(
        "xsimd=" + detail::unquote(std::string(QUOTE(XSIMD_VERSION_MAJOR))) + "." +
        detail::unquote(std::string(QUOTE(XSIMD_VERSION_MINOR))) + "." +
        detail::unquote(std::string(QUOTE(XSIMD_VERSION_PATCH))));
#endif
 
#ifdef XTL_VERSION_MAJOR
    ret.push_back(
        "xtl=" + detail::unquote(std::string(QUOTE(XTL_VERSION_MAJOR))) + "." +
        detail::unquote(std::string(QUOTE(XTL_VERSION_MINOR))) + "." +
        detail::unquote(std::string(QUOTE(XTL_VERSION_PATCH))));
#endif
 
#if defined(XTENSOR_PYTHON_VERSION_MAJOR)
    ret.push_back(
        "xtensor-python=" + detail::unquote(std::string(QUOTE(XTENSOR_PYTHON_VERSION_MAJOR))) +
        "." + detail::unquote(std::string(QUOTE(XTENSOR_PYTHON_VERSION_MINOR))) + "." +
        detail::unquote(std::string(QUOTE(XTENSOR_PYTHON_VERSION_PATCH))));
#endif
 
    std::sort(ret.begin(), ret.end(), std::greater<std::string>());
 
    return ret;
}
 
namespace iterator {
 
template <class It, class T, class R = size_t>
inline R
lower_bound(const It first, const It last, const T& value, R guess = 0, size_t proximity = 10)
{
    if (proximity == 0) {
        if (value <= *(first)) {
            return 0;
        }
        return std::lower_bound(first, last, value) - first - 1;
    }
 
    if (*(first + guess) < value && value <= *(first + guess + 1)) {
        return guess;
    }
 
    size_t l = guess > proximity ? guess - proximity : 0;
    size_t r = std::min(guess + proximity, static_cast<size_t>(last - first - 1));
 
    if (*(first + l) < value && *(first + r) >= value) {
        return std::lower_bound(first + l, first + r, value) - first - 1;
    }
    else if (value <= *(first)) {
        return 0;
    }
    else {
        return std::lower_bound(first, last, value) - first - 1;
    }
}
 
} // namespace iterator
 
namespace inplace {
 
template <class T, class V, class R>
inline void lower_bound(const T& matrix, const V& value, R& index, size_t proximity = 10)
{
    QPOT_ASSERT(value.dimension() == matrix.dimension() - 1);
    QPOT_ASSERT(value.dimension() == index.dimension());
 
    auto nd = value.dimension();
    auto stride = matrix.shape(nd);
    auto n = value.size();
 
#ifdef QPOT_ENABLE_ASSERT
    for (decltype(nd) i = 0; i < nd; ++i) {
        QPOT_ASSERT(matrix.shape(i) == value.shape(i));
        QPOT_ASSERT(matrix.shape(i) == index.shape(i));
    }
#endif
 
    for (decltype(n) i = 0; i < n; ++i) {
        index.flat(i) = iterator::lower_bound(
            &matrix.flat(i * stride),
            &matrix.flat(i * stride) + stride,
            value.flat(i),
            index.flat(i),
            proximity);
    }
}
 
template <class V, class I>
inline void cumsum_chunk(V& cumsum, const V& delta, const I& shift)
{
    QPOT_ASSERT(cumsum.dimension() >= 1);
    QPOT_ASSERT(cumsum.dimension() == delta.dimension());
 
    if (delta.size() == 0) {
        return;
    }
 
    size_t dim = cumsum.dimension();
    size_t n = cumsum.shape(dim - 1);
    size_t ndelta = delta.shape(dim - 1);
 
    for (size_t i = 0; i < shift.size(); ++i) {
 
        if (shift.flat(i) == 0) {
            continue;
        }
 
        auto* d = &delta.flat(i * ndelta);
        auto* c = &cumsum.flat(i * n);
 
        if (shift.flat(i) > 0) {
            QPOT_ASSERT(shift.flat(i) <= n);
            QPOT_ASSERT(ndelta >= shift.flat(i));
            size_t nadd = static_cast<size_t>(shift.flat(i));
            size_t nkeep = n - nadd;
            auto offset = *(c + n - 1);
            std::copy(c + n - nkeep, c + n, c);
            std::copy(d, d + nadd, c + nkeep);
            *(c + nkeep) += offset;
            std::partial_sum(c + nkeep, c + n, c + nkeep);
        }
        else {
            QPOT_ASSERT(-shift.flat(i) < n);
            QPOT_ASSERT(ndelta > -shift.flat(i));
            size_t nadd = static_cast<size_t>(-shift.flat(i));
            size_t nkeep = n - nadd;
            auto offset = *(c);
            std::copy(c, c + nkeep, c + nadd);
            std::copy(d, d + nadd + 1, c);
            std::partial_sum(c, c + nadd + 1, c);
            offset -= *(c + nadd);
            std::transform(c, c + nadd + 1, c, [&](auto& v) { return v + offset; });
        }
    }
}
 
} // namespace inplace
 
template <class T, class V, class R>
inline R lower_bound(const T& matrix, const V& value, const R& index, size_t proximity = 10)
{
    R ret = index;
    inplace::lower_bound(matrix, value, ret, proximity);
    return ret;
}
 
template <class T, class V, class R>
inline R lower_bound(const T& matrix, const V& value)
{
    R ret = xt::zeros<typename R::value_type>(value.shape());
    inplace::lower_bound(matrix, value, ret, 0);
    return ret;
}
 
template <class V, class I>
inline V cumsum_chunk(const V& cumsum, const V& delta, const I& shift)
{
    V ret = cumsum;
    inplace::cumsum_chunk(ret, delta, shift);
    return ret;
}
 
} // namespace QPot
 
#endif