Commit 1495fbeb authored by chbauman's avatar chbauman
Browse files

saved

parent 6c913aef
......@@ -38,14 +38,13 @@ int main(){
WeightVector likelihoods(POP_SIZE);
WeightVector likelihoods_old(POP_SIZE);
CountVector resample_count(POP_SIZE);
std::vector<WeightVector> best_thetas;
CovMat sample_cov_mat;
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));
......@@ -55,6 +54,7 @@ int main(){
/// Iterate
index_t step_count = 0;
do{
/// Compute the next rho
rho_next = next_rho(rho_curr, likelihoods_old);
/// Compute weights
......@@ -72,13 +72,15 @@ int main(){
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, 2> indices_for_mcmc(nnonzero, 2);
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){
......@@ -86,13 +88,15 @@ int main(){
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;
}
}
std::cout << resample_count.transpose() << "<- res count\n";
std::cout << indices_for_mcmc.transpose() << "<- indices\n";
*/
#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;
......@@ -105,24 +109,20 @@ 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};
#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,
#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);
}
ct += rck;
}
}
}
#else
index_t ct = 0;
#endif
assert(likelihoods.minCoeff() >= 0);
/// Prepare next iteration
......
......@@ -41,7 +41,7 @@ struct PotentialName<Potential::BUCKINGHAM> {
};
constexpr index_t POP_SIZE = 400;
constexpr index_t POP_SIZE = 100000;
constexpr index_t N_DIM = (POTENTIAL == Potential::LENNARD_JONES? 2: 3) + 1;// + 1 for the unknown standard deviation
......
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