Ensemble.hpp

Go to the documentation of this file.

#pragma once
#ifdef GPU
    #include "Ensemble.cuh"
#else
    #include "Headers/mytypes.hpp"
    #include "Headers/HilbertSpace/ParityTransSector.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_parity = 1;
 
using RealScalar = Real_t;
using Scalar     = Complex_t<RealScalar>;
 
class SubSpaceManager {
    public:
        using SubHilbertSpace = ParityTransSector<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]
                        = ParityTransSector<RealScalar>(g_parity, 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 const& 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) {
            std::stringstream buff("");
            buff << std::setprecision(3);
            buff << "Hamiltonian_B" << class_B << "_h" << class_h;
            return buff.str();
        };
        static std::string usage(size_t nargs) {
            std::string usage = "HamiltonianEnsemble Usage: ";
            usage += std::to_string(nargs + 0) + ".(Transverse field) ";
            usage += std::to_string(nargs + 1) + ".(Longitudinal field) ";
            return usage;
        };
        static constexpr size_t nargs(void) { return 2; };
 
        static StateSpace& stateSpace(size_t L) {
            SubSpaceManager& subSpaces = SubSpaceManager::get_instance();
            return subSpaces[L];
        }
 
    private:
        static double class_B;  // Transverse field
        static double class_h;  // Longitudinal field
        MatrixType    m_locMat;
 
        HamiltonianEnsemble(char** CL_argv) {
            class_B        = std::atof(CL_argv[0]);
            class_h        = std::atof(CL_argv[1]);
            size_t     dim = g_dimLoc * g_dimLoc;
            MatrixType sX  = MatrixType::Zero(g_dimLoc, g_dimLoc);
            MatrixType sY  = MatrixType::Zero(g_dimLoc, g_dimLoc);
            MatrixType sZ  = MatrixType::Zero(g_dimLoc, g_dimLoc);
 
            double s = (g_dimLoc - 1) / 2.0;
            for(auto j = 0; j != g_dimLoc - 1; ++j) {
                double m     = s - j;
                sZ(j, j)     = m;
                sX(j, j + 1) = std::sqrt((s - m + 1) * (s + m));
                sY(j + 1, j) = std::sqrt((s + m) * (s - m + 1));
            }
            sZ(g_dimLoc - 1, g_dimLoc - 1) = -s;
            sX                             = (sX + sY).eval() / 2.0;
            sY                             = Scalar(0, 1) * (sY - sX).eval();
 
            auto Identity = MatrixType::Identity(g_dimLoc, g_dimLoc);
 
            m_locMat = Eigen::kroneckerProduct(sZ, sZ);
            m_locMat
                += class_B
                   * (Eigen::kroneckerProduct(sX, Identity) + Eigen::kroneckerProduct(Identity, sX))
                   / 2.0;
            m_locMat
                += class_h
                   * (Eigen::kroneckerProduct(sZ, Identity) + Eigen::kroneckerProduct(Identity, sZ))
                   / 2.0;
            std::cout << m_locMat << std::endl;
        }
 
    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_locMat, g_dimLoc); }
};
double HamiltonianEnsemble::class_B = 0;
double HamiltonianEnsemble::class_h = 0;
 
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