From 489522226ec39a1353727f2b3eb704337c31c17b Mon Sep 17 00:00:00 2001 From: Keith O'Hara Date: Fri, 2 Feb 2018 02:12:30 -0500 Subject: [PATCH] mostly size_t changes and formatting changes --- src/aees.cpp | 48 ++++++++++++++++++++++++++++++------------------ src/de.cpp | 16 +++++++++------- src/hmc.cpp | 30 +++++++++++++++++++----------- src/mala.cpp | 35 +++++++++++++++++++---------------- src/rmhmc.cpp | 37 ++++++++++++++++++++++--------------- src/rwmh.cpp | 28 +++++++++++++++++++--------- 6 files changed, 118 insertions(+), 76 deletions(-) diff --git a/src/aees.cpp b/src/aees.cpp index ca0369a..a801e90 100644 --- a/src/aees.cpp +++ b/src/aees.cpp @@ -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 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 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); } }; @@ -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; @@ -123,27 +128,32 @@ 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); @@ -151,7 +161,7 @@ mcmc::aees_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio // 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; } @@ -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); } @@ -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) { diff --git a/src/de.cpp b/src/de.cpp index 41ef245..61586d1 100644 --- a/src/de.cpp +++ b/src/de.cpp @@ -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); @@ -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)) { @@ -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; @@ -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; } @@ -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)); } } diff --git a/src/hmc.cpp b/src/hmc.cpp index 50778b0..8c20680 100644 --- a/src/hmc.cpp +++ b/src/hmc.cpp @@ -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); } }; @@ -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); @@ -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; @@ -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 @@ -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; @@ -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(); } @@ -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)); } } diff --git a/src/mala.cpp b/src/mala.cpp index 34b5b03..a51a8f3 100644 --- a/src/mala.cpp +++ b/src/mala.cpp @@ -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); } }; @@ -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); @@ -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; } }; @@ -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(); } @@ -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)); } } diff --git a/src/rmhmc.cpp b/src/rmhmc.cpp index 6557ed9..354b09e 100644 --- a/src/rmhmc.cpp +++ b/src/rmhmc.cpp @@ -56,13 +56,18 @@ mcmc::rmhmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functi // // lambda function for box constraints - std::function box_log_kernel = [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) { + std::function box_log_kernel \ + = [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) + { 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); } }; @@ -83,7 +88,8 @@ mcmc::rmhmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functi // - for (int i=0; i < n_vals; i++) { + for (size_t i=0; i < n_vals; i++) + { arma::mat tmp_mat = inv_tensor_mat * tensor_deriv.slice(i); grad_obj(i) = - grad_obj(i) + 0.5*( arma::trace(tmp_mat) - arma::as_scalar(mntm_inp.t() * tmp_mat * inv_tensor_mat * mntm_inp) ); @@ -107,7 +113,7 @@ mcmc::rmhmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functi // - for (int i=0; i < n_vals; i++) { + for (size_t i=0; i < n_vals; i++) { arma::mat tmp_mat = inv_tensor_mat * tensor_deriv.slice(i); grad_obj(i) = - grad_obj(i) + 0.5*( arma::trace(tmp_mat) - arma::as_scalar(mntm_inp.t() * tmp_mat * inv_tensor_mat * mntm_inp) ); @@ -175,19 +181,19 @@ mcmc::rmhmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functi int n_accept = 0; double comp_val; - 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 = arma::chol(prev_tensor,"lower")*arma::randn(n_vals,1); prev_K = arma::dot(new_mntm,inv_prev_tensor*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 arma::vec prop_mntm = new_mntm; - for (int kk=0; kk < n_fp_steps; kk++) { + for (size_t kk=0; kk < n_fp_steps; kk++) { prop_mntm = new_mntm + mntm_update_fn(new_draw,prop_mntm,target_data,step_size,inv_prev_tensor,prev_deriv_cube,nullptr); // half-step } @@ -198,7 +204,8 @@ mcmc::rmhmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functi arma::vec prop_draw = new_draw; // new_draw += step_size * inv_prev_tensor * new_mntm; - for (int kk=0; kk < n_fp_steps; kk++) { + for (size_t kk=0; kk < n_fp_steps; kk++) + { inv_new_tensor = arma::inv(box_tensor_fn(prop_draw,nullptr,tensor_data)); prop_draw = new_draw + 0.5 * step_size * ( inv_prev_tensor + inv_new_tensor ) * new_mntm; @@ -212,7 +219,6 @@ mcmc::rmhmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functi // new_mntm += mntm_update_fn(new_draw,new_mntm,target_data,step_size,inv_new_tensor,new_deriv_cube,nullptr); // half-step - } prop_U = cons_term - box_log_kernel(new_draw, nullptr, target_data) + 0.5*std::log(arma::det(new_tensor)); @@ -228,7 +234,7 @@ mcmc::rmhmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functi comp_val = std::min(0.0, - (prop_U + prop_K) + (prev_U + prev_K)); double z_rand = arma::as_scalar(arma::randu(1)); - if (z_rand < std::exp(comp_val)) + if (z_rand < std::exp(comp_val)) { prev_draw = new_draw; @@ -239,7 +245,8 @@ mcmc::rmhmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functi inv_prev_tensor = inv_new_tensor; prev_deriv_cube = new_deriv_cube; - if (jj >= n_draws_burnin) { + if (jj >= n_draws_burnin) + { draws_out.row(jj - n_draws_burnin) = new_draw.t(); n_accept++; } @@ -260,7 +267,7 @@ mcmc::rmhmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functi #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)); } } diff --git a/src/rwmh.cpp b/src/rwmh.cpp index 20dfe93..5eb513b 100644 --- a/src/rwmh.cpp +++ b/src/rwmh.cpp @@ -55,13 +55,18 @@ mcmc::rwmh_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio // lambda function for box constraints - std::function 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 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); } }; @@ -91,7 +96,8 @@ mcmc::rwmh_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++) + { new_draw = prev_draw + cov_mcmc_chol * krand.randn(); @@ -106,15 +112,19 @@ mcmc::rwmh_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio double comp_val = std::min(0.0,prop_LP - prev_LP); 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_LP = prop_LP; - if (jj >= n_draws_burnin) { + if (jj >= n_draws_burnin) + { 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(); } @@ -129,7 +139,7 @@ mcmc::rwmh_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)); } }