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>
 
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;
            using RandomMatrixGenerators
                = std::vector< GaussianRandomMatrixGenerator<OperatorSpace<Scalar>> >;
            RandomMatrixGenerators& m_GRME;
            std::vector<bool>&      m_isConstructed;
 
        public:
            ScalableOperator(RandomMatrixGenerators& GRME, std::vector<bool>& isConstructed)
                : m_GRME{GRME}, m_isConstructed{isConstructed} {}
 
            MatrixType construct_globalOp(size_t L) {
                auto currentSize = m_GRME.size();
                if(m_isConstructed.size() < L || !m_isConstructed[L]) {
                    m_GRME.resize(L + 1);
                    m_isConstructed.resize(L + 1, false);
                    SubSpaceManager& subSpaces = SubSpaceManager::get_instance();
                    auto             dim       = subSpaces[L].dim();
                    auto             HSpace    = HilbertSpace<int>(dim);
                    auto             OpSpace   = OperatorSpace<Scalar>(HSpace);
                    m_GRME[L] = GaussianRandomMatrixGenerator< decltype(OpSpace) >(OpSpace, 0);
                    m_isConstructed[L] = true;
                }
                return m_GRME[L]();
            };
    };
}  // 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: ";
            return usage;
        };
        static constexpr size_t nargs(void) { return 0; };
 
        static StateSpace& stateSpace(size_t L) {
            SubSpaceManager& subSpaces = SubSpaceManager::get_instance();
            return subSpaces[L];
        }
 
    private:
        std::vector< GaussianRandomMatrixGenerator<OperatorSpace<Scalar>> > m_GRME;
        std::vector<bool>                                                   m_isConstructed;
 
        HamiltonianEnsemble(char** CL_argv) {}
 
    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) {}
 
        ScalableOperator sample() { return ScalableOperator(m_GRME, m_isConstructed); }
};
 
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