globals.hpp

Go to the documentation of this file.

//=======================================================================
// Copyright (c) 2014-2023 Baptiste Wicht
// Distributed under the terms of the MIT License.
// (See accompanying file LICENSE or copy at
//  http://opensource.org/licenses/MIT)
//=======================================================================

#pragma once

#include "etl/temporary.hpp"

#include "etl/impl/decomposition.hpp"
#include "etl/impl/det.hpp"

#include "etl/impl/egblas/shuffle.hpp"

namespace etl {

template <etl_expr E>
bool is_square(E&& expr) {
    return is_2d<E> && etl::dim<0>(expr) == etl::dim<1>(expr);
}

template <etl_expr E>
bool is_real_matrix([[maybe_unused]] E&& expr) {
    return !is_complex<E>;
}

template <etl_expr E>
bool is_complex_matrix([[maybe_unused]] E&& expr) {
    return is_complex<E>;
}

template <etl_expr E>
bool is_rectangular(E&& expr) {
    return is_2d<E> && etl::dim<0>(expr) != etl::dim<1>(expr);
}

template <etl_expr E>
bool is_sub_square(E&& expr) {
    return is_3d<E> && etl::dim<1>(expr) == etl::dim<2>(expr);
}

template <etl_expr E>
bool is_sub_rectangular(E&& expr) {
    return is_3d<E> && etl::dim<1>(expr) != etl::dim<2>(expr);
}

template <typename E>
bool is_uniform(E&& expr) {
    if (!etl::size(expr)) {
        return false;
    }

    auto first = *expr.begin();

    for (auto v : expr) {
        if (v != first) {
            return false;
        }
    }

    return true;
}

template <typename E>
bool is_permutation_matrix(E&& expr) {
    if (!is_square(expr)) {
        return false;
    }

    //Conditions:
    //a) Must be a square matrix
    //b) Every row must have one 1
    //c) Every column must have one 1

    for (size_t i = 0; i < etl::dim<0>(expr); ++i) {
        auto sum = value_t<E>(0);
        for (size_t j = 0; j < etl::dim<0>(expr); ++j) {
            if (expr(i, j) != value_t<E>(0) && expr(i, j) != value_t<E>(1)) {
                return false;
            }

            sum += expr(i, j);
        }

        if (sum != value_t<E>(1)) {
            return false;
        }
    }

    for (size_t j = 0; j < etl::dim<0>(expr); ++j) {
        auto sum = value_t<E>(0);
        for (size_t i = 0; i < etl::dim<0>(expr); ++i) {
            sum += expr(i, j);
        }

        if (sum != value_t<E>(1)) {
            return false;
        }
    }

    return true;
}

template <etl_expr L, etl_expr R>
bool operator==(L&& lhs, R&& rhs) {
    // Both expressions must have the same number of dimensions
    if (etl::dimensions(lhs) != etl::dimensions(rhs)) {
        return false;
    }

    // The dimensions must be the same
    for (size_t i = 0; i < etl::dimensions(rhs); ++i) {
        if (etl::dim(lhs, i) != etl::dim(rhs, i)) {
            return false;
        }
    }

    // At this point, the values are necessary for the comparison
    force(lhs);
    force(rhs);

    // Note: Ideally, we should use std::equal, but this is significantly
    // faster to compile

    for (size_t i = 0; i < etl::size(lhs); ++i) {
        if (lhs[i] != rhs[i]) {
            return false;
        }
    }

    return true;
}

template <etl_expr L, etl_expr R>
bool operator!=(L&& lhs, R&& rhs) {
    return !(lhs == rhs);
}

template <std::floating_point T, std::floating_point TE = T>
inline bool approx_equals_float(T a, T b, TE epsilon) {
    using std::fabs;

    const auto abs_a    = fabs(a);
    const auto abs_b    = fabs(b);
    const auto abs_diff = fabs(a - b);

    // Note: min for floating points is the min normalized value
    static constexpr T min_normal = std::numeric_limits<T>::min();
    static constexpr T max        = std::numeric_limits<T>::max();

    if (a == b) {
        // This should handle infinities properly
        return true;
    } else if (a == 0 || b == 0 || abs_diff < min_normal) {
        // a or b is zero or both are extremely close to it
        // relative error is less meaningful here
        return abs_diff < epsilon;
    } else { // use relative error
        return (abs_diff / std::min(abs_a + abs_b, max)) < epsilon;
    }
}

template <etl_expr L, etl_expr E>
bool approx_equals(L&& lhs, E&& rhs, value_t<L> eps) {
    // Both expressions must have the same number of dimensions
    if constexpr (etl::dimensions<L>() != etl::dimensions<E>()) {
        return false;
    } else {
        // The dimensions must be the same
        for (size_t i = 0; i < etl::dimensions(rhs); ++i) {
            if (etl::dim(lhs, i) != etl::dim(rhs, i)) {
                return false;
            }
        }

        // At this point, the values are necessary for the comparison
        force(lhs);
        force(rhs);

        for (size_t i = 0; i < etl::size(lhs); ++i) {
            if (!approx_equals_float(lhs[i], rhs[i], eps)) {
                return false;
            }
        }

        return true;
    }
}

template <etl_expr E>
value_t<E> trace(E&& expr) {
    assert_square(expr);

    auto value = value_t<E>();

    for (size_t i = 0; i < etl::dim<0>(expr); ++i) {
        value += expr(i, i);
    }

    return value;
}

template <etl_expr E>
value_t<E> determinant(E&& expr) {
    assert_square(expr);

    return detail::det_impl::apply(expr);
}

template <etl_expr AT, etl_expr LT, etl_expr UT, etl_expr PT>
bool lu(const AT& A, LT& L, UT& U, PT& P) {
    // All matrices must be square
    if (!is_square(A) || !is_square(L) || !is_square(U) || !is_square(P)) {
        return false;
    }

    // All matrices must be of the same dimension
    if (etl::dim(A, 0) != etl::dim(L, 0) || etl::dim(A, 0) != etl::dim(U, 0) || etl::dim(A, 0) != etl::dim(P, 0)) {
        return true;
    }

    detail::lu_impl::apply(A, L, U, P);

    return true;
}

template <etl_expr AT, etl_expr QT, etl_expr RT>
bool qr(AT& A, QT& Q, RT& R) {
    // A and R have the same dimensions
    if (etl::dim(A, 0) != etl::dim(R, 0) || etl::dim(A, 1) != etl::dim(R, 1)) {
        return false;
    }

    // A and Q have corresponding first dimensions and Q is square
    if (etl::dim(A, 0) != etl::dim(Q, 0) || etl::dim(A, 0) != etl::dim(Q, 1)) {
        return false;
    }

    detail::qr_impl::apply(A, Q, R);

    return true;
}

template <etl_expr T, typename G>
void shuffle_flat(T& vector, G&& g) {
    const auto n = etl::size(vector);

    if (n < 2) {
        return;
    }

    if constexpr (impl::egblas::has_shuffle_seed) {
        std::uniform_int_distribution<size_t> seed_dist;

        vector.ensure_gpu_up_to_date();

        impl::egblas::shuffle_seed(n, vector.gpu_memory(), sizeof(etl::value_t<T>), seed_dist(g));

        vector.invalidate_cpu();
    } else {
        using distribution_t = typename std::uniform_int_distribution<size_t>;
        using param_t        = typename distribution_t::param_type;

        distribution_t dist;

        for (auto i = n - 1; i > 0; --i) {
            auto new_i = dist(g, param_t(0, i));

            using std::swap;
            swap(vector[i], vector[new_i]);
        }
    }
}

template <etl_expr T>
void shuffle_flat(T& vector) {
    static std::random_device rd;
    static etl::random_engine g(rd());

    shuffle_flat(vector, g);
}

template <etl_expr T, typename G>
void shuffle_first(T& matrix, G&& g) {
    const auto n = etl::dim<0>(matrix);

    if (n < 2) {
        return;
    }

    if constexpr (impl::egblas::has_shuffle_seed) {
        std::uniform_int_distribution<size_t> seed_dist;

        matrix.ensure_gpu_up_to_date();

        impl::egblas::shuffle_seed(n, matrix.gpu_memory(), sizeof(etl::value_t<T>) * (etl::size(matrix) / n), seed_dist(g));

        matrix.invalidate_cpu();
    } else {
        using distribution_t = typename std::uniform_int_distribution<size_t>;
        using param_t        = typename distribution_t::param_type;

        distribution_t dist;

        auto temp = etl::force_temporary(matrix(0));

        for (auto i = n - 1; i > 0; --i) {
            auto new_i = dist(g, param_t(0, i));

            temp          = matrix(i);
            matrix(i)     = matrix(new_i);
            matrix(new_i) = temp;
        }
    }
}

template <etl_expr T>
void shuffle_first(T& matrix) {
    static std::random_device rd;
    static etl::random_engine g(rd());

    shuffle_first(matrix, g);
}

template <etl_expr T>
void shuffle(T& vector) {
    if constexpr (is_1d<T>) {
        shuffle_flat(vector);
    } else {
        shuffle_first(vector);
    }
}

template <etl_expr T, typename G>
void shuffle(T& vector, G&& g) {
    if constexpr (is_1d<T>) {
        shuffle_flat(vector, g);
    } else {
        shuffle_first(vector, g);
    }
}

template <etl_expr T1, same_dimensions<T1> T2, typename G>
void parallel_shuffle_flat(T1& v1, T2& v2, G&& g) {
    cpp_assert(etl::size(v1) == etl::size(v2), "Impossible to shuffle vector of different dimensions");

    const auto n = etl::size(v1);

    if (n < 2) {
        return;
    }

    if constexpr (impl::egblas::has_par_shuffle_seed) {
        std::uniform_int_distribution<size_t> seed_dist;

        v1.ensure_gpu_up_to_date();
        v2.ensure_gpu_up_to_date();

        impl::egblas::par_shuffle_seed(n, v1.gpu_memory(), sizeof(etl::value_t<T1>), v2.gpu_memory(), sizeof(etl::value_t<T2>), seed_dist(g));

        v1.invalidate_cpu();
        v2.invalidate_cpu();
    } else {
        using distribution_t = typename std::uniform_int_distribution<size_t>;
        using param_t        = typename distribution_t::param_type;

        distribution_t dist;

        for (auto i = n - 1; i > 0; --i) {
            auto new_i = dist(g, param_t(0, i));

            using std::swap;
            swap(v1[i], v1[new_i]);
            swap(v2[i], v2[new_i]);
        }
    }
}

template <etl_expr T1, etl_expr T2>
void parallel_shuffle_flat(T1& v1, T2& v2) {
    static std::random_device rd;
    static etl::random_engine g(rd());

    parallel_shuffle_flat(v1, v2, g);
}

template <etl_expr T1, etl_expr T2, typename G>
void parallel_shuffle_first(T1& m1, T2& m2, G&& g) {
    cpp_assert(etl::dim<0>(m1) == etl::dim<0>(m2), "Impossible to shuffle together matrices of different first dimension");

    const auto n = etl::dim<0>(m1);

    if (n < 2) {
        return;
    }

    if constexpr (impl::egblas::has_par_shuffle_seed) {
        std::uniform_int_distribution<size_t> seed_dist;

        m1.ensure_gpu_up_to_date();
        m2.ensure_gpu_up_to_date();

        impl::egblas::par_shuffle_seed(n, m1.gpu_memory(), sizeof(etl::value_t<T1>) * (etl::size(m1) / n), m2.gpu_memory(),
                                       sizeof(etl::value_t<T2>) * (etl::size(m2) / n), seed_dist(g));

        m1.invalidate_cpu();
        m2.invalidate_cpu();
    } else {
        using distribution_t = typename std::uniform_int_distribution<size_t>;
        using param_t        = typename distribution_t::param_type;

        distribution_t dist;

        auto t1 = etl::force_temporary(m1(0));
        auto t2 = etl::force_temporary(m2(0));

        for (auto i = n - 1; i > 0; --i) {
            auto new_i = dist(g, param_t(0, i));

            t1        = m1(i);
            m1(i)     = m1(new_i);
            m1(new_i) = t1;

            t2        = m2(i);
            m2(i)     = m2(new_i);
            m2(new_i) = t2;
        }
    }
}

template <etl_expr T1, etl_expr T2>
void parallel_shuffle_first(T1& m1, T2& m2) {
    static std::random_device rd;
    static etl::random_engine g(rd());

    parallel_shuffle_first(m1, m2, g);
}

template <etl_expr T1, etl_expr T2>
void parallel_shuffle(T1& v1, T2& v2) {
    static std::random_device rd;
    static etl::random_engine g(rd());

    parallel_shuffle(v1, v2, g);
}

template <etl_expr T>
void shuffle_swap(T& v1, size_t i, size_t new_i) {
    if constexpr (is_1d<T>) {
        auto t    = v1(i);
        v1(i)     = v1(new_i);
        v1(new_i) = t;
    } else {
        auto s1 = v1(i);
        auto s2 = v1(new_i);

        for (size_t index = 0 ; index < etl::subsize(v1); ++index) {
            auto t    = s1[index];
            s1[index] = s2[index];
            s2[index] = t;
        }
    }
}

template <etl_expr T1, etl_expr T2, typename G>
void parallel_shuffle(T1& v1, T2& v2, G&& g) {
    cpp_assert(etl::dim<0>(v1) == etl::dim<0>(v2), "Impossible to shuffle together matrices of different first dimension");

    const auto n = etl::dim<0>(v1);

    if (n < 2) {
        return;
    }

    if constexpr (impl::egblas::has_par_shuffle_seed) {
        std::uniform_int_distribution<size_t> seed_dist;

        v1.ensure_gpu_up_to_date();
        v2.ensure_gpu_up_to_date();

        impl::egblas::par_shuffle_seed(n, v1.gpu_memory(), sizeof(etl::value_t<T1>) * (etl::size(v1) / n), v2.gpu_memory(),
                                       sizeof(etl::value_t<T2>) * (etl::size(v2) / n), seed_dist(g));

        v1.invalidate_cpu();
        v2.invalidate_cpu();
    } else {
        using distribution_t = typename std::uniform_int_distribution<size_t>;
        using param_t        = typename distribution_t::param_type;

        distribution_t dist;

        // Note: We must handle non-homogeneous matrices here
        // For instance, a matrix and a vector hence, the code in shuffle_swap

        for (auto i = n - 1; i > 0; --i) {
            auto new_i = dist(g, param_t(0, i));

            shuffle_swap(v1, i, new_i);
            shuffle_swap(v2, i, new_i);
        }
    }
}

template <etl_expr M, etl_expr N>
M& merge(M& merged, const N& sub, size_t index) {
    const size_t s = etl::size(sub);
    const size_t n = index * s;

    memory_slice(merged, n, n + s) = sub;

    return merged;
}

template <etl_expr M, etl_expr N>
M& batch_merge(M& merged, const N& sub, size_t index) {
    const size_t s = etl::size(sub(0));
    const size_t n = index * s;

    for (size_t b = 0; b < etl::dim<0>(merged); ++b) {
        memory_slice(merged(b), n, n + s) = sub(b);
    }

    return merged;
}

template <etl_expr M, etl_expr N>
M& dispatch(M& dispatched, const N& merged, size_t index) {
    const size_t s = etl::size(dispatched);
    const size_t n = index * s;

    dispatched = memory_slice(merged, n, n + s);

    return dispatched;
}

template <etl_expr M, etl_expr N>
M& batch_dispatch(M& dispatched, const N& merged, size_t index) {
    const size_t s = etl::size(dispatched(0));
    const size_t n = index * s;

    for (size_t b = 0; b < etl::dim<0>(merged); ++b) {
        dispatched(b) = memory_slice(merged(b), n, n + s);
    }

    return dispatched;
}

template <etl_expr M, typename T>
void binarize(M& matrix, T b) {
    using VT = value_t<M>;

    for (auto& value : matrix) {
        value = value > b ? VT(1) : VT(0);
    }
}

template <etl_expr M>
void normalize_flat(M& matrix) {
    using VT = value_t<M>;

    auto m = mean(matrix);

    matrix = matrix - m;

    auto s = stddev(matrix, 0.0);

    if (s != VT(0)) {
        matrix = matrix / s;
    }
}

template <etl_expr M>
void normalize_sub(M& matrix) {
    using VT = value_t<M>;

#ifdef ETL_CUDA
    // The case of CUDA must be handled differently to avoid going
    // back and forth between both memories
    // While this could simply use egblas_normalize_sub, this would
    // be much slower

    auto mm          = matrix.memory_start();
    const auto n     = etl::dim<0>(matrix);
    const auto sub_n = etl::size(matrix) / n;

    for (size_t i = 0; i < n; ++i) {
        VT m(0);

        for (size_t j = 0; j < sub_n; ++j) {
            m += mm[i * sub_n + j];
        }

        m /= VT(sub_n);

        VT s(0);

        for (size_t j = 0; j < sub_n; ++j) {
            s += (mm[i * sub_n + j] - m) * (mm[i * sub_n + j] - m);
        }

        s = std::sqrt(s / VT(sub_n));

        for (size_t j = 0; j < sub_n; ++j) {
            mm[i * sub_n + j] = (mm[i * sub_n + j] - m) / s;
        }
    }
#else
    for (size_t i = 0; i < etl::dim<0>(matrix); ++i) {
        auto m = mean(matrix(i));

        matrix(i) = matrix(i) - m;

        auto s = stddev(matrix(i), 0.0);

        if (s != VT(0)) {
            matrix(i) = matrix(i) / s;
        }
    }
#endif
}

} //end of namespace etl