
.. _program_listing_file__home_david_GitProjects_h5pp_h5pp_include_h5pp_details_h5ppEigen.h:

Program Listing for File h5ppEigen.h
====================================

|exhale_lsh| :ref:`Return to documentation for file <file__home_david_GitProjects_h5pp_h5pp_include_h5pp_details_h5ppEigen.h>` (``/home/david/GitProjects/h5pp/h5pp/include/h5pp/details/h5ppEigen.h``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #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
   
   }
