Program Listing for File h5ppEigen.h

Return to documentation for file (/home/david/GitProjects/h5pp/h5pp/include/h5pp/details/h5ppEigen.h)

#pragma once

#if __has_include(<Eigen/Core>)
    #define H5PP_EIGEN3
    #include <Eigen/Core>
    #include <unsupported/Eigen/CXX11/Tensor>
#endif

#include <iterator>
#include <numeric>
#if !defined(_MSC_VER)
    #define H5PP_CONSTEXPR constexpr
#else
    #define H5PP_CONSTEXPR
#endif

namespace h5pp {
#ifdef H5PP_EIGEN3
    namespace eigen {
        template<typename Scalar>
        using MatrixType = Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
        template<typename Scalar>
        using VectorType = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
        template<auto rank>
        using array = Eigen::array<Eigen::Index, rank>;
        template<auto rank>
        using DSizes = Eigen::DSizes<Eigen::Index, rank>;
        using array8 = array<8>;
        using array7 = array<7>;
        using array6 = array<6>;
        using array5 = array<5>;
        using array4 = array<4>;
        using array3 = array<3>;
        using array2 = array<2>;
        using array1 = array<1>;

        // Handy functions to copy lists of dimensions
        template<typename T, auto rank, typename Container>
        void copy_dims(Eigen::DSizes<T, rank> dsizes, const Container &container) {
            if constexpr(std::is_array_v<Container>) {
                std::copy_n(std::begin(container), rank, dsizes.begin());
            } else {
                if(container.size() != rank) throw std::runtime_error("copy_dims: Wrong container size, can't copy dimensions.");
                std::copy(std::begin(container), std::end(container), dsizes.begin());
            }
        }
        template<auto rank, typename Container>
        [[nodiscard]] H5PP_CONSTEXPR Eigen::DSizes<Eigen::Index, rank> copy_dims(const Container &container) {
            Eigen::DSizes<Eigen::Index, rank> dsizes;
            if constexpr(std::is_array_v<Container>) {
                std::copy_n(std::begin(container), rank, dsizes.begin());
            } else {
                if(container.size() != rank) throw std::runtime_error("copy_dims: Wrong container size, can't copy dimensions.");
                std::copy(std::begin(container), std::end(container), dsizes.begin());
            }
            return dsizes;
        }

        template<auto rank, typename Container>
        [[nodiscard]] H5PP_CONSTEXPR Eigen::DSizes<Eigen::Index, rank> copy_dims(const Container *container) {
            Eigen::DSizes<Eigen::Index, rank> dsizes;
            std::copy_n(container, rank, dsizes.begin());
            return dsizes;
        }

        // Shorthand for the list of index pairs.
        template<auto N>
        using idxlistpair = Eigen::array<Eigen::IndexPair<Eigen::Index>, N>;

        inline H5PP_CONSTEXPR idxlistpair<0> idx() { return {}; }

        template<std::size_t N, typename idxType>
        H5PP_CONSTEXPR idxlistpair<N> idx(const idxType (&list1)[N], const idxType (&list2)[N]) {
            // Use numpy-style indexing for contraction. Each list contains a list of indices to be contracted for the respective
            // tensors. This function zips them together into pairs as used in Eigen::Tensor module. This does not sort the indices in decreasing order.
            static_assert(std::is_integral_v<idxType>);
            idxlistpair<N> pairlistOut;
            for(size_t i = 0; i < N; i++) { pairlistOut[i] = Eigen::IndexPair<Eigen::Index>{static_cast<Eigen::Index>(list1[i]), static_cast<Eigen::Index>(list2[i])}; }
            return pairlistOut;
        }

        struct idx_dim_pair {
            Eigen::Index idxA;
            Eigen::Index idxB;
            Eigen::Index dimB;
        };

        template<std::size_t NB, std::size_t N>
        H5PP_CONSTEXPR idxlistpair<N> sortIdx(const Eigen::array<Eigen::Index, NB> &dimensions, const Eigen::Index (&idx_ctrct_A)[N], const Eigen::Index (&idx_ctrct_B)[N]) {
            // When doing contractions, some indices may be larger than others. For performance, you want to
            // contract the largest indices first. This will return a sorted index list in decreasing order.
            Eigen::array<idx_dim_pair, N> idx_dim_pair_list;
            for(size_t i = 0; i < N; i++) { idx_dim_pair_list[i] = {idx_ctrct_A[i], idx_ctrct_B[i], dimensions[idx_ctrct_B[i]]}; }
            std::sort(idx_dim_pair_list.begin(), idx_dim_pair_list.end(), [](const auto &i, const auto &j) { return i.dimB > j.dimB; });
            idxlistpair<N> pairlistOut;
            for(size_t i = 0; i < N; i++) { pairlistOut[i] = Eigen::IndexPair<long>{idx_dim_pair_list[i].idxA, idx_dim_pair_list[i].idxB}; }
            return pairlistOut;
        }

        //
        //    //***************************************//
        //    //Different views for rank 1 and 2 tensors//
        //    //***************************************//
        //

        template<typename Scalar>
        constexpr Eigen::Tensor<Scalar, 1> extractDiagonal(const Eigen::Tensor<Scalar, 2> &tensor) {
            auto rows = tensor.dimension(0);
            auto cols = tensor.dimension(1);
            if(tensor.dimension(0) != tensor.dimension(1)) throw std::runtime_error("extractDiagonal expects a square tensor");

            Eigen::Tensor<Scalar, 1> diagonals(rows);
            for(auto i = 0; i < rows; i++) { diagonals(i) = tensor(i, i); }
            return diagonals;
        }

        template<typename Scalar>
        constexpr auto asDiagonal(const Eigen::Tensor<Scalar, 1> &tensor) {
            return tensor.inflate(array1{tensor.size() + 1}).reshape(array2{tensor.size(), tensor.size()});
        }

        template<typename Scalar>
        constexpr auto asDiagonalSquared(const Eigen::Tensor<Scalar, 1> &tensor) {
            return tensor.square().inflate(array1{tensor.size() + 1}).reshape(array2{tensor.size(), tensor.size()});
        }

        template<typename Scalar>
        constexpr auto asDiagonalInversed(const Eigen::Tensor<Scalar, 1> &tensor) {
            return tensor.inverse().inflate(array1{tensor.size() + 1}).reshape(array2{tensor.size(), tensor.size()});
        }

        template<typename Scalar>
        constexpr auto asDiagonalInversed(const Eigen::Tensor<Scalar, 2> &tensor) {
            if(tensor.dimension(0) != tensor.dimension(1)) throw std::runtime_error("Textra::asDiagonalInversed expects a square tensor");
            Eigen::Tensor<Scalar, 2> inversed = asDiagonalInversed(extractDiagonal(tensor));
            return inversed;
        }

        template<typename Scalar>
        constexpr auto asNormalized(const Eigen::Tensor<Scalar, 1> &tensor) {
            Eigen::Map<const VectorType<Scalar>> map(tensor.data(), tensor.size());
            return Eigen::TensorMap<Eigen::Tensor<const Scalar, 1>>(map.normalized().eval().data(), array1{map.size()});
        }

        //    //****************************//
        //    //Matrix to tensor conversions//
        //    //****************************//

        // Detects if Derived is a plain object, like "MatrixXd" or similar.
        // std::decay removes pointer or ref qualifiers if present
        template<typename Derived>
        using is_plainObject = std::is_base_of<Eigen::PlainObjectBase<std::decay_t<Derived>>, std::decay_t<Derived>>;
        template<typename Derived>
        using is_matrixObject = std::is_base_of<Eigen::MatrixBase<std::decay_t<Derived>>, std::decay_t<Derived>>;
        template<typename Derived>
        using is_arrayObject = std::is_base_of<Eigen::ArrayBase<std::decay_t<Derived>>, std::decay_t<Derived>>;

        template<typename Derived, auto rank>
        constexpr Eigen::Tensor<typename Derived::Scalar, rank> Matrix_to_Tensor(const Eigen::EigenBase<Derived> &matrix, const array<rank> &dims) {
            if constexpr(is_plainObject<Derived>::value) {
                // Return map from raw input.
                return Eigen::TensorMap<const Eigen::Tensor<const typename Derived::Scalar, rank>>(matrix.derived().eval().data(), dims);
            } else {
                // Create a temporary
                MatrixType<typename Derived::Scalar> matref = matrix;
                return Eigen::TensorMap<Eigen::Tensor<typename Derived::Scalar, rank>>(matref.data(), dims);
            }
        }

        // Helpful overload
        template<typename Derived, typename... Dims>
        constexpr Eigen::Tensor<typename Derived::Scalar, sizeof...(Dims)> Matrix_to_Tensor(const Eigen::EigenBase<Derived> &matrix, const Dims... dims) {
            return Matrix_to_Tensor(matrix, array<sizeof...(Dims)>{dims...});
        }
        // Helpful overload
        template<typename Derived, auto rank>
        constexpr Eigen::Tensor<typename Derived::Scalar, rank> Matrix_to_Tensor(const Eigen::EigenBase<Derived> &matrix, const DSizes<rank> &dims) {
            array<rank> dim_array = dims;
            std::copy(std::begin(dims), std::end(dims), std::begin(dim_array));
            return Matrix_to_Tensor(matrix, dim_array);
        }

        template<typename Derived>
        constexpr auto Matrix_to_Tensor1(const Eigen::EigenBase<Derived> &matrix) {
            return Matrix_to_Tensor(matrix, matrix.size());
        }
        template<typename Derived>
        constexpr auto Matrix_to_Tensor2(const Eigen::EigenBase<Derived> &matrix) {
            return Matrix_to_Tensor(matrix, matrix.rows(), matrix.cols());
        }

        //****************************//
        // Tensor to matrix conversions//
        //****************************//

        template<typename Scalar>
        constexpr MatrixType<Scalar> Tensor2_to_Matrix(const Eigen::Tensor<Scalar, 2> &tensor) {
            return Eigen::Map<const MatrixType<Scalar>>(tensor.data(), tensor.dimension(0), tensor.dimension(1));
        }

        template<typename Scalar>
        constexpr MatrixType<Scalar> Tensor1_to_Vector(const Eigen::Tensor<Scalar, 1> &tensor) {
            return Eigen::Map<const VectorType<Scalar>>(tensor.data(), tensor.size());
        }

        template<typename Scalar, auto rank, typename sizeType>
        constexpr MatrixType<Scalar> Tensor_to_Matrix(const Eigen::Tensor<Scalar, rank> &tensor, const sizeType rows, const sizeType cols) {
            return Eigen::Map<const MatrixType<Scalar>>(tensor.data(), rows, cols);
        }

        //************************//
        // change storage layout //
        //************************//
        template<typename Derived>
        auto to_RowMajor(const Eigen::TensorBase<Derived, Eigen::ReadOnlyAccessors> &tensor) {
            if constexpr(Eigen::RowMajor == static_cast<Eigen::StorageOptions>(Derived::Layout))
                return tensor;
            else {
                array<Derived::NumIndices> neworder;
                std::iota(std::begin(neworder), std::end(neworder), 0);
                std::reverse(neworder.data(), neworder.data() + neworder.size());
                return Eigen::Tensor<typename Derived::Scalar, Derived::NumIndices, Eigen::RowMajor>(tensor.swap_layout().shuffle(neworder));
            }
        }

        template<typename Derived>
        auto to_ColMajor(const Eigen::TensorBase<Derived, Eigen::ReadOnlyAccessors> &tensor) {
            if constexpr(Eigen::ColMajor == static_cast<Eigen::StorageOptions>(Derived::Layout))
                return tensor;
            else {
                array<Derived::NumIndices> neworder;
                std::iota(std::begin(neworder), std::end(neworder), 0);
                std::reverse(neworder.data(), neworder.data() + neworder.size());
                return Eigen::Tensor<typename Derived::Scalar, Derived::NumIndices, Eigen::ColMajor>(tensor.swap_layout().shuffle(neworder));
            }
        }

        template<typename Derived>
        auto to_RowMajor(const Eigen::DenseBase<Derived> &dense) {
            if constexpr(Derived::IsRowMajor) { return dense; }
            if constexpr(is_matrixObject<Derived>::value) {
                return Eigen::Matrix<typename Derived::Scalar, Derived::RowsAtCompileTime, Derived::ColsAtCompileTime, Eigen::RowMajor>(dense);
            } else if constexpr(is_arrayObject<Derived>::value) {
                return Eigen::Array<typename Derived::Scalar, Derived::RowsAtCompileTime, Derived::ColsAtCompileTime, Eigen::RowMajor>(dense);
            }
            throw std::runtime_error("Wrong dense type?? Report this bug!");
        }

        template<typename Derived>
        auto to_ColMajor(const Eigen::DenseBase<Derived> &dense) {
            if constexpr(not Derived::IsRowMajor) { return dense; }
            if constexpr(is_matrixObject<Derived>::value) {
                return Eigen::Matrix<typename Derived::Scalar, Derived::RowsAtCompileTime, Derived::ColsAtCompileTime, Eigen::ColMajor>(dense);
            } else if constexpr(is_arrayObject<Derived>::value) {
                return Eigen::Array<typename Derived::Scalar, Derived::RowsAtCompileTime, Derived::ColsAtCompileTime, Eigen::ColMajor>(dense);
            }
            throw std::runtime_error("Wrong dense type?? Report this bug!");
        }
    }

#endif

}