Skip to content

Commit

Permalink
mostly size_t changes and formatting changes
Browse files Browse the repository at this point in the history
  • Loading branch information
kthohr committed Feb 2, 2018
1 parent ccb086d commit 4895222
Show file tree
Hide file tree
Showing 6 changed files with 118 additions and 76 deletions.
48 changes: 30 additions & 18 deletions src/aees.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -77,13 +77,18 @@ mcmc::aees_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio

// lambda function for box constraints

std::function<double (const arma::vec& vals_inp, void* target_data)> box_log_kernel = [target_log_kernel, vals_bound, bounds_type, lower_bounds, upper_bounds] (const arma::vec& vals_inp, void* target_data) -> double {
//
if (vals_bound) {
std::function<double (const arma::vec& vals_inp, void* target_data)> box_log_kernel \
= [target_log_kernel, vals_bound, bounds_type, lower_bounds, upper_bounds] (const arma::vec& vals_inp, void* target_data) \
-> double
{
if (vals_bound)
{
arma::vec vals_inv_trans = inv_transform(vals_inp, bounds_type, lower_bounds, upper_bounds);

return target_log_kernel(vals_inv_trans, target_data) + log_jacobian(vals_inp, bounds_type, lower_bounds, upper_bounds);
} else {
}
else
{
return target_log_kernel(vals_inp, target_data);
}
};
Expand Down Expand Up @@ -112,7 +117,7 @@ mcmc::aees_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio

double val_out; // holds kernel value

for (int n=0; n < (total_draws); n++) {
for (size_t n=0; n < total_draws; n++) {

arma::mat X_prev = X_new;
kernel_vals_prev = kernel_vals_new;
Expand All @@ -123,35 +128,40 @@ mcmc::aees_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio

// loop down temperature vector

for (int j=1; j < K; j++) {
if (n > j*(n_initial_draws + n_burnin)) {

for (size_t j=1; j < K; j++)
{
if (n > j*(n_initial_draws + n_burnin))
{
double z_eps = arma::as_scalar(arma::randu(1,1));

if (z_eps > ee_prob_par) {
if (z_eps > ee_prob_par)
{
X_new.col(j) = single_step_mh(X_prev.col(j), temper_vec(j), sqrt_cov_mcmc, box_log_kernel, target_data, &val_out);

kernel_vals_new(0,j) = val_out / temper_vec(j-1);
kernel_vals_new(1,j) = val_out / temper_vec(j);
} else {

}
else
{
int initial_j = (j-1)*(n_initial_draws + n_burnin);

int ring_ind_spacing = std::floor( (n - initial_j + 1) / n_rings);

if (ring_ind_spacing == 0) {
if (ring_ind_spacing == 0)
{
X_new.col(j) = X_prev.col(j);
kernel_vals_new.col(j) = kernel_vals_prev.col(j);
} else {

}
else
{
arma::vec past_kernel_vals = arma::trans(kernel_vals(arma::span(j-1,j-1),arma::span(initial_j,n)));

arma::uvec sort_ind = arma::sort_index(past_kernel_vals);
past_kernel_vals = past_kernel_vals(sort_ind);

// construct rings

for (int i=0; i < (n_rings-1); i++) {
for (size_t i=0; i < (n_rings-1); i++) {
ring_vals(j-1,i) = (past_kernel_vals((i+1)*ring_ind_spacing) + past_kernel_vals((i+1)*ring_ind_spacing-1)) / 2.0;
}

Expand Down Expand Up @@ -183,7 +193,8 @@ mcmc::aees_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio
double comp_val = std::min(0.0, comp_val_1 + comp_val_2);
double z = arma::as_scalar(arma::randu(1,1));

if (z > std::exp(comp_val)) {
if (z > std::exp(comp_val))
{
X_new.col(j) = X_prev.col(j);
kernel_vals_new.col(j) = kernel_vals_prev.col(j);
}
Expand All @@ -203,7 +214,8 @@ mcmc::aees_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio

draws_out.set_size(total_draws,n_vals);

for (int i=0; i < total_draws; i++) {
for (size_t i=0; i < total_draws; i++)
{
arma::vec tmp_vec = X_out(arma::span(),arma::span(K-1,K-1),arma::span(i,i));

if (vals_bound) {
Expand Down
16 changes: 9 additions & 7 deletions src/de.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -79,7 +79,7 @@ mcmc::de_int(const arma::vec& initial_vals, arma::cube& draws_out, std::function
#ifdef MCMC_USE_OMP
#pragma omp parallel for
#endif
for (int i=0; i < n_pop; i++) {
for (size_t i=0; i < n_pop; i++) {
X.row(i) = par_initial_lb.t() + (par_initial_ub.t() - par_initial_lb.t())%arma::randu(1,n_vals);

double prop_kernel_val = box_log_kernel(X.row(i).t(),target_data);
Expand All @@ -99,7 +99,7 @@ mcmc::de_int(const arma::vec& initial_vals, arma::cube& draws_out, std::function
int n_accept = 0;
double par_gamma_run = par_gamma;

for (int j=0; j < n_gen + n_burnin; j++) {
for (size_t j=0; j < n_gen + n_burnin; j++) {
double temperature_j = de_cooling_schedule(j,n_gen);

if (jumps && ((j+1) % 10 == 0)) {
Expand All @@ -109,7 +109,7 @@ mcmc::de_int(const arma::vec& initial_vals, arma::cube& draws_out, std::function
#ifdef MCMC_USE_OMP
#pragma omp parallel for
#endif
for (int i=0; i < n_pop; i++) {
for (size_t i=0; i < n_pop; i++) {

int R_1, R_2;

Expand Down Expand Up @@ -148,11 +148,13 @@ mcmc::de_int(const arma::vec& initial_vals, arma::cube& draws_out, std::function
}
}
}

//
if(j >= n_burnin){

if (j >= n_burnin) {
draws_out.slice(j-n_burnin) = X;
}
//

if (jumps && ((j+1) % 10 == 0)) {
par_gamma_run = par_gamma;
}
Expand All @@ -166,8 +168,8 @@ mcmc::de_int(const arma::vec& initial_vals, arma::cube& draws_out, std::function
#ifdef MCMC_USE_OMP
#pragma omp parallel for
#endif
for (int ii = 0; ii < n_gen; ii++) {
for (int jj = 0; jj < n_pop; jj++) {
for (size_t ii = 0; ii < n_gen; ii++) {
for (size_t jj = 0; jj < n_pop; jj++) {
draws_out.slice(ii).row(jj) = arma::trans(inv_transform(draws_out.slice(ii).row(jj).t(), bounds_type, lower_bounds, upper_bounds));
}
}
Expand Down
30 changes: 19 additions & 11 deletions src/hmc.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -63,11 +63,14 @@ mcmc::hmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function
= [target_log_kernel, vals_bound, bounds_type, lower_bounds, upper_bounds] (const arma::vec& vals_inp, arma::vec* grad_out, void* target_data) \
-> double
{
if (vals_bound) {
if (vals_bound)
{
arma::vec vals_inv_trans = inv_transform(vals_inp, bounds_type, lower_bounds, upper_bounds);

return target_log_kernel(vals_inv_trans, nullptr, target_data) + log_jacobian(vals_inp, bounds_type, lower_bounds, upper_bounds);
} else {
}
else
{
return target_log_kernel(vals_inp, nullptr, target_data);
}
};
Expand All @@ -80,8 +83,8 @@ mcmc::hmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function
const int n_vals = pos_inp.n_elem;
arma::vec grad_obj(n_vals);

if (vals_bound) {

if (vals_bound)
{
arma::vec pos_inv_trans = inv_transform(pos_inp, bounds_type, lower_bounds, upper_bounds);

target_log_kernel(pos_inv_trans,&grad_obj,target_data);
Expand All @@ -97,7 +100,9 @@ mcmc::hmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function
//

return mntm_inp + step_size * jacob_matrix * grad_obj / 2.0;
} else {
}
else
{
target_log_kernel(pos_inp,&grad_obj,target_data);

return mntm_inp + step_size * grad_obj / 2.0;
Expand Down Expand Up @@ -129,14 +134,14 @@ mcmc::hmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function

int n_accept = 0;

for (int jj = 0; jj < n_draws_keep + n_draws_burnin; jj++) {

for (size_t jj = 0; jj < n_draws_keep + n_draws_burnin; jj++)
{
new_mntm = sqrt_precond_matrix*arma::randn(n_vals,1);
prev_K = arma::dot(new_mntm,inv_precond_matrix*new_mntm) / 2.0;

new_draw = prev_draw;

for (int k = 0; k < n_leap_steps; k++)
for (size_t k = 0; k < n_leap_steps; k++)
{ // begin leap frog steps

new_mntm = mntm_update_fn(new_draw,new_mntm,target_data,step_size,nullptr); // half-step
Expand All @@ -163,7 +168,8 @@ mcmc::hmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function
double comp_val = std::min(0.0,- prop_U - prop_K + prev_U + prev_K);
double z = arma::as_scalar(arma::randu(1));

if (z < std::exp(comp_val)) {
if (z < std::exp(comp_val))
{
prev_draw = new_draw;
prev_U = prop_U;
prev_K = prop_K;
Expand All @@ -172,7 +178,9 @@ mcmc::hmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function
draws_out.row(jj - n_draws_burnin) = new_draw.t();
n_accept++;
}
} else {
}
else
{
if (jj >= n_draws_burnin) {
draws_out.row(jj - n_draws_burnin) = prev_draw.t();
}
Expand All @@ -187,7 +195,7 @@ mcmc::hmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function
#ifdef MCMC_USE_OMP
#pragma omp parallel for
#endif
for (int jj = 0; jj < n_draws_keep; jj++) {
for (size_t jj = 0; jj < n_draws_keep; jj++) {
draws_out.row(jj) = arma::trans(inv_transform(draws_out.row(jj).t(), bounds_type, lower_bounds, upper_bounds));
}
}
Expand Down
35 changes: 19 additions & 16 deletions src/mala.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -61,12 +61,14 @@ mcmc::mala_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio
= [target_log_kernel, vals_bound, bounds_type, lower_bounds, upper_bounds] (const arma::vec& vals_inp, arma::vec* grad_out, void* target_data) \
-> double
{
//
if (vals_bound) {
if (vals_bound)
{
arma::vec vals_inv_trans = inv_transform(vals_inp, bounds_type, lower_bounds, upper_bounds);

return target_log_kernel(vals_inv_trans, nullptr, target_data) + log_jacobian(vals_inp, bounds_type, lower_bounds, upper_bounds);
} else {
}
else
{
return target_log_kernel(vals_inp, nullptr, target_data);
}
};
Expand All @@ -79,8 +81,8 @@ mcmc::mala_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio
const int n_vals = vals_inp.n_elem;
arma::vec grad_obj(n_vals);

if (vals_bound) {

if (vals_bound)
{
arma::vec vals_inv_trans = inv_transform(vals_inp, bounds_type, lower_bounds, upper_bounds);

target_log_kernel(vals_inv_trans,&grad_obj,target_data);
Expand All @@ -95,15 +97,13 @@ mcmc::mala_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio

//

arma::vec prop_vals_out = vals_inp + step_size * step_size * jacob_matrix * precond_matrix * grad_obj / 2.0;

return prop_vals_out;
} else {
return vals_inp + step_size * step_size * jacob_matrix * precond_matrix * grad_obj / 2.0;
}
else
{
target_log_kernel(vals_inp,&grad_obj,target_data);

arma::vec prop_vals_out = vals_inp + step_size * step_size * precond_matrix * grad_obj / 2.0;

return prop_vals_out;
return vals_inp + step_size * step_size * precond_matrix * grad_obj / 2.0;
}
};

Expand All @@ -129,14 +129,17 @@ mcmc::mala_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio
int n_accept = 0;
arma::vec krand(n_vals);

for (int jj = 0; jj < n_draws_keep + n_draws_burnin; jj++)
for (size_t jj = 0; jj < n_draws_keep + n_draws_burnin; jj++)
{
if (vals_bound) {
if (vals_bound)
{
arma::mat jacob_matrix;
arma::vec mean_vec = mala_mean_fn(prev_draw, target_data, step_size, &jacob_matrix);

new_draw = mean_vec + step_size * arma::chol(jacob_matrix,"lower") * sqrt_precond_matrix * krand.randn();
} else {
}
else
{
new_draw = mala_mean_fn(prev_draw, target_data, step_size, nullptr) + step_size * sqrt_precond_matrix * krand.randn();
}

Expand Down Expand Up @@ -179,7 +182,7 @@ mcmc::mala_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio
#ifdef MCMC_USE_OMP
#pragma omp parallel for
#endif
for (int jj = 0; jj < n_draws_keep; jj++) {
for (size_t jj = 0; jj < n_draws_keep; jj++) {
draws_out.row(jj) = arma::trans(inv_transform(draws_out.row(jj).t(), bounds_type, lower_bounds, upper_bounds));
}
}
Expand Down
Loading

0 comments on commit 4895222

Please sign in to comment.