Commit 63b96028 authored by stefanow's avatar stefanow
Browse files

log_normal

parent 1b2a3918
......@@ -8,7 +8,7 @@ optimizationFlags := -O3 -march=native
macroDefines := -DNDEBUG
else
$(info Debug build)
optimizationFlags := -O0 -g
optimizationFlags := -O0
endif
ifdef UNIF_PRIOR
......@@ -18,6 +18,13 @@ else
$(info Gaussian prior)
endif
ifdef USE_FLOATS
$(info using floats)
macroDefines += -DUSE_FLOATS
else
$(info using numeric_ts)
endif
ifdef USE_BRENTS_METHOD
$(info Using Brent's method to find the tho values, needs boost)
......@@ -41,7 +48,7 @@ boostIncludes :=
endif
CFLAGS := -fopenmp -std=c++14 $(optimizationFlags) $(macroDefines)
CFLAGS := -g -fopenmp -std=c++14 $(optimizationFlags) $(macroDefines) $(CFLAGS)
eigenIncludes := $(shell pkg-config --cflags eigen3)
......
CXX ?= g++
macroDefines :=
ifdef NDEBUG
$(info Release build)
optimizationFlags := -O3 -march=native
macroDefines := -DNDEBUG
else
$(info Debug build)
optimizationFlags := -O0 -g
endif
ifdef UNIF_PRIOR
$(info Uniform prior)
macroDefines += -DUNIF_PRIOR
else
$(info Gaussian prior)
endif
ifdef USE_FLOATS
$(info using floats)
macroDefines += -DUSE_FLOATS
else
$(info using numeric_ts)
endif
ifdef USE_BRENTS_METHOD
$(info Using Brent's method to find the tho values, needs boost)
macroDefines += -DUSE_BRENTS_METHOD
hostname = $(shell hostname)
ifeq ($(hostname),swmb)
$(info Compiling on Stefano's laptop)
boostIncludes := -I/opt/local/include
endif
ifneq ($(findstring euler,$(hostname)),)
$(info Compiling on EULER, wish me luck!)
$(shell module load boost/1.57.0)
boostIncludes := -I/cluster/apps/boost/1.57.0/x86_64/gcc_4.8.2/serial/include
endif
else
$(info Not using Brent's method)
boostIncludes :=
endif
CFLAGS := -fopenmp -std=c++14 $(optimizationFlags) $(macroDefines)
eigenIncludes := $(shell pkg-config --cflags eigen3)
ifeq ($(eigenIncludes),)
$(shell module load eigen)
eigenIncludes := -I/cluster/apps/eigen/3.2.1/x86_64/gcc_4.8.2/serial/include
endif
sources := $(wildcard *.cpp)
headers := $(wildcard *.hpp)
executable := execLJ execBuck execLJ.dSYM execBuck.dSYM
toClean := $(executable)
.PHONY: all
all: run
execLJ: $(sources) $(headers)
$(CXX) $(CFLAGS) $(eigenIncludes) $(boostIncludes) -DUSE_LENNARD_JONES_POTENTIAL -o execLJ $(sources)
execBuck: $(sources) $(headers)
$(CXX) $(CFLAGS) $(eigenIncludes) $(boostIncludes) -DUSE_BUCKINGHAM_POTENTIAL -o execBuck $(sources)
.PHONY: run
run: execLJ execBuck
./execLJ
./execBuck
.PHONY: clean
clean:
-rm -Rf $(toClean)
CXX ?= g++
macroDefines :=
ifdef NDEBUG
$(info Release build)
optimizationFlags := -O3 -march=native
macroDefines := -DNDEBUG
else
$(info Debug build)
optimizationFlags := -O0 -g
endif
ifdef UNIF_PRIOR
$(info Uniform prior)
macroDefines += -DUNIF_PRIOR
else
$(info Gaussian prior)
endif
ifdef USE_FLOATS
$(info using floats)
macroDefines += -DUSE_FLOATS
else
$(info using numeric_ts)
endif
ifdef USE_BRENTS_METHOD
$(info Using Brent's method to find the tho values, needs boost)
macroDefines += -DUSE_BRENTS_METHOD
hostname = $(shell hostname)
ifeq ($(hostname),swmb)
$(info Compiling on Stefano's laptop)
boostIncludes := -I/opt/local/include
endif
ifneq ($(findstring euler,$(hostname)),)
$(info Compiling on EULER, wish me luck!)
$(shell module load boost/1.57.0)
boostIncludes := -I/cluster/apps/boost/1.57.0/x86_64/gcc_4.8.2/serial/include
endif
else
$(info Not using Brent's method)
boostIncludes :=
endif
CFLAGS := -fopenmp -std=c++14 $(optimizationFlags) $(macroDefines)
eigenIncludes := $(shell pkg-config --cflags eigen3)
ifeq ($(eigenIncludes),)
$(shell module load eigen)
eigenIncludes := -I/cluster/apps/eigen/3.2.1/x86_64/gcc_4.8.2/serial/include
endif
sources := $(wildcard *.cpp)
headers := $(wildcard *.hpp)
executable := execLJ execBuck execLJ.dSYM execBuck.dSYM
toClean := $(executable)
.PHONY: all
all: run
execLJ: $(sources) $(headers)
$(CXX) $(CFLAGS) $(eigenIncludes) $(boostIncludes) -DUSE_LENNARD_JONES_POTENTIAL -o execLJ $(sources)
execBuck: $(sources) $(headers)
$(CXX) $(CFLAGS) $(eigenIncludes) $(boostIncludes) -DUSE_BUCKINGHAM_POTENTIAL -o execBuck $(sources)
.PHONY: run
run: execLJ execBuck
./execLJ
./execBuck
.PHONY: clean
clean:
-rm -Rf $(toClean)
#pragma once
#include "typedefs_and_globals.hpp"
static std::normal_distribution<> dist_mvn;
struct normal_random_variable
{
normal_random_variable(Eigen::MatrixXd const& covar)
normal_random_variable(Matrix const& covar)
{
Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigenSolver(covar);
Eigen::SelfAdjointEigenSolver<Matrix> eigenSolver(covar);
transform = eigenSolver.eigenvectors() * eigenSolver.eigenvalues().cwiseSqrt().asDiagonal();
}
Eigen::MatrixXd transform;
Matrix transform;
Eigen::VectorXd operator()() const
Vector operator()() const
{
const int numberOfEntries = transform.cols();
Eigen::VectorXd randomStandardNormal(numberOfEntries);
Vector randomStandardNormal(numberOfEntries);
for (int i = 0; i < numberOfEntries; ++i) {
randomStandardNormal(i) = dist_mvn(GEN);
}
......
......@@ -8,6 +8,8 @@
#include <string>
#include <functional>
#include <cstddef>
#include <cstdint>
static inline numeric_t lennard_jones_potential_function(const numeric_t r, const ThetaVector& theta) {
......@@ -31,7 +33,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 +77,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;
......
......@@ -2,7 +2,14 @@
#include <Eigen/Dense>
#include "xoroshiro.h"
#ifdef USE_FLOATS
#pragma message("Using floats instead of doubles")
typedef float numeric_t;
#else
#pragma message("Using doubles")
typedef double numeric_t;
#endif
typedef unsigned int index_t;
typedef Eigen::Matrix<numeric_t, Eigen::Dynamic, 1> Vector;
typedef Eigen::Matrix<numeric_t, Eigen::Dynamic, Eigen::Dynamic> Matrix;
......@@ -53,7 +60,7 @@ 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;
......
......@@ -10,11 +10,15 @@ template <typename T>
T normal_pdf_1D(T x, T m, T s){
static const T inv_sqrt_2pi = 0.3989422804014327;
T a = (x - m) / s;
if (std::fabs(a) > 10) {
return 1e-16;
}
return inv_sqrt_2pi / s * std::exp(-T(0.5) * a * a);
}
template <typename T>
T log_normal_pdf_1D(T x, T m, T s) {
static const T log_inv_sqrt_2pi = -0.918938533204673;
const T a = (x - m) / s;
return log_inv_sqrt_2pi - std::log(s) - 0.5 * a * a;
}
/// PDF of a Gamma distribution with \alpha = 'a' and \beta = 'b'
template <typename T>
T gamma_pdf_1D(T x, T a, T b){
......@@ -25,7 +29,7 @@ inline static numeric_t mvn_iso(const Vector &x, const Vector &meanVec, const nu
return std::exp(x.binaryExpr(
meanVec,
[=](const numeric_t point, const numeric_t mean) {
return std::log(normal_pdf_1D(point, mean, st_dev));
return log_normal_pdf_1D(point, mean, st_dev);
}
).sum());
}
......@@ -34,7 +38,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<Matrix> Chol;
Chol chol(covMat);
// Handle non positive definite covariance somehow:
#ifndef NDEBUG
......
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