Commit 1c374b8e authored by stefanow's avatar stefanow
Browse files

Merge branch 'chris' of gitlab.ethz.ch:stefanow/tmcmc into chris

parents 52c5ad30 3361430c
......@@ -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 that returns a sample of a 'N_DIM'-dimensional Gaussian with mean zero and covariance matrix
......
......@@ -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) {
const numeric_t sigma_div_r_to_the_sixth = std::pow(theta(1) / r, 6);
......
......@@ -38,7 +38,7 @@ int main(){
WeightVector likelihoods(POP_SIZE);
WeightVector likelihoods_old(POP_SIZE);
CountVector resample_count(POP_SIZE);
std::vector<WeightVector> best_thetas;
std::vector<ThetaVector> best_thetas;
CovMat sample_cov_mat;
ThetaVector sample_mean;
numeric_t evidence_s = 1.0;
......
......@@ -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;
......@@ -38,17 +45,17 @@ struct PotentialName<Potential::BUCKINGHAM> {
static constexpr const char* value = "Buckingham";
};
constexpr index_t POP_SIZE = 100000;
constexpr index_t POP_SIZE = 100;
constexpr index_t N_DIM = (POTENTIAL == Potential::LENNARD_JONES? 2: 3) + 1;// + 1 for the unknown standard deviation
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);
constexpr index_t EIGEN_NUM_ROWS_FOR_POP = (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, EIGEN_NUM_ROWS_FOR_POP> PopMatrix;
typedef Eigen::Matrix<numeric_t, N_DIM, -1> PartPopMatrix;
typedef Eigen::Matrix<numeric_t, N_DIM, EIGEN_NUM_ROWS_FOR_POP> PartPopMatrix;
typedef numeric_t (*PotentialFunction)(const numeric_t r, const ThetaVector& theta);
......
......@@ -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());
}
......
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