Commit 599d4c39 authored by chbauman's avatar chbauman
Browse files

less likelihood

parent 4fc3734f
......@@ -32,6 +32,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;
......@@ -39,17 +40,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)
......@@ -88,16 +91,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;
best_thetas.push_back(get_most_likely_theta(thetas, weights));
} while(rho_curr < 1.0);
......
......@@ -41,7 +41,7 @@ struct PotentialName<Potential::BUCKINGHAM> {
};
constexpr index_t POP_SIZE = 100000;
constexpr index_t POP_SIZE = 400000;
constexpr index_t N_DIM = (POTENTIAL == Potential::LENNARD_JONES? 2: 3) + 1;// + 1 for the unknown standard deviation
constexpr int MAX_POP_SIZE_FOR_STACK = 5000;
......@@ -67,7 +67,7 @@ extern XoroshiroRandomNumberEngine GEN;
//#define UNIF_PRIOR
//#define USE_BRENTS_METHOD
constexpr numeric_t BOUND = 1.5e2;
constexpr numeric_t BOUND = 1.e2;
#pragma omp threadprivate(GEN)
#define BURNIN 5
//#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