RandomMatrix.hpp

Go to the documentation of this file.

#pragma once
 
#include "Headers/HilbertSpace/OperatorSpace.hpp"
#include <iostream>
#include <random>
 
template<class Derived, typename std::enable_if_t<
                            std::is_convertible_v<Derived, OperatorSpaceBase<Derived>> >* = nullptr>
class GaussianRandomMatrixGenerator {
    private:
        using Scalar     = typename Derived::Scalar;
        using RealScalar = typename Eigen::NumTraits<typename Derived::Scalar>::Real;
        Derived                          m_opSpace;
        int                              m_seed;
        std::mt19937                     mt;
        std::normal_distribution<double> Gaussian;
        mutable unsigned long long       count;
 
    public:
        GaussianRandomMatrixGenerator(Derived const& opSpace, int seed = 0, RealScalar stddev = 1.0)
            : m_opSpace{opSpace}, m_seed(seed), mt(seed), Gaussian{0.0, stddev}, count{0} {}
        GaussianRandomMatrixGenerator() = default;
        GaussianRandomMatrixGenerator(GaussianRandomMatrixGenerator const& other) = delete;
        GaussianRandomMatrixGenerator(GaussianRandomMatrixGenerator&& other) = default;
        ~GaussianRandomMatrixGenerator() = default;
 
        /* @{ */
        __host__ __device__ GaussianRandomMatrixGenerator& operator=(
            GaussianRandomMatrixGenerator const& other)
            = delete;
        __host__ __device__ GaussianRandomMatrixGenerator& operator=(
            GaussianRandomMatrixGenerator&& other)
            = default;
        /* @} */
 
        void reset() { mt = std::mt19937(m_seed); }
 
        void discard(unsigned long long num) {
            for(auto j = 0; j < num * m_opSpace.dim(); ++j) Gaussian(mt);
            count += num * m_opSpace.dim();
        }
 
        Eigen::MatrixX<Scalar> operator()() {
            Eigen::VectorXd coeff;
            coeff = coeff.NullaryExpr(m_opSpace.dim(), [&]() { return Gaussian(mt); });
            // std::cout << coeff << std::endl;
 
            Eigen::MatrixX<Scalar> sample
                = Eigen::MatrixX<Scalar>::Zero(m_opSpace.baseDim(), m_opSpace.baseDim());
// #pragma omp declare reduction(+ : Eigen::MatrixX<Scalar> : omp_out=omp_out+omp_in) initializer(omp_priv = omp_orig)
// #pragma omp parallel for reduction(+ : sample)
            for(auto j = 0; j < m_opSpace.dim(); ++j) { sample += coeff(j) * m_opSpace.basisOp(j); }
 
            count += m_opSpace.dim();
            return sample;
        }
};
 
#ifdef GPU
    #include "RandomMatrix.cuh"
#endif