Commit 6486dcf6 authored by chbauman's avatar chbauman
Browse files

small changes

parent 1495fbeb
static std::normal_distribution<> dist_mvn;
struct normal_random_variable
{
/// Struct that returns a sample of a 'N_DIM'-dimensional Gaussian with mean zero and covariance matrix
/// 'transform'
struct normal_random_variable {
normal_random_variable(CovMat const& covar)
{
normal_random_variable(CovMat const& covar){
Eigen::SelfAdjointEigenSolver<CovMat> eigenSolver(covar);
transform = eigenSolver.eigenvectors() * eigenSolver.eigenvalues().cwiseSqrt().asDiagonal();
}
CovMat transform;
ThetaVector operator()() const
{
ThetaVector operator()() const {
ThetaVector randomStandardNormal;
for (int i = 0; i < N_DIM; ++i) {
for(index_t i = 0; i < N_DIM; ++i) {
randomStandardNormal(i) = dist_mvn(GEN);
}
return transform * randomStandardNormal;
......
......@@ -10,7 +10,6 @@
#include <cstddef>
static inline numeric_t lennard_jones_potential_function(const numeric_t r, const ThetaVector& theta) {
const numeric_t sigma_div_r_to_the_sixth = std::pow(theta(1) / r, 6);
return 4 * theta(0) * (sigma_div_r_to_the_sixth * sigma_div_r_to_the_sixth - sigma_div_r_to_the_sixth);
}
......@@ -19,7 +18,6 @@ static inline numeric_t buckingham_potential_function(const numeric_t r, const T
return std::exp(theta(0) - theta(1) * r) - std::exp(theta(2)) / std::pow(r, 6);
}
static inline Data parse_data(const char* filename) {
std::ifstream data_stream(filename);
......
......@@ -122,6 +122,21 @@ int main(){
}
#else
index_t ct = 0;
#pragma omp parallel shared(sample_mean, sample_cov_mat, rho_curr, thetas)
{
#pragma omp single
for(index_t k = 0; k < POP_SIZE; ++k){
const index_t rck = resample_count(k);
if(rck > 0){
#pragma omp task firstprivate(ct, rck)
{
thetas_new.middleCols(ct, rck) = mcmc(sample, thetas.col(k), rck,
sample_mean, sample_cov_mat, rho_curr, likelihoods, likelihoods_old, likelihood, k, ct);
}
ct += rck;
}
}
}
#endif
assert(likelihoods.minCoeff() >= 0);
......
......@@ -10,7 +10,6 @@ typedef Eigen::Matrix<numeric_t, Eigen::Dynamic, 2> Data;
enum class Potential {LENNARD_JONES, BUCKINGHAM};
#if defined(USE_LENNARD_JONES_POTENTIAL) && defined(USE_BUCKINGHAM_POTENTIAL)
#error "Can't specify both potentials!"
#elif !defined(USE_LENNARD_JONES_POTENTIAL) && !defined(USE_BUCKINGHAM_POTENTIAL)
......@@ -18,7 +17,6 @@ enum class Potential {LENNARD_JONES, BUCKINGHAM};
#define USE_LENNARD_JONES_POTENTIAL
#endif
#if defined(USE_LENNARD_JONES_POTENTIAL)
constexpr Potential POTENTIAL = Potential::LENNARD_JONES;
#elif defined(USE_BUCKINGHAM_POTENTIAL)
......@@ -40,10 +38,8 @@ struct PotentialName<Potential::BUCKINGHAM> {
static constexpr const char* value = "Buckingham";
};
constexpr index_t POP_SIZE = 100000;
constexpr index_t N_DIM = (POTENTIAL == Potential::LENNARD_JONES? 2: 3) + 1;// + 1 for the unknown standard deviation
constexpr index_t MAX_POP_SIZE_FOR_STACK = 5000;
constexpr index_t EIGEN_NUM_ROWS_FOR_POP = -1;//(POP_SIZE > MAX_POP_SIZE_FOR_STACK? Eigen::Dynamic: POP_SIZE);
......@@ -57,9 +53,11 @@ typedef Eigen::Matrix<numeric_t, N_DIM, -1> PartPopMatrix;
typedef numeric_t (*PotentialFunction)(const numeric_t r, const ThetaVector& theta);
/// 'Relative' Standard Deviation for Gaussian Prior
constexpr numeric_t NORMAL_PRIOR_STDDEV = 0.5;
constexpr numeric_t NORMAL_PRIOR_STDDEV = 0.01;
/// Bound for uniform Prior, [0, BOUND] in each dimension
constexpr numeric_t BOUND = 1.e2;
/// Beta
constexpr numeric_t BETA = 0.2;
......@@ -67,10 +65,10 @@ constexpr numeric_t BETA = 0.2;
/// RNG
extern XoroshiroRandomNumberEngine GEN;
/// Defines
//#define UNIF_PRIOR
#define USE_BRENTS_METHOD
constexpr numeric_t BOUND = 1.e2;
#define USE_FOR_VERSION
#define BURNIN 10
#pragma omp threadprivate(GEN)
#define BURNIN 10
......@@ -5,7 +5,7 @@
#include <Eigen/Eigenvalues>
#include <cmath>
/// PDF of a Gaussian distribution with mean 'm' and std 's'
/// PDF of a 1D-Gaussian distribution with mean 'm' and std 's'
template <typename T>
T normal_pdf_1D(T x, T m, T s){
static const T inv_sqrt_2pi = 0.3989422804014327;
......@@ -20,7 +20,7 @@ template <typename T>
T gamma_pdf_1D(T x, T a, T b){
return std::exp(-x / b) * std::pow(x, a - 1.) / std::pow(b, a) / tgamma(a);
}
/// PDF of an isotropic Gaussian distribution with std 'st_dev' in each dimension
/// PDF of an isotropic multi-variate Gaussian distribution with std 'st_dev' in each dimension
inline static numeric_t mvn_iso(const Vector &x, const Vector &meanVec, const numeric_t st_dev) {
return std::exp(x.binaryExpr(
meanVec,
......@@ -30,24 +30,7 @@ inline static numeric_t mvn_iso(const Vector &x, const Vector &meanVec, const nu
).sum());
}
/*
/// From StackOverflow
/// https://stackoverflow.com/questions/41538095/evaluate-multivariate-normal-gaussian-density-in-c
inline static numeric_t mvn(const Vector &x, const Vector &meanVec, const Matrix &covMat){
const numeric_t logSqrt2Pi = 0.5*std::log(2*M_PI);
typedef Eigen::LLT<Matrix> Chol;
Chol chol(covMat);
// Handle non positive definite covariance somehow:
#ifndef NDEBUG
//std::cerr << "Eigen values " << covMat.eigenvalues() << "\n";
#endif
if(chol.info()!=Eigen::Success) throw std::runtime_error("decomposition failed!");
const Chol::Traits::MatrixL& L = chol.matrixL();
numeric_t quadform = (L.solve(x - meanVec)).squaredNorm();
return std::exp(-x.rows()*logSqrt2Pi - 0.5*quadform) / L.determinant();
}*/
/// Returns the "exact" parameters depending on the potential
Eigen::Matrix<numeric_t, N_DIM - 1, 1> get_exact_solution() {
Eigen::Matrix<numeric_t, N_DIM - 1, 1> exact_solution;
switch (POTENTIAL) {
......@@ -63,12 +46,10 @@ Eigen::Matrix<numeric_t, N_DIM - 1, 1> get_exact_solution() {
break;
}
}
return exact_solution;
}
/// Returns sample from prior f(\Theta)
/// Returns 'POP_SIZE' samples from the prior
constexpr numeric_t var_mean = (POTENTIAL == Potential::LENNARD_JONES? 4.823: 1.26075);
static std::normal_distribution<numeric_t> normal_prior(1.0, NORMAL_PRIOR_STDDEV); // 0.2 * 0.2
static std::gamma_distribution<numeric_t> var_dist(var_mean * var_mean, var_mean);
......@@ -97,16 +78,14 @@ inline static PopMatrix prior_sample(){
return prior;
}
/// Returns the value of the prior pdf at theta
/// Returns the value of the prior pdf at theta if the prior is Gaussian
inline static numeric_t prior(const ThetaVector & theta){
numeric_t prior_pdf = 1.0;
const auto exact_solution = get_exact_solution();
for (index_t i = 0; i < N_DIM - 1; ++i) {
prior_pdf *= normal_pdf_1D<numeric_t>(theta(i), exact_solution(i), exact_solution(i) * normal_prior.stddev());
}
prior_pdf *= gamma_pdf_1D<numeric_t>(theta(N_DIM - 1), var_dist.alpha(), var_dist.beta());
return prior_pdf;
}
......@@ -165,7 +144,6 @@ inline static PartPopMatrix mcmc(const normal_random_variable & sample, const Th
#else
new_thetas.col(i) = particle;
new_likelihoods(ct + i) = particle_likelihood;
//std::cout << particle_likelihood << "<-part. lh\n";
#endif
/*
std::cout << "Cov mat " << sample_cov_mat << "\n"
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment