# Sample R program to simulate from normal model, known mean by slice sampling.

# Setup
library('coda')

rm(list=ls())
setwd("~/Data/ScytheMCMC/examples")

X <- rnorm(1000, .25, 1)

# Save table in Scythe format (dense matrix with spaces)
write.table(X, 'build/X.txt', quote=FALSE, row.names=FALSE, col.names=FALSE, sep=" ")

# Compile model, with FAST flag, using g++
system("g++ -DFAST -O3 -funroll-loops -I /usr/local/include/boost-1_39/ -I ../src normal.cpp -o build/normal_fast")
system("g++ -O3 -funroll-loops -I /usr/local/include/boost-1_39/ -I ../src normal.cpp -o build/normal_slow")

# Compile model, with FAST flag, using Intel C++ compiler (does not use MKL). 
system("icpc -DFAST -I /usr/local/include/boost-1_39/ -I ../src normal.cpp -o build/normal_mkl")

# 30 times slower than MKL (FAST)
system("./build/normal_slow --config-file=normal.ini --chains=1 --sample-size=5000 --burnin=5000 --thin=0 --random-seed=123456 --out-file=build/out.txt")
# 3 times slower than MKL (FAST)
system("./build/normal_fast --config-file=normal.ini --chains=1 --sample-size=5000 --burnin=5000 --thin=0 --random-seed=123456 --out-file=build/out.txt")
# Fastest (Still not vectorized, which should lead to 2-4 time speed increase).
system("./build/normal_mkl --config-file=normal.ini --chains=1 --sample-size=5000 --burnin=5000 --thin=0 --random-seed=123456 --out-file=build/out.txt")

result <- read.csv('build/out.txt')
resultmc <- as.mcmc(result)
summary(resultmc)
plot(resultmc)
# Slice sampling, while slower per iteration, often has great mixing properties compared to MH.
autocorr.diag(resultmc)
rejectionRate(resultmc)

# Perfect mixing!
effectiveSize(resultmc)

# Check result:
hist(result$mean)
abline(v=.25)