Tensor.hpp

Go to the documentation of this file.

// This file is part of GQCG-GQCP.
//
// Copyright (C) 2017-2020  the GQCG developers
//
// GQCG-GQCP is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// GQCG-GQCP is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with GQCG-GQCP.  If not, see <http://www.gnu.org/licenses/>.
 
#pragma once
 
 
#include "Mathematical/Representation/Matrix.hpp"
#include "Utilities/type_traits.hpp"
 
#include <boost/algorithm/string.hpp>
 
#include <unsupported/Eigen/CXX11/Tensor>
 
#include <algorithm>
#include <iostream>
#include <string>
 
 
namespace GQCP {
 
 
template <typename _Scalar, int _Rank>
class Tensor:
    public Eigen::Tensor<_Scalar, _Rank> {
 
public:
    // The scalar type of one of the elements of the tensor.
    using Scalar = _Scalar;
 
    // The rank of the tensor, i.e. the number of axes.
    static constexpr auto Rank = _Rank;
 
    // The type of 'this'.
    using Self = Tensor<Scalar, Rank>;
 
    // The Eigen base type.
    using Base = Eigen::Tensor<Scalar, Rank>;
 
 
public:
    /*
     *  MARK: Constructors
     */
 
    // Inherit `Eigen::Tensor`'s constructors.
    using Eigen::Tensor<Scalar, _Rank>::Tensor;
 
 
    /*
     *  MARK: Named constructors
     */
 
    template <int Z = Rank>
    static enable_if_t<Z == 4, Self> FromBlock(const Self& T, const size_t i, const size_t j, const size_t k, const size_t l, const size_t cutoff = 0) {
 
        Tensor<double, Rank> T_result {static_cast<long>(T.dimension(0) - i - cutoff),
                                       static_cast<long>(T.dimension(1) - j - cutoff),
                                       static_cast<long>(T.dimension(2) - k - cutoff),
                                       static_cast<long>(T.dimension(3) - l - cutoff)};
        T_result.setZero();
 
        for (size_t p = 0; p < T_result.dimension(0); p++) {
            for (size_t q = 0; q < T_result.dimension(1); q++) {
                for (size_t r = 0; r < T_result.dimension(2); r++) {
                    for (size_t s = 0; s < T_result.dimension(3); s++) {
                        T_result(p, q, r, s) = T(i + p, j + q, k + r, l + s);
                    }
                }
            }
        }
 
        return T_result;
    }
 
 
    /*
     *  MARK: Conversions
     */
 
    const Base& Eigen() const { return static_cast<const Base&>(*this); }
 
    Base& Eigen() { return static_cast<Base&>(*this); }
 
 
    /*
     *  MARK: Tensor operations
     */
 
    template <size_t r, size_t s, int Z = Rank>
    enable_if_t<Z == 4, Self&> addBlock(const MatrixX<Scalar>& M, const size_t i, const size_t j, const size_t k, const size_t l) {
 
        // Initialize series of arrays with 1 or 0 values, so that the correct tensor indices given by the template argument correspond to the matrix indices
        size_t ia[4] = {1, 0, 0, 0};
        size_t ja[4] = {0, 1, 0, 0};
        size_t ka[4] = {0, 0, 1, 0};
        size_t la[4] = {0, 0, 0, 1};
 
        for (size_t x = 0; x < M.rows(); x++) {
            for (size_t y = 0; y < M.cols(); y++) {
 
                size_t i_effective = i + x * ia[r] + y * ia[s];
                size_t j_effective = j + x * ja[r] + y * ja[s];
                size_t k_effective = k + x * ka[r] + y * ka[s];
                size_t l_effective = l + x * la[r] + y * la[s];
 
                this->operator()(i_effective, j_effective, k_effective, l_effective) += M(x, y);
            }
        }
 
        return (*this);
    }
 
 
    template <int Z = Rank>
    enable_if_t<Z == 4, Self&> addBlock(const Self& T, const size_t i, const size_t j, const size_t k, const size_t l) {
 
        for (size_t p = 0; p < T.dimension(0); p++) {
            for (size_t q = 0; q < T.dimension(1); q++) {
                for (size_t r = 0; r < T.dimension(2); r++) {
                    for (size_t s = 0; s < T.dimension(3); s++) {
                        this->operator()(i + p, j + q, k + r, l + s) += T(p, q, r, s);
                    }
                }
            }
        }
 
        return (*this);
    }
 
 
    template <int Z = Rank>
    enable_if_t<Z == 4, Matrix<Scalar>> pairWiseReduced(const size_t start_i = 0, const size_t start_j = 0, const size_t start_k = 0, const size_t start_l = 0) const {
 
        // Initialize the resulting matrix
        const auto dims = this->dimensions();
        Matrix<Scalar> M {(dims[0] - start_i) * (dims[1] - start_j),
                          (dims[2] - start_k) * (dims[3] - start_l)};
 
        // Calculate the compound indices and bring the elements from the tensor over into the matrix
        size_t row_index = 0;
        for (size_t j = start_j; j < dims[1]; j++) {      // "column major" ordering for row_index<-i,j so we do j first, then i
            for (size_t i = start_i; i < dims[0]; i++) {  // in column major indices, columns are contiguous, so the first of two indices changes more rapidly
 
                size_t column_index = 0;
                for (size_t l = start_l; l < dims[3]; l++) {      // "column major" ordering for column_index<-k,l so we do l first, then k
                    for (size_t k = start_k; k < dims[2]; k++) {  // in column major indices, columns are contiguous, so the first of two indices changes more rapidly
 
                        M(row_index, column_index) = this->operator()(i, j, k, l);
 
                        column_index++;
                    }
                }
 
                row_index++;
            }
        }
 
        return M;
    }
 
 
    template <int N, int LHSRank = Rank, int RHSRank>
    Tensor<Scalar, LHSRank + RHSRank - 2 * N> einsum(const Tensor<Scalar, RHSRank>& rhs, const std::string& lhs_labels, const std::string& rhs_labels, const std::string& output_labels) const {
 
        // Check the length of the indices.
        constexpr auto ResultRank = LHSRank + RHSRank - 2 * N;
 
        if (lhs_labels.size() != LHSRank) {
            throw std::invalid_argument("Tensor.einsum(const Tensor<Scalar, RHSRank>&, const std::string&, const std::string&, const std::string&): The number of indices for the left-hand side of the contraction does not match the rank of the left-hand side tensor.");
        }
 
        if (rhs_labels.size() != RHSRank) {
            throw std::invalid_argument("Tensor.einsum(const Tensor<Scalar, RHSRank>&, const std::string&, const std::string&, const std::string&): The number of indices for the right-hand side of the contraction does not match the rank of the right-hand side tensor.");
        }
 
        if (output_labels.size() != ResultRank) {
            throw std::invalid_argument("Tensor.einsum(const Tensor<Scalar, RHSRank>&, const std::string&, const std::string&, const std::string&): The number of output indices does not match the number of axes that should be contracted over.");
        }
 
 
        // Find the indices that should be contracted over, these are the positions of the axis labels that are both in the lhs- and rhs-indices.
        Eigen::array<Eigen::IndexPair<int>, N> contraction_pairs {};
        size_t array_position = 0;              // The index at which an `Eigen::IndexPair` should be placed.
        for (size_t i = 0; i < LHSRank; i++) {  // 'i' loops over the lhs-indices
            const auto current_label = lhs_labels[i];
            const auto match = rhs_labels.find(current_label);
            if (match != std::string::npos) {  // an actual match was found
                contraction_pairs[array_position] = Eigen::IndexPair<int>(i, match);
                array_position++;
            }
        }
 
        // Perform the contraction. It is an intermediate result, because we still have to align the tensor axes that Eigen's contraction module produces with the user's requested axes.
        Tensor<Scalar, ResultRank> T_intermediate = this->contract(rhs.Eigen(), contraction_pairs);
 
        // Finally, we should find the shuffle indices that map the obtained intermediate axes to the requested axes.
        // The intermediate axes are just concatenated from left to right, removing any duplicates.
        auto intermediate_indices = lhs_labels + rhs_labels;
 
        for (size_t i = 1; i < LHSRank + RHSRank; i++) {  // We should skip the first iteration, since there wouldn't be any duplicates anyways.
 
            // Check if the current index label appears in the substring before it.
            const auto current_label = intermediate_indices[i];
            const auto match = intermediate_indices.substr(0, i).find(current_label);
 
            if (match != std::string::npos) {  // an actual match was found
                intermediate_indices.erase(match, 1);
                intermediate_indices.erase(i - 1, 1);  // i-1 is used because the length of the parameter ìntermediate_indices`changed.
 
                // We have to set back the current index to account for the removed indices.
                i -= 2;
            }
        }
 
        // Eigen3 does not document its tensor contraction clearly, so see the accepted answer on stackoverflow (https://stackoverflow.com/a/47558349/7930415):
        //      Eigen3 does not accept a way to specify the output axes: instead, it retains the order from left to right of the axes that survive the contraction.
        //      This means that, in order to get the right ordering of the axes, we will have to swap axes.
 
        // Find the position of the intermediate axes' labels in the requested output labels in order to set up the required shuffle indices.
        // This is only necessary when not contracting over all axes, in other words, when the string of output labels is not empty.
        if (output_labels != "") {
            Eigen::array<int, 4> shuffle_indices {};
 
            for (size_t i = 0; i < 4; i++) {
                const auto current_label = intermediate_indices[i];
                shuffle_indices[i] = output_labels.find(current_label);
            }
 
            return T_intermediate.shuffle(shuffle_indices);
        }
 
        return T_intermediate;
    }
 
 
    template <int N, int LHSRank = Rank, int RHSRank>
    Tensor<Scalar, LHSRank + RHSRank - 2 * N> einsum(std::string contraction_string, const Tensor<Scalar, RHSRank>& rhs) const {
        // Remove unnecessary symbols from string.
        boost::erase_all(contraction_string, " ");  // Remove all spaces.
        boost::erase_all(contraction_string, ">");
        boost::replace_all(contraction_string, "-", " ");
        boost::replace_all(contraction_string, ",", " ");
 
        // Split the stringstream in the necessary components, 3 in most cases.
        std::vector<std::string> segment_list;
        boost::split(segment_list, contraction_string, boost::is_any_of(" "));
 
        // If a contraction over all indices is done, only 2 sets of labels are saved in the vector. The last set of labels is an empty string, and is added manually.
        if (segment_list.size() == 2) {
            segment_list.push_back("");
        }
 
        // The following code can be used to calculate template parameter `N` at runtime. At the moment however, this is not used.
        // Determine number of axes to contract over
        //
        // std::vector<char> segment_1(segment_list[0].begin(), segment_list[0].end());
        // std::vector<char> segment_2(segment_list[1].begin(), segment_list[1].end());
        // std::vector<char> intersection;
        //
        // std::set_intersection(segment_1.begin(), segment_1.end(),
        //                       segment_2.begin(), segment_2.end(), std::back_inserter(intersection));
        //
        // const int Z = intersection.size();
 
        return this->einsum<N>(rhs, segment_list[0], segment_list[1], segment_list[2]);
    }
 
 
    template <int N, int LHSRank = Rank>
    Tensor<Scalar, LHSRank + 2 - 2 * N> einsum(std::string contraction_string, const Matrix<Scalar> rhs) const {
        // Convert the given `Matrix` to its equivalent rank-2 `Tensor`.
        const auto tensor_from_matrix = Tensor<Scalar, 2>(Eigen::TensorMap<Eigen::Tensor<const Scalar, 2>>(rhs.data(), rhs.rows(), rhs.cols()));
        return this->einsum<N>(contraction_string, tensor_from_matrix);
    }
 
 
    template <int Z = Rank>
    const enable_if_t<Z == 2, GQCP::Matrix<Scalar>> asMatrix() const {
 
        const auto rows = this->dimension(0);
        const auto cols = this->dimension(1);
 
        return GQCP::Matrix<Scalar>(Eigen::Map<const Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>>(this->data(), rows, cols));
    }
 
    GQCP::Matrix<Scalar> reshape(const size_t rows, const size_t cols) const {
 
        // The dimensions of the tensor have to match the dimensions of the target matrix.
        if (rows * cols != this->numberOfElements()) {
            throw std::invalid_argument("Tensor.reshape(const size_t rows, const size_t cols): The tensor cannot be reshaped into a matrix with the specified dimensions.");
        }
 
        // Initialize a matrix with the wanted dimensions.
        GQCP::Matrix<Scalar> M {rows, cols};
 
        // We create an intermediate vector to contain all the values of the tensor.
        std::vector<Scalar> vectorized;
 
        // We reserve the amount of data we expect the tensor to take up.
        vectorized.reserve(rows * cols);
 
        // Fill the vector with the elements of the tensor.
        for (int i = 0; i < this->dimension(0); i++) {
            for (int j = 0; j < this->dimension(1); j++) {
                for (int k = 0; k < this->dimension(2); k++) {
                    for (int l = 0; l < this->dimension(3); l++) {
                        vectorized.push_back(this->operator()(i, j, k, l));
                    }
                }
            }
        }
 
        // Fill the matrix with the elements of the intermediate vector.
        for (int r = 0; r < rows; r++) {
            for (int c = 0; c < cols; c++) {
 
                // The columns are looped first. Hence, the position within the row is determined by c (+c).
                // The rows are determined by the outer loop. When you start filling a new row, you have to skip all the vector elements used for the previous row (r * cols).
                M(r, c) = vectorized[r * cols + c];
            }
        }
 
        return M;
    }
 
 
    /*
     *  MARK: General information
     */
 
    bool hasEqualDimensionsAs(const Self& other) const {
 
        for (size_t i = 0; i < Rank; i++) {
            if (this->dimension(i) != other.dimension(i)) {
                return false;
            }
        }
 
        return true;
    }
 
 
    template <int Z = Rank>
    enable_if_t<Z == 4, bool> isApprox(const Self& other, const double tolerance = 1.0e-12) const {
 
        if (!this->hasEqualDimensionsAs(other)) {
            throw std::invalid_argument("RankFourTensor<Scalar>::isApprox(Self, double): the tensors have different dimensions");
        }
 
        // Check every pair of values
        for (size_t i = 0; i < this->dimension(0); i++) {
            for (size_t j = 0; j < this->dimension(1); j++) {
                for (size_t k = 0; k < this->dimension(2); k++) {
                    for (size_t l = 0; l < this->dimension(3); l++) {
                        if (std::abs(this->operator()(i, j, k, l) - other(i, j, k, l)) > tolerance) {
                            return false;
                        }
                    }
                }
            }
        }  // rank-4 tensor traversing
 
        return true;
    }
 
 
    size_t numberOfElements() const {
        return this->dimension(0) * this->dimension(1) * this->dimension(2) * this->dimension(3);
    }
 
 
    template <int Z = Rank>
    enable_if_t<Z == 4> print(std::ostream& output_stream = std::cout) const {
 
        for (size_t i = 0; i < this->dimension(0); i++) {
            for (size_t j = 0; j < this->dimension(1); j++) {
                for (size_t k = 0; k < this->dimension(2); k++) {
                    for (size_t l = 0; l < this->dimension(3); l++) {
                        output_stream << i << ' ' << j << ' ' << k << ' ' << l << "  " << this->operator()(i, j, k, l) << std::endl;
                    }
                }
            }
        }
    }
};
 
 
}  // namespace GQCP