blitz.hh

#ifndef VTESTBED_BASE_BLITZ_HH
#define VTESTBED_BASE_BLITZ_HH

#include <array>
#include <cmath>
#include <functional>
#include <limits>
#include <ostream>

#include <blitz/array.h>

#include <vtestbed/base/bstream.hh>
#include <vtestbed/base/constants.hh>

namespace blitz {

    inline bool
    isfinite(float rhs) noexcept {
        return ::std::isfinite(rhs);
    }

    inline bool
    isfinite(double rhs) noexcept {
        return ::std::isfinite(rhs);
    }

    BZ_DECLARE_FUNCTION(isfinite);

    inline int
    is_power_of_two(int n) {
        return n > 0 && (n & (n-1)) == 0;
    }

    BZ_DECLARE_FUNCTION(is_power_of_two);

    inline int
    next_power_of_two(int n) {
        if (n > std::numeric_limits<int>::max()/2) {
            throw std::overflow_error(__func__);
        }
        int m = 1;
        while (m < n) {
            m *= 2;
        }
        return m;
    }

    BZ_DECLARE_FUNCTION(next_power_of_two);

    inline float radians_to_degrees(float x) {
        return x*360.f/vtb::base::constants<float>::pi(2);
    }

    inline double radians_to_degrees(double x) {
        return x*360.0/vtb::base::constants<double>::pi(2);
    }

    BZ_DECLARE_FUNCTION(radians_to_degrees);

    template <class T, int N>
    inline vtb::core::bstream&
    operator<<(vtb::core::bstream& out, const blitz::TinyVector<T,N>& rhs) {
        for (int i=0; i<N; ++i) {
            out << rhs(i);
        }
        return out;
    }

    template <class T, int N>
    inline vtb::core::bstream&
    operator>>(vtb::core::bstream& in, blitz::TinyVector<T,N>& rhs) {
        for (int i=0; i<N; ++i) {
            in >> rhs(i);
        }
        return in;
    }

    template <class T, int N, int M>
    inline vtb::core::bstream&
    operator<<(vtb::core::bstream& out, const blitz::TinyMatrix<T,N,M>& rhs) {
        for (int i=0; i<N; ++i) { for (int j=0; j<M; ++j) { out << rhs(i,j); } }
        return out;
    }

    template <class T, int N, int M>
    inline vtb::core::bstream&
    operator>>(vtb::core::bstream& in, blitz::TinyMatrix<T,N,M>& rhs) {
        for (int i=0; i<N; ++i) { for (int j=0; j<M; ++j) { in >> rhs(i,j); } }
        return in;
    }

    template <class T, int N>
    inline vtb::core::bstream&
    operator<<(vtb::core::bstream& out, const blitz::Array<T,N>& rhs) {
        blitz::TinyVector<uint32_t,N> shape = rhs.shape();
        out << shape;
        const T* data = rhs.data();
        const int n = rhs.numElements();
        for (int i=0; i<n; ++i) {
            out << data[i];
        }
        return out;
    }

    template <class T, int N>
    inline vtb::core::bstream&
    operator>>(vtb::core::bstream& in, blitz::Array<T,N>& rhs) {
        blitz::TinyVector<uint32_t,N> shape;
        shape = 0;
        in >> shape;
        rhs.resize(shape);
        T* data = rhs.data();
        const int n = rhs.numElements();
        for (int i=0; i<n; ++i) {
            in >> data[i];
        }
        return in;
    }

    template <class T, int N>
    inline blitz::TinyVector<T,N>
    zero(blitz::TinyVector<T,N>) {
        return blitz::TinyVector<T,N>{T{0}};
    }

    template <class T, int N>
    inline blitz::TinyVector<T,N>
    one(blitz::TinyVector<T,N>) {
        return blitz::TinyVector<T,N>{T{1}};
    }

    template <class T, int N, int M>
    inline TinyVector<T,N>
    product(const TinyMatrix<T,N,M>& lhs, const TinyVector<T,M>& rhs) {
        TinyVector<T,N> result;
        for (int i=0; i<N; ++i) {
            T sum{0};
            for (int j=0; j<M; ++j) {
                sum += lhs(i,j)*rhs(j);
            }
            result(i) = sum;
        }
        return result;
    }

    template <class T, int N, int M>
    inline blitz::TinyVector<T,N>
    operator*(const TinyMatrix<T,N,M>& lhs, const TinyVector<T,M>& rhs) {
        return product(lhs, rhs);
    }

    template <class T, int N, int M, int K>
    inline TinyMatrix<T,N,M>
    product(const TinyMatrix<T,N,K>& lhs, const TinyMatrix<T,K,M>& rhs) {
        TinyMatrix<T,N,M> result;
        for (int i=0; i<N; ++i) {
            for (int j=0; j<M; ++j) {
                T sum{};
                for (int k=0; k<K; ++k) { sum += lhs(i,k)*rhs(k,j); }
                result(i,j) = sum;
            }
        }
        return result;
    }

    template <class T, int N, int M, int K>
    inline TinyMatrix<T,N,M>
    operator*(const TinyMatrix<T,N,K>& lhs, const TinyMatrix<T,K,M>& rhs) {
        return product(lhs, rhs);
    }

    template <class T, int N>
    inline bool
    near(
        const blitz::TinyVector<T,N>& lhs,
        const blitz::TinyVector<T,N>& rhs,
        const T eps
    ) {
        using blitz::all;
        using blitz::abs;
        return !any(abs(lhs - rhs) > eps);
    }

    template <class T, int N>
    inline bool
    near(
        const blitz::TinyVector<T,N>& lhs,
        const T rhs,
        const T eps
    ) {
        using blitz::all;
        using blitz::abs;
        return !any(abs(lhs - rhs) > eps);
    }

    template <int N>
    inline bool
    empty(const RectDomain<N>& rect){
        return any(rect.lbound() > rect.ubound());
    }

    template <class T,int N,int M>
    MinMaxValue<TinyVector<T,N>>
    minmax(Array<TinyVector<T,N>,M> x) {
        MinMaxValue<TinyVector<T,N>> result;
        result.min = std::numeric_limits<T>::max();
        result.max = std::numeric_limits<T>::lowest();
        if (x.size() == 0) {
            return result;
        }
        for (const auto& v : x) {
            for (int i=0; i<N; ++i) {
                T m = v(i);
                if (m < result.min(i)) {
                    result.min(i) = m;
                }
                if (result.max(i) < m) {
                    result.max(i) = m;
                }
            }
        }
        return result;
    }

    template <class T, int N>
    inline blitz::TinyVector<T,N>
    min_vector(Array<blitz::TinyVector<T,N>,1> x) {
        blitz::TinyVector<T,N> m{std::numeric_limits<T>::max()};
        for (const auto& v : x) {
            m = blitz::min(m, v);
        }
        return m;
    }

    template <class T, int N>
    inline blitz::TinyVector<T,N>
    max_vector(Array<blitz::TinyVector<T,N>,1> x) {
        blitz::TinyVector<T,N> m{std::numeric_limits<T>::min()};
        for (const auto& v : x) {
            m = blitz::max(m, v);
        }
        return m;
    }

    template <class T, int N, int M>
    inline TinyVector<T,N*M>
    flatten(const TinyVector<T,N> x[M]) {
        TinyVector<T,N*M> result;
        for (int i=0; i<M; ++i) {
            for (int j=0; j<N; ++j) {
                result(i*N + j) = x[i][j];
            }
        }
        return result;
    }

    template <class T, int N, class ... Args>
    inline TinyVector<T,N*(sizeof...(Args)+1)>
    flatten(const TinyVector<T,N>& head, const Args& ... tail) {
        const TinyVector<T,N> args[sizeof...(Args)+1] = {head, tail...};
        return flatten<T,N,sizeof...(Args)+1>(args);
    }

    namespace bits {

        template <class T, int N, int M>
        struct Tie {
            std::reference_wrapper<TinyVector<T,N>> arr[M];
            inline void
            operator=(const TinyVector<T,N*M>& rhs) {
                for (int i=0; i<M; ++i) {
                    for (int j=0; j<N; ++j) {
                        arr[i].get()[j] = rhs(i*N + j);
                    }
                }
            }
        };

    }

    template <class T, int N, class ... Args>
    inline bits::Tie<T,N,sizeof...(Args)+1>
    tie(TinyVector<T,N>& head, Args& ... tail) {
        return {{head, tail...}};
    }

    template <class T, int N>
    inline blitz::TinyMatrix<T,N,N>
    identity_matrix() {
        blitz::TinyMatrix<T,N,N> result;
        result = 0;
        for (int i=0; i<N; ++i) {
            result(i,i) = 1;
        }
        return result;
    }

    template <class T, int N>
    inline blitz::TinyMatrix<T,N,N>
    const_matrix(const std::array<T,N*N>& rhs) {
        blitz::TinyMatrix<T,N,N> m;
        for (int i=0; i<N*N; ++i) { m.data()[i] = rhs[i]; }
        return m;
    }

    template <class T, int N>
    inline blitz::TinyVector<T,N>
    const_vector(T value) {
        blitz::TinyVector<T,N> result;
        result = value;
        return result;
    }

    template <class T, int N>
    inline blitz::TinyVector<T,1>
    slice(const blitz::TinyVector<T,N>& x, int i) {
        return {x(i)};
    }

    template <class T, int N>
    inline blitz::TinyVector<T,2>
    slice(const blitz::TinyVector<T,N>& x, int i, int j) {
        return {x(i), x(j)};
    }

    template <class T, int N>
    inline blitz::TinyVector<T,3>
    slice(const blitz::TinyVector<T,N>& x, int i, int j, int k) {
        return {x(i), x(j), x(k)};
    }

    template <class T, int N>
    inline blitz::TinyVector<T,4>
    slice(const blitz::TinyVector<T,N>& x, int i, int j, int k, int l) {
        return {x(i), x(j), x(k), x(l)};
    }

    template <class T>
    inline T
    clamp(T x, T x_min, T x_max) {
        if (x < x_min) x = x_min;
        if (x > x_max) x = x_max;
        return x;
    }

    template <class T, int N>
    inline blitz::TinyVector<T,N>
    clamp(
        const blitz::TinyVector<T,N>& x,
        const blitz::TinyVector<T,N>& x_min,
        const blitz::TinyVector<T,N>& x_max
    ) {
        return blitz::max(x_min, blitz::min(x_max, x));
    }

    inline int
    div_ceil(int lhs, int rhs) noexcept {
        return lhs/rhs + (lhs%rhs == 0 ? 0 : 1);
    }

    BZ_DECLARE_FUNCTION2(div_ceil);

    template <class T>
    class Statistics {

    private:
        T _sum1{0};
        T _sum2{0};
        int _count = 0;

    public:

        inline int
        count() const noexcept {
            return this->_count;
        }

        inline const T&
        mean() const noexcept {
            return this->_sum1;
        }

        inline T
        variance(T addon) const noexcept {
            return this->_sum2 / (this->_count + addon);
        }

        inline T
        population_variance() const noexcept {
            return this->variance(0);
        }

        inline T
        unbiased_variance() const noexcept {
            return this->variance(-1);
        }

        inline T
        variance() const noexcept {
            return this->unbiased_variance();
        }

        inline void
        update(T x) {
            ++_count;
            const T delta = x - _sum1;
            _sum1 += delta / _count;
            const T delta2 = x - _sum1;
            _sum2 += delta*delta2;
        }

        inline void
        clear() {
            _sum1 = 0;
            _sum2 = 0;
            _count = 0;
        }

        inline friend std::ostream&
        operator<<(std::ostream& out, const Statistics& rhs) {
            return out
                << "mean=" << rhs.mean()
                << ",variance=" << rhs.variance();
        }

    };

    template <class T, int N>
    inline Statistics<T>
    statistics(const Array<T,N>& rhs) {
        const int n = rhs.numElements();
        const T* data = rhs.data();
        Statistics<T> stats;
        for (int i=0; i<n; ++i) {
            stats.update(data[i]);
        }
        return stats;
    }

    template <class T, int N>
    inline T
    variance(const Array<T,N>& rhs) {
        return statistics<T,N>(rhs).variance();
    }

    template <class T, int N>
    inline const Array<T,N>&
    real(const Array<T,N>& x) {
        return x;
    }

    template <class T, int N>
    inline T
    length(const blitz::TinyVector<T,N>& x) {
        using blitz::abs;
        using blitz::max;
        using blitz::pow2;
        using blitz::sum;
        using std::sqrt;
        const T s = max(abs(x));
        if (s == T{}) {
            return s;
        }
        return s*sqrt(sum(pow2(x/s)));
    }

    template <class T>
    inline T
    length(const blitz::TinyVector<T,2>& x) {
        return blitz::hypot(x(0), x(1));
    }

    template <class T>
    inline T
    length(const blitz::TinyVector<T,1>& x) {
        return std::abs(x(0));
    }

    template <class T, int N>
    inline T
    sum_of_squares(const blitz::TinyVector<T,N>& x) {
        return sum(pow2(x));
    }

    template <class T>
    inline T
    sum_of_squares(const blitz::TinyVector<T,2>& x) {
        return pow2(x(0)) + pow2(x(1));
    }

    template <class T>
    inline T
    sum_of_squares(const blitz::TinyVector<T,3>& x) {
        return pow2(x(0)) + pow2(x(1)) + pow2(x(2));
    }

    template <class T>
    inline T
    sum_of_squares(const blitz::TinyVector<T,4>& x) {
        return pow2(x(0)) + pow2(x(1)) + pow2(x(2)) + pow2(x(3));
    }

    template <class T, int N>
    inline T
    distance(const blitz::TinyVector<T,N>& v1, const blitz::TinyVector<T,N>& v2) {
        return length(blitz::TinyVector<T,N>(v2 - v1));
    }

    template <class T, int N>
    inline void
    transpose(TinyMatrix<T,N,N>& m) {
        using std::swap;
        for (int i=0; i<N; ++i) {
            for (int j=i+1; j<N; ++j) {
                swap(m(i,j), m(j,i));
            }
        }
    }

    template <class T, int N>
    inline TinyVector<T,N-1>
    puncture(const TinyVector<T,N>& v, int j) {
        TinyVector<T,N-1> r;
        for (int i=0; i<j; ++i) { r(i) = v(i); }
        for (int i=j+1; i<N; ++i) { r(i-1) = v(i); }
        return r;
    }

    template <class T, int N>
    inline TinyVector<T,N+1>
    restore(const TinyVector<T,N>& v, int j) {
        TinyVector<T,N+1> r;
        for (int i=0; i<j; ++i) { r(i) = v(i); }
        r(j) = T{};
        for (int i=j; i<N; ++i) { r(i+1) = v(i); }
        return r;
    }

    template <class T>
    inline TinyVector<T,1>
    cross(const TinyVector<T,2>& a, const TinyVector<T,2>& b) {
        return a(0)*b(1) - a(1)*b(0);
    }

}

#endif // vim:filetype=cpp