MultiVariateGaussian.h

Go to the documentation of this file.

#ifndef IPPL_MULTI_VARIATE_GAUSSIAN_H
#define IPPL_MULTI_VARIATE_GAUSSIAN_H
 
#include "Distribution.h"
#include "SamplingBase.hpp"
#include <Kokkos_Random.hpp>
#include "Ippl.h"
#include "Utilities/Options.h"
#include "OPALTypes.h"
#include <memory>
#include <cmath>
 
using ParticleContainer_t = ParticleContainer<double, 3>;
using FieldContainer_t = FieldContainer<double, 3>;
using Distribution_t = Distribution;
using GeneratorPool = typename Kokkos::Random_XorShift64_Pool<>;
using Dist_t = ippl::random::NormalDistribution<double, 3>;
using Matrix_t = ippl::Vector<ippl::Vector<double, 6>, 6>;

class MultiVariateGaussian : public SamplingBase {
public:
    MultiVariateGaussian(std::shared_ptr<ParticleContainer_t> &pc, 
                         std::shared_ptr<FieldContainer_t> &fc, 
                         std::shared_ptr<Distribution_t> &opalDist);

    MultiVariateGaussian(std::shared_ptr<ParticleContainer_t> pc,
                   const Vector_t<double, 3>& meanR,
                   const Vector_t<double, 3>& meanP,
                   const Vector_t<double, 3>& sigmaR,
                   const Vector_t<double, 3>& sigmaP,
                   const Vector_t<double, 3>& cutoffR,
                   const Vector_t<double, 3>& cutoffP,
                   bool fixMeanR=true,
                   bool fixMeanP=true);

    MultiVariateGaussian(std::shared_ptr<ParticleContainer_t> pc,
                   const Vector_t<double, 3>& meanR,
                   const Vector_t<double, 3>& meanP,
                   const Matrix_t &cov,
                   const Vector_t<double, 3>& cutoffR,
                   const Vector_t<double, 3>& cutoffP,
                   bool fixMeanR=true,
                   bool fixMeanP=true);
    void ComputeCholeskyFactorization();

    void ComputeCenteredBounds();

    void generateParticles(size_t &numberOfParticles, Vector_t<double, 3> nr) override;

    IpplTimings::TimerRef samplerTimer_m;
 
    void setMeanR(const Vector_t<double, 3>& meanR) {
        meanR_m = meanR;
    }
 
    void setMeanP(const Vector_t<double, 3>& meanP) {
        meanP_m = meanP;
    }
 
    void setCutoffR(const Vector_t<double, 3>& cutoffR) {
        cutoffR_m = cutoffR;
    }
 
    void setCutoffP(const Vector_t<double, 3>& cutoffP) {
        cutoffP_m = cutoffP;
    }
 
    void setFixMeanR(bool fixMeanR) {
        fixMeanR_m = fixMeanR;
    }
 
    void setFixMeanP(bool fixMeanP) {
        fixMeanP_m = fixMeanP;
    }
 
    void setSigmaR(const Vector_t<double, 3>& sigmaR) {
        sigmaR_m = sigmaR;
    }
    void setSigmaP(const Vector_t<double, 3>& sigmaP) {
        sigmaP_m = sigmaP;
    }
 
    void setCovarianceMatrix(const Matrix_t &cov) {
        cov_m = cov;
    }
 
    void setL(const Matrix_t &L) {
        L_m = L;
    }
 
private:
    /*
    * @brief Mean vectors for position and momentum.
    */
    Vector_t<double, 3> meanR_m, meanP_m;
 
    /*
    * @brief Covariance matrix for the Gaussian distribution.
    */
    Matrix_t cov_m;
 
    /*
    * @brief Cholesky cov=LT*L factorized matrix for the Gaussian distribution.
    */
    Matrix_t L_m;
 
    /*
    * @brief Sampling bounds for the particle distribution.
    */
    Vector_t<double, 3> rmin_m, rmax_m, pmin_m, pmax_m;
 
    /*
    * @brief Normalized sampling bounds for the particle distribution.
    */
    Vector_t<double, 3> normRmin_m, normRmax_m, normPmin_m, normPmax_m;

    Vector_t<double, 6> min_m, max_m, normMin_m, normMax_m;

    Vector_t<double, 3> cutoffR_m;
    Vector_t<double, 3> cutoffP_m;

    void initRandomPool();

    GeneratorPool randPool_m;

    Vector_t<double, 3> sigmaR_m;
    Vector_t<double, 3> sigmaP_m;

    bool fixMeanR_m;
    bool fixMeanP_m;
};
 
#endif // IPPL_MULTI_VARIATE_GAUSSIAN_H