update base library and build settings

kthohr · Feb 6, 2018 · 96ceb54 · 96ceb54
1 parent 3bf587f
commit 96ceb54
Show file tree

Hide file tree

Showing 241 changed files with 4,857 additions and 2,321 deletions.
diff --git a/.gitignore b/.gitignore
@@ -6,3 +6,4 @@
 *.vscode
 /src/bmlib/include/armadillo*
 .Rproj.user
+*.tmp
diff --git a/src/Makevars b/src/Makevars
@@ -15,7 +15,7 @@ else
     BMR_OPENMP=
 endif
 
-PKG_CPPFLAGS= $(CXX11STD) $(BMR_OPENMP) -DUSE_RCPP_ARMADILLO -DARMA_NO_DEBUG -DSTATSLIB_GO_INLINE -I./modules -I./bmlib/include -I./bmlib/include/mcmc -I./bmlib/include/optim
+PKG_CPPFLAGS= $(CXX11STD) $(BMR_OPENMP) -DUSE_RCPP_ARMADILLO -DARMA_NO_DEBUG -DSTATS_GO_INLINE -I./modules -I./bmlib/include -I./bmlib/include/mcmc -I./bmlib/include/optim
 PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
 
 CXX_STD=CXX11STD

diff --git a/src/Makevars.win b/src/Makevars.win
@@ -8,7 +8,7 @@ BM_HEADER_DIR = $(BM_DIR)/include
 #BMR_OPENMP=$(SHLIB_OPENMP_CXXFLAGS) -DARMA_USE_OPENMP
 BMR_OPENMP=$(SHLIB_OPENMP_CXXFLAGS)
 
-PKG_CPPFLAGS= $(CXX11STD) $(BMR_OPENMP) -DUSE_RCPP_ARMADILLO -DARMA_NO_DEBUG -DSTATSLIB_GO_INLINE -I./modules -I./bmlib/include -I./bmlib/include/mcmc -I./bmlib/include/optim
+PKG_CPPFLAGS= $(CXX11STD) $(BMR_OPENMP) -DUSE_RCPP_ARMADILLO -DARMA_NO_DEBUG -DSTATS_GO_INLINE -I./modules -I./bmlib/include -I./bmlib/include/mcmc -I./bmlib/include/optim
 PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS)
 
 CXX_STD=CXX11STD

diff --git a/src/bmlib/include/dsge/dsge_class.hpp b/src/bmlib/include/dsge/dsge_class.hpp
@@ -95,10 +95,10 @@ class dsge
 
         arma::vec estim_mode(const arma::vec& initial_vals);
         arma::vec estim_mode(const arma::vec& initial_vals, arma::mat& vcov_mat);
-        arma::vec estim_mode(const arma::vec& initial_vals, arma::mat* vcov_mat, optim::opt_settings* settings_inp);
+        arma::vec estim_mode(const arma::vec& initial_vals, arma::mat* vcov_mat, optim::algo_settings* settings_inp);
 
         void estim_mcmc(const arma::vec& initial_vals);
-        void estim_mcmc(const arma::vec& initial_vals, mcmc::mcmc_settings* settings_inp);
+        void estim_mcmc(const arma::vec& initial_vals, mcmc::algo_settings* settings_inp);
 
         arma::cube IRF(const int n_irf_periods, const bool observ_irfs) const;
 

diff --git a/src/bmlib/include/dsge/dsge_class.tpp b/src/bmlib/include/dsge/dsge_class.tpp
@@ -201,15 +201,15 @@ dsge<T>::estim_mode(const arma::vec& initial_vals, arma::mat& vcov_mat)
 
 template<typename T>
 arma::vec
-dsge<T>::estim_mode(const arma::vec& initial_vals, arma::mat* vcov_mat, optim::opt_settings* settings_inp)
+dsge<T>::estim_mode(const arma::vec& initial_vals, arma::mat* vcov_mat, optim::algo_settings* settings_inp)
 {
     dsge_estim_data<T> mode_data;
     mode_data.dsge_obj = *this;
 
     //
     // optimization settings
 
-    optim::opt_settings settings;
+    optim::algo_settings settings;
 
     if (settings_inp) {
         settings = *settings_inp;
@@ -256,7 +256,7 @@ dsge<T>::mcmc_objfn(const arma::vec& pars_inp, void* mcmc_data)
 
 template<typename T>
 void
-dsge<T>::estim_mcmc(const arma::vec& initial_vals, mcmc::mcmc_settings* settings_inp)
+dsge<T>::estim_mcmc(const arma::vec& initial_vals, mcmc::algo_settings* settings_inp)
 {
 
     dsge_estim_data<T> mcmc_data;
@@ -265,7 +265,7 @@ dsge<T>::estim_mcmc(const arma::vec& initial_vals, mcmc::mcmc_settings* settings
     //
     // MCMC settings
 
-    mcmc::mcmc_settings settings;
+    mcmc::algo_settings settings;
 
     if (settings_inp) {
         settings = *settings_inp;

diff --git a/src/bmlib/include/dsge/dsgevar_class.hpp b/src/bmlib/include/dsge/dsgevar_class.hpp
@@ -88,9 +88,9 @@ class dsgevar
 
         arma::vec estim_mode(const arma::vec& initial_vals);
         arma::vec estim_mode(const arma::vec& initial_vals, arma::mat& vcov_mat);
-        arma::vec estim_mode(const arma::vec& initial_vals, arma::mat* vcov_mat, optim::opt_settings* settings_inp);
+        arma::vec estim_mode(const arma::vec& initial_vals, arma::mat* vcov_mat, optim::algo_settings* settings_inp);
 
-        void estim_mcmc(const arma::vec& initial_vals, mcmc::mcmc_settings* settings_inp);
+        void estim_mcmc(const arma::vec& initial_vals, mcmc::algo_settings* settings_inp);
 
         arma::cube IRF(const int n_irf_periods) const;
 

diff --git a/src/bmlib/include/dsge/dsgevar_class.tpp b/src/bmlib/include/dsge/dsgevar_class.tpp
@@ -269,15 +269,15 @@ dsgevar<T>::estim_mode(const arma::vec& initial_vals, arma::mat& vcov_mat)
 
 template<typename T>
 arma::vec
-dsgevar<T>::estim_mode(const arma::vec& initial_vals, arma::mat* vcov_mat, optim::opt_settings* settings_inp)
+dsgevar<T>::estim_mode(const arma::vec& initial_vals, arma::mat* vcov_mat, optim::algo_settings* settings_inp)
 {
     dsgevar_estim_data<T> mode_data;
     mode_data.dsgevar_obj = *this;
 
     //
     // optimization settings
 
-    optim::opt_settings settings;
+    optim::algo_settings settings;
 
     if (settings_inp) {
         settings = *settings_inp;
@@ -324,7 +324,7 @@ dsgevar<T>::mcmc_objfn(const arma::vec& pars_inp, void* mcmc_data)
 
 template<typename T>
 void
-dsgevar<T>::estim_mcmc(const arma::vec& initial_vals, mcmc::mcmc_settings* settings_inp)
+dsgevar<T>::estim_mcmc(const arma::vec& initial_vals, mcmc::algo_settings* settings_inp)
 {
 
     dsgevar_estim_data<T> mcmc_data;
@@ -333,7 +333,7 @@ dsgevar<T>::estim_mcmc(const arma::vec& initial_vals, mcmc::mcmc_settings* setti
     //
     // MCMC settings
 
-    mcmc::mcmc_settings settings;
+    mcmc::algo_settings settings;
 
     if (settings_inp) {
         settings = *settings_inp;

diff --git a/src/bmlib/include/mcmc/aees.hpp b/src/bmlib/include/mcmc/aees.hpp
@@ -0,0 +1,64 @@
+/*################################################################################
+  ##
+  ##   Copyright (C) 2011-2018 Keith O'Hara
+  ##
+  ##   This file is part of the MCMC C++ library.
+  ##
+  ##   MCMC is free software: you can redistribute it and/or modify
+  ##   it under the terms of the GNU General Public License as published by
+  ##   the Free Software Foundation, either version 2 of the License, or
+  ##   (at your option) any later version.
+  ##
+  ##   MCMC is distributed in the hope that it will be useful,
+  ##   but WITHOUT ANY WARRANTY; without even the implied warranty of
+  ##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  ##   GNU General Public License for more details.
+  ##
+  ################################################################################*/
+
+/*
+ * Adaptive Equi-Energy Sampler
+ */
+
+#ifndef _mcmc_aees_HPP
+#define _mcmc_aees_HPP
+
+bool aees_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, void* target_data)> target_log_kernel, void* target_data, algo_settings* settings_inp);
+
+bool aees(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, void* target_data)> target_log_kernel, void* target_data);
+bool aees(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, void* target_data)> target_log_kernel, void* target_data, algo_settings& settings);
+
+// single-step Metropolis-Hastings for tempered distributions
+inline
+arma::vec
+single_step_mh(const arma::vec& X_prev, const double temper_val, const arma::mat& sqrt_cov_mcmc, std::function<double (const arma::vec& vals_inp, void* target_data)> target_log_kernel, void* target_data, double* val_out)
+{
+    const int n_vals = X_prev.n_elem;
+
+    arma::vec X_new = X_prev + std::sqrt(temper_val)*sqrt_cov_mcmc*arma::randn(n_vals,1);
+
+    //
+
+    double val_new  = target_log_kernel(X_new, target_data);
+    double val_prev = target_log_kernel(X_prev, target_data);
+
+    double comp_val = std::min(0.0, (val_new - val_prev) / temper_val );
+
+    double z = arma::as_scalar(arma::randu(1,1));
+
+    if (z < std::exp(comp_val)){
+        if (val_out) {
+            *val_out = val_new;
+        }
+
+        return X_new;
+    } else {
+        if (val_out) {
+            *val_out = val_prev;
+        }
+
+        return X_prev;
+    }
+}
+
+#endif
diff --git a/src/bmlib/include/mcmc/de.hpp b/src/bmlib/include/mcmc/de.hpp
@@ -18,24 +18,18 @@
 
 /*
  * Differential Evolution (DE) MCMC
- *
- * Keith O'Hara
- * 03/01/2016
- *
- * This version:
- * 08/12/2017
  */
 
 #ifndef _mcmc_de_HPP
 #define _mcmc_de_HPP
 
-bool de_int(const arma::vec& initial_vals, arma::cube& draws_out, std::function<double (const arma::vec& vals_inp, void* target_data)> target_log_kernel, void* target_data, mcmc_settings* settings_inp);
+bool de_int(const arma::vec& initial_vals, arma::cube& draws_out, std::function<double (const arma::vec& vals_inp, void* target_data)> target_log_kernel, void* target_data, algo_settings* settings_inp);
 
 bool de(const arma::vec& initial_vals, arma::cube& draws_out, std::function<double (const arma::vec& vals_inp, void* target_data)> target_log_kernel, void* target_data);
-bool de(const arma::vec& initial_vals, arma::cube& draws_out, std::function<double (const arma::vec& vals_inp, void* target_data)> target_log_kernel, void* target_data, mcmc_settings& settings);
+bool de(const arma::vec& initial_vals, arma::cube& draws_out, std::function<double (const arma::vec& vals_inp, void* target_data)> target_log_kernel, void* target_data, algo_settings& settings);
 
 inline
-double 
+double
 de_cooling_schedule(const int s, const int N_gen)
 {
     return 1.0;

diff --git a/src/bmlib/include/mcmc/determine_bounds_type.hpp b/src/bmlib/include/mcmc/determine_bounds_type.hpp
@@ -18,12 +18,6 @@
 
 /*
  * Determine the upper-lower bounds combo type
- *
- * Keith O'Hara
- * 05/01/2012
- *
- * This version:
- * 08/12/2017
  */
 
 // note: std::isfinite is not true for: NaN, -Inf, or +Inf

diff --git a/src/bmlib/include/mcmc/hmc.hpp b/src/bmlib/include/mcmc/hmc.hpp
@@ -0,0 +1,34 @@
+/*################################################################################
+  ##
+  ##   Copyright (C) 2011-2018 Keith O'Hara
+  ##
+  ##   This file is part of the MCMC C++ library.
+  ##
+  ##   MCMC is free software: you can redistribute it and/or modify
+  ##   it under the terms of the GNU General Public License as published by
+  ##   the Free Software Foundation, either version 2 of the License, or
+  ##   (at your option) any later version.
+  ##
+  ##   MCMC is distributed in the hope that it will be useful,
+  ##   but WITHOUT ANY WARRANTY; without even the implied warranty of
+  ##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  ##   GNU General Public License for more details.
+  ##
+  ################################################################################*/
+
+/*
+ * Hamiltonian Monte Carlo
+ */
+
+ #ifndef _mcmc_hmc_HPP
+ #define _mcmc_hmc_HPP
+
+ bool hmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* target_data)> target_log_kernel, void* target_data, algo_settings* settings_inp);
+
+ bool hmc(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* target_data)> target_log_kernel, void* target_data);
+ bool hmc(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* target_data)> target_log_kernel, void* target_data, algo_settings& settings);
+
+ //
+
+ #endif
+
diff --git a/src/bmlib/include/mcmc/inv_jacobian_adjust.hpp b/src/bmlib/include/mcmc/inv_jacobian_adjust.hpp
@@ -0,0 +1,48 @@
+/*################################################################################
+  ##
+  ##   Copyright (C) 2011-2018 Keith O'Hara
+  ##
+  ##   This file is part of the MCMC C++ library.
+  ##
+  ##   MCMC is free software: you can redistribute it and/or modify
+  ##   it under the terms of the GNU General Public License as published by
+  ##   the Free Software Foundation, either version 2 of the License, or
+  ##   (at your option) any later version.
+  ##
+  ##   MCMC is distributed in the hope that it will be useful,
+  ##   but WITHOUT ANY WARRANTY; without even the implied warranty of
+  ##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  ##   GNU General Public License for more details.
+  ##
+  ################################################################################*/
+
+/*
+ * inverse Jacobian adjustment
+ */
+
+inline
+arma::mat
+inv_jacobian_adjust(const arma::vec& vals_trans_inp, const arma::uvec& bounds_type, const arma::vec& lower_bounds, const arma::vec& upper_bounds)
+{
+    const int n_vals = bounds_type.n_elem;
+
+    arma::mat ret_mat = arma::eye(n_vals,n_vals);
+
+    for (int i=0; i < n_vals; i++) {
+        switch (bounds_type(i)) {
+            case 2: // lower bound only
+                ret_mat(i,i) = 1.0 / std::exp(vals_trans_inp(i));
+                break;
+            case 3: // upper bound only
+                ret_mat(i,i) = 1.0 / std::exp(-vals_trans_inp(i));
+                break;
+            case 4: // upper and lower bounds
+                ret_mat(i,i) = 1.0 / ( std::exp(vals_trans_inp(i))*(upper_bounds(i) - lower_bounds(i)) / std::pow(std::exp(vals_trans_inp(i)) + 1,2) );
+                break;
+        }
+    }
+
+    //
+
+    return ret_mat;
+}
diff --git a/src/bmlib/include/mcmc/log_jacobian.hpp b/src/bmlib/include/mcmc/log_jacobian.hpp
@@ -18,12 +18,6 @@
 
 /*
  * log Jacobian adjustment
- *
- * Keith O'Hara
- * 05/01/2012
- *
- * This version:
- * 08/12/2017
  */
 
 inline

diff --git a/src/bmlib/include/mcmc/mala.hpp b/src/bmlib/include/mcmc/mala.hpp
@@ -0,0 +1,64 @@
+/*################################################################################
+  ##
+  ##   Copyright (C) 2011-2018 Keith O'Hara
+  ##
+  ##   This file is part of the MCMC C++ library.
+  ##
+  ##   MCMC is free software: you can redistribute it and/or modify
+  ##   it under the terms of the GNU General Public License as published by
+  ##   the Free Software Foundation, either version 2 of the License, or
+  ##   (at your option) any later version.
+  ##
+  ##   MCMC is distributed in the hope that it will be useful,
+  ##   but WITHOUT ANY WARRANTY; without even the implied warranty of
+  ##   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  ##   GNU General Public License for more details.
+  ##
+  ################################################################################*/
+
+/*
+ * Metropolis-adjusted Langevin algorithm
+ */
+
+#ifndef _mcmc_mala_HPP
+#define _mcmc_mala_HPP
+
+bool mala_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* target_data)> target_log_kernel, void* target_data, algo_settings* settings_inp);
+
+bool mala(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* target_data)> target_log_kernel, void* target_data);
+bool mala(const arma::vec& initial_vals, arma::mat& draws_out, std::function<double (const arma::vec& vals_inp, arma::vec* grad_out, void* target_data)> target_log_kernel, void* target_data, algo_settings& settings);
+
+//
+
+inline
+double
+mala_prop_adjustment(const arma::vec& prop_vals, const arma::vec& prev_vals, const double step_size, const bool vals_bound, const arma::mat& precond_mat, std::function<arma::vec (const arma::vec& vals_inp, void* target_data, const double step_size, arma::mat* jacob_matrix_out)> mala_mean_fn, void* target_data)
+{
+
+    double ret_val = 0;
+    const double step_size_sq = step_size*step_size;
+
+    //
+
+    if (vals_bound) {
+
+        arma::mat prop_inv_jacob, prev_inv_jacob;
+
+        arma::vec prop_mean = mala_mean_fn(prop_vals, target_data, step_size, &prop_inv_jacob);
+        arma::vec prev_mean = mala_mean_fn(prev_vals, target_data, step_size, &prev_inv_jacob);
+
+        ret_val = stats_mcmc::dmvnorm(prev_vals, prop_mean, step_size_sq*prop_inv_jacob*precond_mat, true) - stats_mcmc::dmvnorm(prop_vals, prev_mean, step_size_sq*prop_inv_jacob*precond_mat, true);
+
+    } else {
+        arma::vec prop_mean = mala_mean_fn(prop_vals, target_data, step_size, nullptr);
+        arma::vec prev_mean = mala_mean_fn(prev_vals, target_data, step_size, nullptr);
+
+        ret_val = stats_mcmc::dmvnorm(prev_vals, prop_mean, step_size_sq*precond_mat, true) - stats_mcmc::dmvnorm(prop_vals, prev_mean, step_size_sq*precond_mat, true);
+    }
+
+    //
+
+    return ret_val;
+}
+
+#endif