Commit f2e650e5 authored by stefanow's avatar stefanow
Browse files

Other pop size

parents 5b7b587e 2ed08c00
......@@ -2,21 +2,20 @@
#include "typedefs_and_globals.hpp"
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;
......
......@@ -12,7 +12,6 @@
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);
}
......@@ -21,7 +20,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);
......
......@@ -46,8 +46,8 @@ int main(){
ThetaVector sample_mean;
numeric_t evidence_s = 1.0;
const auto likelihood = get_likelihood_function_from_data("data.txt");
/// Compute likelihoods
/// Compute likelihoods of prior samples
#pragma omp parallel for schedule(static)
for(index_t i = 0; i < POP_SIZE; ++i){
likelihoods_old(i) = likelihood(thetas.col(i));
......@@ -57,6 +57,7 @@ int main(){
/// Iterate
index_t step_count = 0;
do{
/// Compute the next rho
rho_next = next_rho(rho_curr, likelihoods_old);
/// Compute weights
......@@ -73,9 +74,32 @@ int main(){
weights /= weight_sum;
best_thetas.push_back(get_most_likely_theta(thetas, weights));
/// Resampling
resample(weights, resample_count);
resample(weights, resample_count);
#ifdef USE_FOR_VERSION
/// Precomputing indices for mcmc chains
index_t nnonzero = 0;
for(index_t i = 0; i < POP_SIZE; ++i){
if(resample_count(i) > 0) ++nnonzero;
}
Eigen::Matrix<index_t, -1, 3> indices_for_mcmc(nnonzero, 3);
index_t local_count = 0;
index_t curr_number = 0;
for(index_t i = 0; i < POP_SIZE; ++i){
if(resample_count(i) > 0) {
index_t temp = resample_count(i);
indices_for_mcmc(local_count, 0) = temp;
indices_for_mcmc(local_count, 1) = curr_number;
indices_for_mcmc(local_count, 2) = i;
curr_number += temp;
++local_count;
}
}
#endif
//std::cout << resample_count.transpose() << "<- res count\n";
//std::cout << indices_for_mcmc.transpose() << "<- indices\n";
/// Compute sample mean and covariance
sample_mean = thetas * weights;
......@@ -88,9 +112,19 @@ int main(){
sample_cov_mat += 0.001 * CovMat::Identity();
sample_cov_mat *= BETA * BETA;
/// MCMC
index_t ct = 0;
/// MCMC
normal_random_variable sample {sample_cov_mat};
#ifdef USE_FOR_VERSION
#pragma omp parallel for schedule(static)
for(index_t i = 0; i < nnonzero; ++i){
const index_t k = indices_for_mcmc(i, 2);
const index_t ct = indices_for_mcmc(i, 1);
const index_t rck = indices_for_mcmc(i, 0);
thetas_new.middleCols(ct, rck) = mcmc(sample, thetas.col(k), rck,
sample_mean, sample_cov_mat, rho_curr, likelihoods, likelihoods_old, likelihood, k, ct);
}
#else
index_t ct = 0;
#pragma omp parallel shared(sample_mean, sample_cov_mat, rho_curr, thetas)
{
#pragma omp single
......@@ -106,6 +140,7 @@ int main(){
}
}
}
#endif
assert(likelihoods.minCoeff() >= 0);
/// Prepare next iteration
......
......@@ -17,7 +17,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)
......@@ -25,7 +24,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)
......@@ -58,13 +56,16 @@ typedef Eigen::Matrix<numeric_t, EIGEN_NUM_ROWS_FOR_POP, 1> WeightVector;
typedef Eigen::Matrix<index_t, EIGEN_NUM_ROWS_FOR_POP, 1> CountVector;
typedef Eigen::Matrix<numeric_t, N_DIM, N_DIM> CovMat;
typedef Eigen::Matrix<numeric_t, N_DIM, EIGEN_NUM_ROWS_FOR_POP> PopMatrix;
typedef Eigen::Matrix<numeric_t, N_DIM, -1> PartPopMatrix;
typedef Eigen::Matrix<numeric_t, N_DIM, EIGEN_NUM_ROWS_FOR_POP> 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 = 1;
/// Bound for uniform Prior, [0, BOUND] in each dimension
constexpr numeric_t BOUND = 1.e2;
/// Beta
constexpr numeric_t BETA = 0.2;
......@@ -72,10 +73,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 0
......@@ -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;
......@@ -24,7 +24,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,
......@@ -34,43 +34,26 @@ 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) {
case Potential::LENNARD_JONES: {
exact_solution << 93.2, 3.372;
//exact_solution << 93.0312, 3.37062;
break;
}
case Potential::BUCKINGHAM: {
exact_solution << std::log(8.508e7), 1 / 0.273, std::log(5.13979e5);
//exact_solution << std::log(1.45709e8), 3.86461, std::log(468683);
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);
......@@ -99,16 +82,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;
}
......@@ -167,7 +148,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