Commit 6c913aef authored by chbauman's avatar chbauman
Browse files

conflicts resolved

parents 8647901d 898ade3d
......@@ -3,19 +3,19 @@ static std::normal_distribution<> dist_mvn;
struct normal_random_variable
{
normal_random_variable(Eigen::Matrix<numeric_t, -1, -1> const& covar)
normal_random_variable(CovMat const& covar)
{
Eigen::SelfAdjointEigenSolver<Eigen::Matrix<numeric_t, -1, -1>> eigenSolver(covar);
Eigen::SelfAdjointEigenSolver<CovMat> eigenSolver(covar);
transform = eigenSolver.eigenvectors() * eigenSolver.eigenvalues().cwiseSqrt().asDiagonal();
}
Eigen::Matrix<numeric_t, -1, -1> transform;
CovMat transform;
Eigen::Matrix<numeric_t, -1, 1> operator()() const
ThetaVector operator()() const
{
const int numberOfEntries = transform.cols();
Eigen::Matrix<numeric_t, -1, 1> randomStandardNormal(numberOfEntries);
for (int i = 0; i < numberOfEntries; ++i) {
ThetaVector randomStandardNormal;
for (int i = 0; i < N_DIM; ++i) {
randomStandardNormal(i) = dist_mvn(GEN);
}
return transform * randomStandardNormal;
......
......@@ -9,13 +9,13 @@
/// Compute the next \rho
static inline numeric_t next_rho(numeric_t prev_rho, const Vector & likelihoods){
static inline numeric_t next_rho(numeric_t prev_rho, const WeightVector & likelihoods){
/// f = (COV - 1)^2
auto f = [&](numeric_t rho_next){
assert(likelihoods.minCoeff() >= 0);
const Vector like_to_p = likelihoods.array().pow(rho_next - prev_rho);
const WeightVector like_to_p = likelihoods.array().pow(rho_next - prev_rho);
const numeric_t sample_mean = like_to_p.sum() / POP_SIZE;
const Vector like_to_p_min_mean = like_to_p.array() - sample_mean;
const WeightVector like_to_p_min_mean = like_to_p.array() - sample_mean;
const numeric_t sample_var = std::sqrt(like_to_p_min_mean.array().pow(2).sum() / (POP_SIZE - 1));
return std::pow(sample_var / sample_mean - 1, 2);
};
......@@ -24,7 +24,6 @@ static inline numeric_t next_rho(numeric_t prev_rho, const Vector & likelihoods)
numeric_t f_argmin;
numeric_t f_min;
#ifdef USE_BRENTS_METHOD
const numeric_t lower_bound_of_search_space = prev_rho;
......
......@@ -32,7 +32,7 @@ int main(){
std::cout << "Initializing ... ";
numeric_t rho_next, rho_curr = 0;
PopMatrix thetas(N_DIM, POP_SIZE);
thetas = prior_sample(POP_SIZE);
thetas = prior_sample();
PopMatrix thetas_new(N_DIM, POP_SIZE);
WeightVector weights(POP_SIZE);
WeightVector likelihoods(POP_SIZE);
......@@ -124,6 +124,7 @@ int main(){
}
}
assert(likelihoods.minCoeff() >= 0);
/// Prepare next iteration
std::cout << "Step " << step_count++ << " done.\n";
thetas.swap(thetas_new);
......
......@@ -41,21 +41,24 @@ struct PotentialName<Potential::BUCKINGHAM> {
};
constexpr index_t POP_SIZE = 6000;
constexpr index_t POP_SIZE = 400;
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;
constexpr int EIGEN_NUM_ROWS_FOR_POP = (POP_SIZE > MAX_POP_SIZE_FOR_STACK? Eigen::Dynamic: POP_SIZE);
constexpr index_t MAX_POP_SIZE_FOR_STACK = 5000;
constexpr index_t EIGEN_NUM_ROWS_FOR_POP = -1;//(POP_SIZE > MAX_POP_SIZE_FOR_STACK? Eigen::Dynamic: POP_SIZE);
typedef Eigen::Matrix<numeric_t, N_DIM, 1> ThetaVector;
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, -1> PopMatrix;
typedef Eigen::Matrix<numeric_t, N_DIM, EIGEN_NUM_ROWS_FOR_POP> PopMatrix;
typedef Eigen::Matrix<numeric_t, N_DIM, -1> PartPopMatrix;
typedef numeric_t (*PotentialFunction)(const numeric_t r, const ThetaVector& theta);
constexpr numeric_t NORMAL_PRIOR_STDDEV = 0.01;
/// Beta
......
......@@ -30,11 +30,12 @@ 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<Eigen::Matrix<numeric_t, -1, -1>> Chol;
typedef Eigen::LLT<Matrix> Chol;
Chol chol(covMat);
// Handle non positive definite covariance somehow:
#ifndef NDEBUG
......@@ -45,7 +46,7 @@ inline static numeric_t mvn(const Vector &x, const Vector &meanVec, const Matrix
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();
}
}*/
Eigen::Matrix<numeric_t, N_DIM - 1, 1> get_exact_solution() {
Eigen::Matrix<numeric_t, N_DIM - 1, 1> exact_solution;
......@@ -73,11 +74,11 @@ static std::normal_distribution<numeric_t> normal_prior(1.0, NORMAL_PRIOR_STDDEV
static std::gamma_distribution<numeric_t> var_dist(var_mean * var_mean, var_mean);
static std::uniform_real_distribution<numeric_t> unif_dist(0.,1.);
inline static PopMatrix prior_sample(const index_t n_samples){
PopMatrix prior(N_DIM, n_samples);
inline static PopMatrix prior_sample(){
PopMatrix prior(N_DIM, POP_SIZE);
#ifdef UNIF_PRIOR
#pragma omp parallel for schedule(static)
for(index_t j = 0; j < n_samples; ++j){
for(index_t j = 0; j < POP_SIZE; ++j){
for(index_t i = 0; i < N_DIM - 1; ++i){
prior(i,j) = unif_dist(GEN) * BOUND;
}
......@@ -86,7 +87,7 @@ inline static PopMatrix prior_sample(const index_t n_samples){
#else
const auto exact_solution = get_exact_solution();
#pragma omp parallel for schedule(static)
for(index_t i = 0; i < n_samples; ++i){
for(index_t i = 0; i < POP_SIZE; ++i){
for (index_t j = 0; j < N_DIM-1; ++j) {
prior(j, i) = normal_prior(GEN) * exact_solution(j);
}
......@@ -125,10 +126,10 @@ inline static void resample(const WeightVector & weights, CountVector & resample
/// Do 'rck' steps of MCMC rck times
template <typename Likelihood>
inline static PopMatrix mcmc(const normal_random_variable & sample, const ThetaVector & theta, const index_t rck, const ThetaVector & sample_mean,
inline static PartPopMatrix mcmc(const normal_random_variable & sample, const ThetaVector & theta, const index_t rck, const ThetaVector & sample_mean,
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);
PartPopMatrix new_thetas(N_DIM, rck);
ThetaVector particle, particle_prop;
std::uniform_real_distribution<numeric_t> acc_dist(0,1);
particle = theta;
......@@ -181,8 +182,6 @@ inline static PopMatrix mcmc(const normal_random_variable & sample, const ThetaV
};
Eigen::Matrix<numeric_t, N_DIM - 1, 1> relative_error_to_exact_solution(const ThetaVector& theta) {
return (theta.head(N_DIM-1) - get_exact_solution()).array().abs().matrix().cwiseQuotient(get_exact_solution().array().abs().matrix());
}
......
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