Commit f407a40c authored by stefanow's avatar stefanow
Browse files

merge

parents cf741a1d 599d4c39
......@@ -36,6 +36,7 @@ int main(){
PopMatrix thetas_new(N_DIM, POP_SIZE);
WeightVector weights(POP_SIZE);
WeightVector likelihoods(POP_SIZE);
WeightVector likelihoods_old(POP_SIZE);
CountVector resample_count(POP_SIZE);
std::vector<WeightVector> best_thetas;
......@@ -43,17 +44,19 @@ int main(){
ThetaVector sample_mean;
numeric_t evidence_s = 1.0;
const auto likelihood = get_likelihood_function_from_data("data.txt");
/// Compute likelihoods
#pragma omp parallel for schedule(static)
for(index_t i = 0; i < POP_SIZE; ++i){
likelihoods_old(i) = likelihood(thetas.col(i));
std::cout << likelihoods_old(i) << "<-likelihoods_old(i) \n";
}
std::cout << "Done.\n";
/// Iterate
index_t step_count = 0;
do{
/// Compute likelihoods
#pragma omp parallel for schedule(static)
for(index_t i = 0; i < POP_SIZE; ++i){
likelihoods(i) = likelihood(thetas.col(i));
}
rho_next = next_rho(rho_curr, likelihoods);
rho_next = next_rho(rho_curr, likelihoods_old);
/// Compute weights
#pragma omp parallel for schedule(static)
......@@ -94,16 +97,20 @@ int main(){
if(rck > 0){
#pragma omp task firstprivate(ct, rck)
{
thetas_new.middleCols(ct, rck) = mcmc(thetas.col(k), rck, sample_mean, sample_cov_mat, rho_curr, likelihoods(k), likelihood);
thetas_new.middleCols(ct, rck) = mcmc(thetas.col(k), rck,
sample_mean, sample_cov_mat, rho_curr, likelihoods, likelihoods_old, likelihood, k, ct);
}
ct += rck;
}
}
}
assert(likelihoods.minCoeff() >= 0);
std::cout << likelihoods.minCoeff() << "<- min likelihood\n";
std::cout << likelihoods.maxCoeff() << "<- max likelihood\n";
/// Prepare next iteration
std::cout << "Step " << step_count++ << " done.\n";
thetas.swap(thetas_new);
likelihoods.swap(likelihoods_old);
rho_curr = rho_next;
} while(rho_curr < 1.0);
......
......@@ -67,7 +67,7 @@ extern XoroshiroRandomNumberEngine GEN;
//#define UNIF_PRIOR
//#define USE_BRENTS_METHOD
constexpr numeric_t BOUND = 1e2;
constexpr numeric_t BOUND = 1.e2;
#pragma omp threadprivate(GEN)
#define BURNIN 2
//#define BURNIN 0
......@@ -63,7 +63,7 @@ Eigen::Matrix<numeric_t, N_DIM - 1, 1> get_exact_solution() {
return exact_solution;
}
/// Returns sample from prior f(\Theta)
constexpr numeric_t var_mean = (POTENTIAL == Potential::LENNARD_JONES? 4.823: 1.26075);
......@@ -128,13 +128,14 @@ inline static void resample(const WeightVector & weights, CountVector & resample
/// Do 'rck' steps of MCMC rck times
template <typename Likelihood>
inline static PopMatrix mcmc(const ThetaVector & theta, const index_t rck, const ThetaVector & sample_mean,
const CovMat & sample_cov_mat, const numeric_t rho_curr, const numeric_t likelihood_theta,
const Likelihood & likelihood){
const CovMat & sample_cov_mat, const numeric_t rho_curr, WeightVector & new_likelihoods,
const WeightVector & likelihoods, const Likelihood & likelihood, const index_t k, const index_t ct){
PopMatrix new_thetas(N_DIM, rck);
ThetaVector particle, particle_prop;
std::uniform_real_distribution<numeric_t> acc_dist(0,1);
particle = theta;
numeric_t particle_likelihood = likelihood_theta;
numeric_t particle_likelihood = likelihoods(k);
std::cout << particle_likelihood << "<-part. lh old\n";
#ifdef BURNIN
for(index_t i = 0; i < rck + BURNIN; ++i){
#else
......@@ -159,9 +160,14 @@ inline static PopMatrix mcmc(const ThetaVector & theta, const index_t rck, const
}
#endif
#ifdef BURNIN
if(i >= BURNIN) new_thetas.col(i - BURNIN) = particle;
if(i >= BURNIN) {
new_thetas.col(i - BURNIN) = particle;
new_likelihoods(ct + i - BURNIN) = particle_likelihood;
}
#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