<commit_before>//
// This file is part of pomerol, an exact diagonalization library aimed at
// solving condensed matter models of interacting fermions.
//
// Copyright (C) 2016-2021 <NAME>, <NAME> and contributors
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

/// \file include/pomerol/MonomialOperator.hpp
/// \brief Storage for an operator that is a product of creation/annihilation operators.
/// \author <NAME> (<EMAIL>)
/// \author <NAME> (<EMAIL>)

#ifndef POMEROL_INCLUDE_MONOMIALOPERATOR_HPP
#define POMEROL_INCLUDE_MONOMIALOPERATOR_HPP

#include "ComputableObject.hpp"
#include "Hamiltonian.hpp"
#include "HilbertSpace.hpp"
#include "Misc.hpp"
#include "MonomialOperatorPart.hpp"
#include "Operators.hpp"
#include "StatesClassification.hpp"

#include "mpi_dispatcher/misc.hpp"
#include "mpi_dispatcher/mpi_skel.hpp"

#include <boost/bimap.hpp>

#include <cstddef>
#include <memory>
#include <stdexcept>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <vector>

namespace Pomerol {

/// \addtogroup ED
///@{

/// A pair of invariant subspace indices.
using BlockMapping = std::pair<BlockNumber, BlockNumber>;

/// \brief Monomial quantum operator.
///
/// This class stores an operator \f$\hat M\f$, which is a monomial, i.e. a product
/// of fermionic and/or bosonic creation/annihilation operators. The operator is stored as
/// a list of matrix blocks (\ref MonomialOperatorPart), each connecting a pair of
/// invariant subspaces of the Hamiltonian. For a given right invariant subspace,
/// there exists at most one part connecting it to a left subspace (and the other way around).
class MonomialOperator : public ComputableObject {
public:
    /// A bi-map container for connections between invariant subspaces established by a monomial operator.
    using BlocksBimap = boost::bimaps::bimap<boost::bimaps::set_of<BlockNumber>, boost::bimaps::set_of<BlockNumber>>;
    /// A single subspace-to-subspace connection established by a monomial operator.
    using BlockMapping = BlocksBimap::value_type;

protected:
    friend class FieldOperatorContainer;

    /// Whether the \ref MOp object is complex-valued.
    bool MOpComplex;
    /// A type-erased real/complex-valued \p libcommute::loperator object.
    std::shared_ptr<void> MOp;

    /// Return a constant reference to the stored \p libcommute::loperator object.
    /// \tparam Complex Request a reference to the complex-valued linear operator.
    /// \pre The compile-time value of \ref Complex must agree with the result of \ref isComplex().
    template <bool Complex> LOperatorTypeRC<Complex> const& getMOp() const {
        assert(MOpComplex == Complex);
        return *std::static_pointer_cast<LOperatorTypeRC<Complex>>(MOp);
    }

    /// Whether the stored parts are complex-valued.
    bool Complex;

    /// Information about invariant subspaces of the Hamiltonian.
    StatesClassification const& S;
    /// The Hamiltonian.
    Hamiltonian const& H;

    /// An automatic space partition object.
    libcommute::space_partition const& Partition;

    /// A map between positions of parts in the \ref parts list and the respective right subspace indices.
    std::unordered_map<std::size_t, BlockNumber> mapPartsFromRight;
    /// A map between positions of parts in the \ref parts list and the respective left subspace indices.
    std::unordered_map<std::size_t, BlockNumber> mapPartsFromLeft;

    /// Left-to-right connections between invariant subspaces established by this monomial operator.
    BlocksBimap LeftRightBlocks;

    /// List of parts (matrix blocks).
    std::vector<MonomialOperatorPart> parts;

public:
    /// Constructor.
    /// \tparam ScalarType Scalar type (either double or std::complex<double>) of the expression \p MO.
    /// \tparam IndexTypes Types of indices carried by operators in the expression \p MO.
    /// \param[in] MO Expression of the monomial operator \f$\hat M\f$.
    /// \param[in] HS Hilbert space.
    /// \param[in] S Information about invariant subspaces of the Hamiltonian.
    /// \param[in] H The Hamiltonian.
    /// \pre \p MO is a monomial operator.
    template <typename ScalarType, typename... IndexTypes>
    MonomialOperator(libcommute::expression<ScalarType, IndexTypes...> const& MO,
                     HilbertSpace<IndexTypes...> const& HS,
                     StatesClassification const& S,
                     Hamiltonian const& H)
        : MOpComplex(std::is_same<ScalarType, ComplexType>::value),
          MOp(std::make_shared<LOperatorType<ScalarType>>(MO, HS.getFullHilbertSpace())),
          Complex(MOpComplex || H.isComplex()),
          S(S),
          H(H),
          Partition(HS.getSpacePartition()) {
        if(MO.size() > 1)
            throw std::runtime_error("Only monomial expressions are supported");
    }

    /// Is the monomial operator a complex-valued matrix?
    bool isComplex() const { return Complex; }

    /// Return a reference to the part by a given left invariant subspace.
    /// \param[in] LeftIndex Index of the left invariant subspace
    /// \pre \ref prepare() has been called.
    MonomialOperatorPart& getPartFromLeftIndex(BlockNumber LeftIndex);
    /// Return a constant reference to the part by a given left invariant subspace.
    /// \param[in] LeftIndex Index of the left invariant subspace
    /// \pre \ref prepare() has been called.
    MonomialOperatorPart const& getPartFromLeftIndex(BlockNumber LeftIndex) const;
    /// Return a reference to the part by a given right invariant subspace.
    /// \param[in] RightIndex Index of the right invariant subspace
    /// \pre \ref prepare() has been called.
    MonomialOperatorPart& getPartFromRightIndex(BlockNumber RightIndex);
    /// Return a constant reference to the part by a given right invariant subspace.
    /// \param[in] RightIndex Index of the right invariant subspace
    /// \pre \ref prepare() has been called.
    MonomialOperatorPart const& getPartFromRightIndex(BlockNumber RightIndex) const;

    /// For a given right invariant subspace, return the corresponding left invariant subspace.
    /// \param[in] RightIndex Index of the right subspace.
    /// \return Index of the left subspace.
    /// \pre \ref prepare() has been called.
    BlockNumber getLeftIndex(BlockNumber RightIndex) const;
    /// For a given left invariant subspace, return the corresponding right invariant subspace.
    /// \param[in] LeftIndex Index of the left subspace.
    /// \return Index of the right subspace.
    /// \pre \ref prepare() has been called.
    BlockNumber getRightIndex(BlockNumber LeftIndex) const;

    /// Return a constant reference to the left-to-right connection map.
    /// \pre \ref prepare() has been called.
    BlocksBimap const& getBlockMapping() const;

    /// Allocate memory for all parts.
    /// \tparam IndexTypes Types of indices carried by operators acting in the Hilbert space \p HS.
    /// \param[in] HS The Hilbert space.
    template <typename... IndexTypes> void prepare(HilbertSpace<IndexTypes...> const& HS) {
        if(getStatus() >= Prepared)
            return;

        auto const& FullHS = HS.getFullHilbertSpace();
        auto Connections = MOpComplex ? Partition.find_connections(getMOp<true>(), FullHS) :
                                        Partition.find_connections(getMOp<false>(), FullHS);

        parts.reserve(Connections.size());
        for(auto const& Conn : Connections) {
            mapPartsFromRight.emplace(Conn.first, parts.size());
            mapPartsFromLeft.emplace(Conn.second, parts.size());
            LeftRightBlocks.insert(BlockMapping(Conn.second, Conn.first));

            if(MOpComplex)
                parts.emplace_back(getMOp<true>(), S, H.getPart(Conn.first), H.getPart(Conn.second));
            else<commit_msg>MonomialOperator: Fix support for unpartitioned Hilbert spaces<commit_after>//
// This file is part of pomerol, an exact diagonalization library aimed at
// solving condensed matter models of interacting fermions.
//
// Copyright (C) 2016-2021 <NAME>, <NAME> and contributors
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

/// \file include/pomerol/MonomialOperator.hpp
/// \brief Storage for an operator that is a product of creation/annihilation operators.
/// \author <NAME> (<EMAIL>)
/// \author <NAME> (<EMAIL>)

#ifndef POMEROL_INCLUDE_MONOMIALOPERATOR_HPP
#define POMEROL_INCLUDE_MONOMIALOPERATOR_HPP

#include "ComputableObject.hpp"
#include "Hamiltonian.hpp"
#include "HilbertSpace.hpp"
#include "Misc.hpp"
#include "MonomialOperatorPart.hpp"
#include "Operators.hpp"
#include "StatesClassification.hpp"

#include "mpi_dispatcher/misc.hpp"
#include "mpi_dispatcher/mpi_skel.hpp"

#include <boost/bimap.hpp>

#include <cstddef>
#include <memory>
#include <stdexcept>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <vector>

namespace Pomerol {

/// \addtogroup ED
///@{

/// A pair of invariant subspace indices.
using BlockMapping = std::pair<BlockNumber, BlockNumber>;

/// \brief Monomial quantum operator.
///
/// This class stores an operator \f$\hat M\f$, which is a monomial, i.e. a product
/// of fermionic and/or bosonic creation/annihilation operators. The operator is stored as
/// a list of matrix blocks (\ref MonomialOperatorPart), each connecting a pair of
/// invariant subspaces of the Hamiltonian. For a given right invariant subspace,
/// there exists at most one part connecting it to a left subspace (and the other way around).
class MonomialOperator : public ComputableObject {
public:
    /// A bi-map container for connections between invariant subspaces established by a monomial operator.
    using BlocksBimap = boost::bimaps::bimap<boost::bimaps::set_of<BlockNumber>, boost::bimaps::set_of<BlockNumber>>;
    /// A single subspace-to-subspace connection established by a monomial operator.
    using BlockMapping = BlocksBimap::value_type;

protected:
    friend class FieldOperatorContainer;

    /// Whether the \ref MOp object is complex-valued.
    bool MOpComplex;
    /// A type-erased real/complex-valued \p libcommute::loperator object.
    std::shared_ptr<void> MOp;

    /// Return a constant reference to the stored \p libcommute::loperator object.
    /// \tparam Complex Request a reference to the complex-valued linear operator.
    /// \pre The compile-time value of \ref Complex must agree with the result of \ref isComplex().
    template <bool Complex> LOperatorTypeRC<Complex> const& getMOp() const {
        assert(MOpComplex == Complex);
        return *std::static_pointer_cast<LOperatorTypeRC<Complex>>(MOp);
    }

    /// Whether the stored parts are complex-valued.
    bool Complex;

    /// Information about invariant subspaces of the Hamiltonian.
    StatesClassification const& S;
    /// The Hamiltonian.
    Hamiltonian const& H;

    /// A map between positions of parts in the \ref parts list and the respective right subspace indices.
    std::unordered_map<std::size_t, BlockNumber> mapPartsFromRight;
    /// A map between positions of parts in the \ref parts list and the respective left subspace indices.
    std::unordered_map<std::size_t, BlockNumber> mapPartsFromLeft;

    /// Left-to-right connections between invariant subspaces established by this monomial operator.
    BlocksBimap LeftRightBlocks;

    /// List of parts (matrix blocks).
    std::vector<MonomialOperatorPart> parts;

public:
    /// Constructor.
    /// \tparam ScalarType Scalar type (either double or std::complex<double>) of the expression \p MO.
    /// \tparam IndexTypes Types of indices carried by operators in the expression \p MO.
    /// \param[in] MO Expression of the monomial operator \f$\hat M\f$.
    /// \param[in] HS Hilbert space.
    /// \param[in] S Information about invariant subspaces of the Hamiltonian.
    /// \param[in] H The Hamiltonian.
    /// \pre \p MO is a monomial operator.
    template <typename ScalarType, typename... IndexTypes>
    MonomialOperator(libcommute::expression<ScalarType, IndexTypes...> const& MO,
                     HilbertSpace<IndexTypes...> const& HS,
                     StatesClassification const& S,
                     Hamiltonian const& H)
        : MOpComplex(std::is_same<ScalarType, ComplexType>::value),
          MOp(std::make_shared<LOperatorType<ScalarType>>(MO, HS.getFullHilbertSpace())),
          Complex(MOpComplex || H.isComplex()),
          S(S),
          H(H) {
        if(MO.size() > 1)
            throw std::runtime_error("Only monomial expressions are supported");
    }

    /// Is the monomial operator a complex-valued matrix?
    bool isComplex() const { return Complex; }

    /// Return a reference to the part by a given left invariant subspace.
    /// \param[in] LeftIndex Index of the left invariant subspace
    /// \pre \ref prepare() has been called.
    MonomialOperatorPart& getPartFromLeftIndex(BlockNumber LeftIndex);
    /// Return a constant reference to the part by a given left invariant subspace.
    /// \param[in] LeftIndex Index of the left invariant subspace
    /// \pre \ref prepare() has been called.
    MonomialOperatorPart const& getPartFromLeftIndex(BlockNumber LeftIndex) const;
    /// Return a reference to the part by a given right invariant subspace.
    /// \param[in] RightIndex Index of the right invariant subspace
    /// \pre \ref prepare() has been called.
    MonomialOperatorPart& getPartFromRightIndex(BlockNumber RightIndex);
    /// Return a constant reference to the part by a given right invariant subspace.
    /// \param[in] RightIndex Index of the right invariant subspace
    /// \pre \ref prepare() has been called.
    MonomialOperatorPart const& getPartFromRightIndex(BlockNumber RightIndex) const;

    /// For a given right invariant subspace, return the corresponding left invariant subspace.
    /// \param[in] RightIndex Index of the right subspace.
    /// \return Index of the left subspace.
    /// \pre \ref prepare() has been called.
    BlockNumber getLeftIndex(BlockNumber RightIndex) const;
    /// For a given left invariant subspace, return the corresponding right invariant subspace.
    /// \param[in] LeftIndex Index of the left subspace.
    /// \return Index of the right subspace.
    /// \pre \ref prepare() has been called.
    BlockNumber getRightIndex(BlockNumber LeftIndex) const;

    /// Return a constant reference to the left-to-right connection map.
    /// \pre \ref prepare() has been called.
    BlocksBimap const& getBlockMapping() const;

    /// Allocate memory for all parts.
    /// \tparam IndexTypes Types of indices carried by operators acting in the Hilbert space \p HS.
    /// \param[in] HS The Hilbert space.
    template <typename... IndexTypes> void prepare(HilbertSpace<IndexTypes...> const& HS) {
        if(getStatus() >= Prepared)
            return;

        if(HS.getStatus() != Computed) { // Hilbert space has not been partitioned
            mapPartsFromRight.emplace(0, 0);
            mapPartsFromLeft.emplace(0, 0);
            LeftRightBlocks.insert(BlockMapping(0, 0));

            if(MOpComplex)
                parts.emplace_back(getMOp<true>(), S, H.getPart(0), H.getPart(0));
            else
                parts.emplace_back(getMOp<false>(), S, H.getPart(0), H.getPart(0));

            Status = Prepared;
            return;
        }

        auto const& FullHS = HS.getFullHilbertSpace();
        auto const& Partition = HS.getSpacePartition();
        auto Connections = MOpComplex ? Partition.find_connections(getMOp<true>(), FullHS) :
                                        Partition.find_connections(getMOp<false>(), FullHS);

        parts.reserve(Connections.size());
        for(auto const& Conn : Connections) {
            mapPartsFromRight.emplace(Conn.first, parts.size());
            mapPartsFromLeft.emplace(Conn.second, parts.size());
            LeftRightBlocks.insert(BlockMapping(Conn.second, Conn.first));

            if(MOpComplex)
                parts.emplace_back(getMOp<true>(), S, H.getPart(Conn.first), H.getPart(Conn.second));
            else
{'fraction_non_alphanumeric': 0.10553284938693344, 'fraction_numerical': 0.002196059293600927, 'mean_word_length': 3.078109452736318, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 85, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3206592', 'n_tokens_mistral': 4673, 'n_tokens_neox': 4428, 'n_words': 1577}