Commit 8647901d authored by chbauman's avatar chbauman
Browse files

save to change branch :)

parent 58b74aff
......@@ -3,18 +3,18 @@ static std::normal_distribution<> dist_mvn;
struct normal_random_variable
{
normal_random_variable(Eigen::MatrixXd const& covar)
normal_random_variable(Eigen::Matrix<numeric_t, -1, -1> const& covar)
{
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigenSolver(covar);
Eigen::SelfAdjointEigenSolver<Eigen::Matrix<numeric_t, -1, -1>> eigenSolver(covar);
transform = eigenSolver.eigenvectors() * eigenSolver.eigenvalues().cwiseSqrt().asDiagonal();
}
Eigen::MatrixXd transform;
Eigen::Matrix<numeric_t, -1, -1> transform;
Eigen::VectorXd operator()() const
Eigen::Matrix<numeric_t, -1, 1> operator()() const
{
const int numberOfEntries = transform.cols();
Eigen::VectorXd randomStandardNormal(numberOfEntries);
Eigen::Matrix<numeric_t, -1, 1> randomStandardNormal(numberOfEntries);
for (int i = 0; i < numberOfEntries; ++i) {
randomStandardNormal(i) = dist_mvn(GEN);
}
......
......@@ -31,7 +31,7 @@ static inline Data parse_data(const char* filename) {
std::vector<numeric_t> raw_data;
std::size_t values_read = 0;
double new_value;
numeric_t new_value;
for (;;) {
data_stream >> new_value;
......@@ -75,7 +75,7 @@ static inline decltype(auto) get_likelihood_function_from_data(const char* filen
return [=] (const ThetaVector & theta) {
const Vector computed_potentials = data.col(0).unaryExpr(
[&] (const double r) {
[&] (const numeric_t r) {
return potential_function(r, theta);
});
const Vector measured_potentials = data.col(1);
......
......@@ -11,7 +11,7 @@
/// Compute the next \rho
static inline numeric_t next_rho(numeric_t prev_rho, const Vector & likelihoods){
/// f = (COV - 1)^2
auto f = [&](double rho_next){
auto f = [&](numeric_t rho_next){
assert(likelihoods.minCoeff() >= 0);
const Vector like_to_p = likelihoods.array().pow(rho_next - prev_rho);
const numeric_t sample_mean = like_to_p.sum() / POP_SIZE;
......@@ -27,8 +27,8 @@ static inline numeric_t next_rho(numeric_t prev_rho, const Vector & likelihoods)
#ifdef USE_BRENTS_METHOD
const double lower_bound_of_search_space = prev_rho;
const double upper_bound_of_search_space = 1;
const numeric_t lower_bound_of_search_space = prev_rho;
const numeric_t upper_bound_of_search_space = 1;
boost::uintmax_t maximum_number_of_iterations = 10000;
const int minimum_number_of_correct_binary_digits_in_mantissa = 9;
......@@ -58,7 +58,7 @@ static inline numeric_t next_rho(numeric_t prev_rho, const Vector & likelihoods)
#endif
std::cout << "Next rho: " << f_argmin << " with COV " << std::sqrt(f_min) + 1.0 << "\n";
const numeric_t new_rho = std::min(1.0, f_argmin);
const numeric_t new_rho = std::min(numeric_t(1.0), f_argmin);
assert(new_rho > prev_rho);
return new_rho;
......
......@@ -71,9 +71,28 @@ int main(){
weights /= weight_sum;
best_thetas.push_back(get_most_likely_theta(thetas, weights));
/// Resampling
resample(weights, resample_count);
/*
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);
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;
curr_number += temp;
++local_count;
}
}
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;
......
......@@ -41,7 +41,7 @@ struct PotentialName<Potential::BUCKINGHAM> {
};
constexpr index_t POP_SIZE = 400000;
constexpr index_t POP_SIZE = 6000;
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;
......@@ -53,10 +53,10 @@ typedef Eigen::Matrix<numeric_t, N_DIM, N_DIM> CovMat;
typedef Eigen::Matrix<numeric_t, N_DIM, -1> PopMatrix;
typedef double (*PotentialFunction)(const double r, const ThetaVector& theta);
typedef numeric_t (*PotentialFunction)(const numeric_t r, const ThetaVector& theta);
constexpr numeric_t NORMAL_PRIOR_STDDEV = 10;
constexpr numeric_t NORMAL_PRIOR_STDDEV = 0.01;
/// Beta
constexpr numeric_t BETA = 0.2;
......@@ -66,8 +66,8 @@ extern XoroshiroRandomNumberEngine GEN;
//#define UNIF_PRIOR
//#define USE_BRENTS_METHOD
#define USE_BRENTS_METHOD
constexpr numeric_t BOUND = 1.e2;
#pragma omp threadprivate(GEN)
//#define BURNIN 0
#define BURNIN 10
......@@ -34,7 +34,7 @@ inline static numeric_t mvn_iso(const Vector &x, const Vector &meanVec, const nu
/// 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<Eigen::MatrixXd> Chol;
typedef Eigen::LLT<Eigen::Matrix<numeric_t, -1, -1>> Chol;
Chol chol(covMat);
// Handle non positive definite covariance somehow:
#ifndef NDEBUG
......@@ -52,11 +52,13 @@ Eigen::Matrix<numeric_t, N_DIM - 1, 1> get_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;
}
}
......
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