MultiVariateGaussian.h

Go to the documentation of this file.

#ifndef IPPL_MULTI_VARIATE_GAUSSIAN_H
#define IPPL_MULTI_VARIATE_GAUSSIAN_H
 
#include <Kokkos_Random.hpp>
#include <cmath>
#include <memory>
#include "Distribution.h"
#include "Ippl.h"
#include "OPALTypes.h"
#include "SamplingBase.hpp"
#include "Utilities/Options.h"
 
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,
            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;
 
    void emitParticles(double t, double dt) 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