Ensemble.hpp

Go to the documentation of this file.

#pragma once
#ifdef GPU
    #include "Ensemble.cuh"
#else
    #include "Headers/mytypes.hpp"
    #include "Headers/HilbertSpace/TransSector.hpp"
    #include "Headers/HilbertSpace/OperatorSpace.hpp"
    #include "Headers/HilbertSpace/RandomMatrix.hpp"
    #include <Eigen/Dense>
    #include <Eigen/Sparse>
    #include <unsupported/Eigen/KroneckerProduct>
 
int g_dimLoc   = 2;
int g_momentum = 0;
 
using RealScalar = Real_t;
using Scalar     = Complex_t<RealScalar>;
 
class SubSpaceManager {
    public:
        using SubHilbertSpace = TransSector<RealScalar>;
        using TotalSpace      = typename SubHilbertSpace::TotalSpace;
 
    private:
        std::vector<SubHilbertSpace> m_sectors;
        std::vector<bool>            m_isConstructed;
        SubSpaceManager() = default;
 
    public:
        SubSpaceManager(const SubSpaceManager&)              = delete;
        SubSpaceManager& operator=(const SubSpaceManager&)   = delete;
        SubSpaceManager(SubSpaceManager&&)                   = delete;
        SubSpaceManager&        operator=(SubSpaceManager&&) = delete;
        static SubSpaceManager& get_instance() {
            static SubSpaceManager instance;
            return instance;
        }
 
        SubHilbertSpace& operator[](size_t L) {
            SubSpaceManager& instance = this->get_instance();
            if(instance.m_sectors.size() <= L) {
                instance.m_sectors.resize(L + 1);
                instance.m_isConstructed.resize(L + 1, false);
 
                if(!instance.m_isConstructed[L]) {
                    instance.m_sectors[L]
                        = TransSector<RealScalar>(g_momentum, ManyBodySpinSpace(L, g_dimLoc));
                    instance.m_isConstructed[L] = true;
                }
            }
            return instance.m_sectors[L];
        }
};
 
namespace Ensemble::internal {
    template<class MatrixType>
    class ScalableOperator {
        private:
            using Scalar = typename MatrixType::Scalar;
            MatrixType m_locOp;
            size_t     m_dimLoc;
            size_t     m_Locality;
 
        public:
            ScalableOperator(MatrixType&& locOp, size_t dimLoc) : m_locOp{locOp}, m_dimLoc{dimLoc} {
                if(m_locOp.rows() != m_locOp.cols()) {
                    std::cerr << "Error at [" << __FILE__ << ":" << __LINE__ << "]\n\t"
                              << __PRETTY_FUNCTION__ << "\nMessage:\t locOp.rows()("
                              << m_locOp.rows() << ") != locOp.cols()(" << locOp.cols() << "): "
                              << "local operator must be a square matrix" << std::endl;
                    std::exit(EXIT_FAILURE);
                }
                size_t dim = m_locOp.rows();
                for(m_Locality = 0; dim != 1; ++m_Locality) {
                    if(dim % dimLoc != 0) {
                        std::cerr << "Error at [" << __FILE__ << ":" << __LINE__ << "]\n\t"
                                  << __PRETTY_FUNCTION__ << "\nMessage:\t locOp.rows()("
                                  << m_locOp.rows() << ") must be a power of dimLoc(" << dimLoc
                                  << ")" << std::endl;
                        std::exit(EXIT_FAILURE);
                    }
                    dim /= dimLoc;
                }
            }
 
            MatrixType construct_globalOp(size_t L) {
                if(L < m_Locality) {
                    std::cerr << "Error at [" << __FILE__ << ":" << __LINE__ << "]\n\t"
                              << __PRETTY_FUNCTION__ << "\n\t L(" << L
                              << ") must not be smaller than m_Locality(" << m_Locality << ")."
                              << std::endl;
                    std::exit(EXIT_FAILURE);
                }
 
                size_t dimCompl = 1;
                for(int j = 0; j < L - m_Locality; ++j) { dimCompl *= g_dimLoc; }
 
                auto             Identity     = MatrixType::Identity(dimCompl, dimCompl);
                SubSpaceManager& subSpaces    = SubSpaceManager::get_instance();
                auto&            subSpaceBase = subSpaces[L].basis();
                MatrixType res = subSpaceBase.adjoint() * Eigen::kroneckerProduct(Identity, m_locOp)
                                 * subSpaceBase;
                return res;
            };
    };
}  // namespace Ensemble::internal
 
class HamiltonianEnsemble {
    public:
        using MatrixType       = Eigen::MatrixX<Scalar>;
        using EigenSolver      = Eigen::SelfAdjointEigenSolver<MatrixType>;
        using StateSpace       = typename SubSpaceManager::SubHilbertSpace;
        using ScalableOperator = Ensemble::internal::ScalableOperator<MatrixType>;
 
        static std::string dataname(size_t rep) { return "Hamiltonian_No" + std::to_string(rep); };
        static std::string usage(size_t nargs) {
            std::string usage = "HamiltonianEnsemble Usage: ";
            usage += std::to_string(nargs + 0) + ".(Locality) ";
            return usage;
        };
        static constexpr size_t nargs(void) { return 1; };
 
        static StateSpace& stateSpace(size_t L) {
            SubSpaceManager& subSpaces = SubSpaceManager::get_instance();
            return subSpaces[L];
        }
 
    private:
        size_t                                                 m_Locality;
        GaussianRandomMatrixGenerator< OperatorSpace<Scalar> > m_GRME;
 
        HamiltonianEnsemble(char** CL_argv) : m_Locality(std::atoi(CL_argv[0])) {
            assert(m_Locality >= 1);
            size_t dim = 1;
            for(int j = 0; j < m_Locality; ++j) { dim *= g_dimLoc; };
            auto HSpace  = HilbertSpace<int>(dim);
            auto OpSpace = OperatorSpace<Scalar>(HSpace);
            m_GRME       = GaussianRandomMatrixGenerator< decltype(OpSpace) >(OpSpace, 0);
        }
 
    public:
        HamiltonianEnsemble(const HamiltonianEnsemble&)              = delete;
        HamiltonianEnsemble& operator=(const HamiltonianEnsemble&)   = delete;
        HamiltonianEnsemble(HamiltonianEnsemble&&)                   = delete;
        HamiltonianEnsemble&        operator=(HamiltonianEnsemble&&) = delete;
        static HamiltonianEnsemble& get_instance(char** CL_argv) {
            static HamiltonianEnsemble instance(CL_argv);
            return instance;
        }
 
        void discard(size_t n) { m_GRME.discard(n); }
 
        ScalableOperator sample() {
            auto hMat = m_GRME();
            // std::cout << hMat << std::endl;
            ScalableOperator res(std::move(hMat), g_dimLoc);
            return res;
        }
};
 
class ObservableEnsemble {
    public:
        using MatrixType    = Eigen::MatrixX<Scalar>;
        using EigenSolver   = Eigen::SelfAdjointEigenSolver<MatrixType>;
        using StateSpace    = ManyBodySpinSpace;
        using OperatorSpace = mBodyOperatorSpace<Scalar>;
 
        static size_t      class_M;
        static std::string dataname(size_t rep) {
            return "m" + std::to_string(class_M) + "_Observable" + "_No" + std::to_string(rep);
        };
        static std::string usage(size_t nargs) {
            std::string usage = "ObservableEnsemble Usage: ";
            usage += std::to_string(nargs + 0) + ".(m) ";
            return usage;
        };
        static constexpr size_t nargs(void) { return 1; };
 
    private:
        size_t                                       m_L;
        size_t                                       m_M;
        StateSpace                                   m_stateSpace;
        OperatorSpace                                m_opSpace;
        GaussianRandomMatrixGenerator<OperatorSpace> m_mbGRME;
 
    public:
        ObservableEnsemble(size_t L, size_t m) : m_L(L), m_M(m) {
            assert(m_L >= 1);
            class_M      = m_M;
            m_stateSpace = ManyBodySpinSpace(m_L, g_dimLoc);
            m_opSpace    = mBodyOperatorSpace<Scalar>(m_M, m_stateSpace);
            m_mbGRME     = GaussianRandomMatrixGenerator< decltype(m_opSpace) >(m_opSpace, m_M);
        }
 
        ObservableEnsemble(size_t L, char** CL_argv)
            : ObservableEnsemble(L, std::atoi(CL_argv[0])) {}
 
        ObservableEnsemble()                                     = default;
        ObservableEnsemble(const ObservableEnsemble&)            = delete;
        ObservableEnsemble& operator=(const ObservableEnsemble&) = delete;
        ObservableEnsemble(ObservableEnsemble&&)                 = default;
        ObservableEnsemble& operator=(ObservableEnsemble&&)      = default;
 
        StateSpace& stateSpace() { return m_stateSpace; }
 
        OperatorSpace& operatorSpace() { return m_opSpace; }
 
        void reset() { m_mbGRME.reset(); }
 
        void discard(size_t n) { m_mbGRME.discard(n); }
 
        MatrixType sample() {
            SubSpaceManager& subSpaces    = SubSpaceManager::get_instance();
            auto&            subSpaceBase = subSpaces[m_L].basis();
            auto             op           = m_mbGRME();
            // std::cout << "observable sample:\n" << op << std::endl;
            MatrixType res = subSpaceBase.adjoint() * op * subSpaceBase;
 
            // std::cout << "observable res:\n" << res << std::endl;
            return res;
        }
};
size_t ObservableEnsemble::class_M = 0;
#endif