MPI warmup auto

src/stan/analyze/mcmc/autocovariance.hpp

@@ -1,7 +1,7 @@
 #ifndef STAN_ANALYZE_MCMC_AUTOCOVARIANCE_HPP
 #define STAN_ANALYZE_MCMC_AUTOCOVARIANCE_HPP
 
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/fun.hpp>
 #include <boost/accumulators/accumulators.hpp>
 #include <boost/accumulators/statistics/stats.hpp>
 #include <boost/accumulators/statistics/variance.hpp>

src/stan/analyze/mcmc/compute_potential_scale_reduction.hpp

@@ -1,7 +1,7 @@
 #ifndef STAN_ANALYZE_MCMC_COMPUTE_POTENTIAL_SCALE_REDUCTION_HPP
 #define STAN_ANALYZE_MCMC_COMPUTE_POTENTIAL_SCALE_REDUCTION_HPP
 
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/fun.hpp>
 #include <stan/analyze/mcmc/autocovariance.hpp>
 #include <stan/analyze/mcmc/split_chains.hpp>
 #include <boost/accumulators/accumulators.hpp>

src/stan/io/dump.hpp

@@ -3,7 +3,7 @@
 
 #include <stan/io/validate_zero_buf.hpp>
 #include <stan/io/var_context.hpp>
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/fun.hpp>
 #include <boost/lexical_cast.hpp>
 #include <boost/throw_exception.hpp>
 #include <boost/type_traits/is_floating_point.hpp>

src/stan/io/reader.hpp

@@ -2,7 +2,7 @@
 #define STAN_IO_READER_HPP
 
 #include <boost/throw_exception.hpp>
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/fun.hpp>
 #include <stdexcept>
 #include <string>
 #include <vector>

src/stan/io/stan_csv_reader.hpp

@@ -2,7 +2,7 @@
 #define STAN_IO_STAN_CSV_READER_HPP
 
 #include <boost/algorithm/string.hpp>
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
 #include <istream>
 #include <iostream>
 #include <sstream>

src/stan/io/writer.hpp

@@ -1,7 +1,7 @@
 #ifndef STAN_IO_WRITER_HPP
 #define STAN_IO_WRITER_HPP
 
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/fun.hpp>
 #include <stdexcept>
 #include <vector>
 

src/stan/mcmc/chains.hpp

@@ -2,7 +2,7 @@
 #define STAN_MCMC_CHAINS_HPP
 
 #include <stan/io/stan_csv_reader.hpp>
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/fun.hpp>
 #include <stan/analyze/mcmc/compute_effective_sample_size.hpp>
 #include <stan/analyze/mcmc/compute_potential_scale_reduction.hpp>
 #include <boost/accumulators/accumulators.hpp>

src/stan/mcmc/covar_adaptation.hpp

@@ -1,7 +1,7 @@
 #ifndef STAN_MCMC_COVAR_ADAPTATION_HPP
 #define STAN_MCMC_COVAR_ADAPTATION_HPP
 
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/fun.hpp>
 #include <stan/mcmc/windowed_adaptation.hpp>
 #include <vector>
 

src/stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp

@@ -0,0 +1,76 @@
+#ifndef STAN_MCMC_HMC_HAMILTONIANS_AUTO_E_METRIC_HPP
+#define STAN_MCMC_HMC_HAMILTONIANS_AUTO_E_METRIC_HPP
+
+#include <stan/callbacks/logger.hpp>
+#include <stan/math/prim/fun.hpp>
+#include <stan/mcmc/hmc/hamiltonians/base_hamiltonian.hpp>
+#include <stan/mcmc/hmc/hamiltonians/auto_e_point.hpp>
+#include <boost/random/variate_generator.hpp>
+#include <boost/random/normal_distribution.hpp>
+#include <Eigen/Cholesky>
+
+namespace stan {
+  namespace mcmc {
+
+    // Euclidean manifold with dense metric
+    template <class Model, class BaseRNG>
+    class auto_e_metric
+      : public base_hamiltonian<Model, auto_e_point, BaseRNG> {
+    public:
+      explicit auto_e_metric(const Model& model)
+        : base_hamiltonian<Model, auto_e_point, BaseRNG>(model) {}
+
+      double T(auto_e_point& z) {
+        return 0.5 * z.p.transpose() * z.inv_e_metric_ * z.p;
+      }
+
+      double tau(auto_e_point& z) {
+        return T(z);
+      }
+
+      double phi(auto_e_point& z) {
+        return this->V(z);
+      }
+
+      double dG_dt(auto_e_point& z, callbacks::logger& logger) {
+        return 2 * T(z) - z.q.dot(z.g);
+      }
+
+      Eigen::VectorXd dtau_dq(auto_e_point& z, callbacks::logger& logger) {
+        return Eigen::VectorXd::Zero(this->model_.num_params_r());
+      }
+
+      Eigen::VectorXd dtau_dp(auto_e_point& z) {
+	if(z.is_diagonal_) {
+	  return z.inv_e_metric_.diagonal().cwiseProduct(z.p);
+	} else {
+	  return z.inv_e_metric_ * z.p;
+	}
+      }
+
+      Eigen::VectorXd dphi_dq(auto_e_point& z, callbacks::logger& logger) {
+        return z.g;
+      }
+
+      void sample_p(auto_e_point& z, BaseRNG& rng) {
+        typedef typename stan::math::index_type<Eigen::VectorXd>::type idx_t;
+        boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
+          rand_gaus(rng, boost::normal_distribution<>());
+
+	if(z.is_diagonal_) {
+	  for (int i = 0; i < z.p.size(); ++i)
+	    z.p(i) = rand_gaus() / sqrt(z.inv_e_metric_(i, i));
+	} else {
+	  Eigen::VectorXd u(z.p.size());
+
+	  for (idx_t i = 0; i < u.size(); ++i)
+	    u(i) = rand_gaus();
+
+	  z.p = z.inv_e_metric_.llt().matrixU().solve(u);
+	}
+      }
+    };
+
+  }  // mcmc
+}  // stan
+#endif

src/stan/mcmc/hmc/hamiltonians/auto_e_point.hpp

@@ -0,0 +1,69 @@
+#ifndef STAN_MCMC_HMC_HAMILTONIANS_AUTO_E_POINT_HPP
+#define STAN_MCMC_HMC_HAMILTONIANS_AUTO_E_POINT_HPP
+
+#include <stan/callbacks/writer.hpp>
+#include <stan/mcmc/hmc/hamiltonians/ps_point.hpp>
+
+namespace stan {
+  namespace mcmc {
+    /**
+     * Point in a phase space with a base
+     * Euclidean manifold with auto metric
+     */
+    class auto_e_point: public ps_point {
+    public:
+      /** 
+       * Inverse mass matrix.
+       */
+      Eigen::MatrixXd inv_e_metric_;
+
+      /**
+       * Is inv_e_metric_ diagonal or not
+       */
+      bool is_diagonal_;
+
+      /**
+       * Construct a auto point in n-dimensional phase space
+       * with identity matrix as inverse mass matrix.
+       *
+       * @param n number of dimensions
+       */
+      explicit auto_e_point(int n)
+        : ps_point(n), inv_e_metric_(n, n), is_diagonal_(true) {
+        inv_e_metric_.setIdentity();
+      }
+
+      /**
+       * Set elements of mass matrix
+       *
+       * @param inv_e_metric initial mass matrix
+       */
+      void
+      set_metric(const Eigen::MatrixXd& inv_e_metric) {
+        inv_e_metric_ = inv_e_metric;
+	is_diagonal_ = false;
+      }
+
+      /**
+       * Write elements of mass matrix to string and handoff to writer.
+       *
+       * @param writer Stan writer callback
+       */
+      inline
+      void
+      write_metric(stan::callbacks::writer& writer) {
+        writer("Elements of inverse mass matrix:");
+        for (int i = 0; i < inv_e_metric_.rows(); ++i) {
+          std::stringstream inv_e_metric_ss;
+          inv_e_metric_ss << inv_e_metric_(i, 0);
+          for (int j = 1; j < inv_e_metric_.cols(); ++j)
+            inv_e_metric_ss << ", " << inv_e_metric_(i, j);
+          writer(inv_e_metric_ss.str());
+        }
+      }
+    };
+
+  }  // mcmc
+}  // stan
+
+#endif

src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp

@@ -2,7 +2,7 @@
 #define STAN_MCMC_HMC_HAMILTONIANS_DENSE_E_METRIC_HPP
 
 #include <stan/callbacks/logger.hpp>
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/fun.hpp>
 #include <stan/mcmc/hmc/hamiltonians/base_hamiltonian.hpp>
 #include <stan/mcmc/hmc/hamiltonians/dense_e_point.hpp>
 #include <boost/random/variate_generator.hpp>

src/stan/mcmc/hmc/nuts/adapt_auto_e_nuts.hpp

@@ -0,0 +1,79 @@
+#ifndef STAN_MCMC_HMC_NUTS_ADAPT_AUTO_E_NUTS_HPP
+#define STAN_MCMC_HMC_NUTS_ADAPT_AUTO_E_NUTS_HPP
+
+#include <stan/callbacks/logger.hpp>
+#include <stan/mcmc/hmc/nuts/auto_e_nuts.hpp>
+#include <stan/mcmc/hmc/mpi_cross_chain_adapter.hpp>
+#include <stan/mcmc/mpi_auto_adaptation.hpp>
+
+namespace stan {
+namespace mcmc {
+/**
+ * The No-U-Turn sampler (NUTS) with multinomial sampling
+ * with a Gaussian-Euclidean disintegration and adaptive
+ * dense or diagonal metric and adaptive step size
+ */
+template <class Model, class BaseRNG>
+class adapt_auto_e_nuts : public auto_e_nuts<Model, BaseRNG>,
+			  public mpi_cross_chain_adapter<adapt_auto_e_nuts<Model, BaseRNG>>,
+			  public stepsize_covar_adapter {
+protected:
+  const Model& model_;
+public:
+  adapt_auto_e_nuts(const Model& model, BaseRNG& rng)
+    : model_(model), auto_e_nuts<Model, BaseRNG>(model, rng),
+    stepsize_covar_adapter(model.num_params_r()) {}
+
+  ~adapt_auto_e_nuts() {}
+
+  sample
+  transition(sample& init_sample, callbacks::logger& logger) {
+    sample s = auto_e_nuts<Model, BaseRNG>::transition(init_sample,
+						       logger);
+
+    if (this->adapt_flag_) {
+      this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_,
+						s.accept_stat());
+
+      if (this -> use_cross_chain_adapt()) {
+        this -> add_cross_chain_sample(s.log_prob());
+        bool update = this -> cross_chain_adaptation(logger);
+        if (this -> is_cross_chain_adapted()) {
+          update = false;
+        }
+
+        if (update) {
+	  this->z_.is_diagonal_ = reinterpret_cast<mpi_auto_adaptation<Model> *>(this->metric_adapt)->is_diagonal_;
+
+          this->init_stepsize(logger);
+
+          this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
+          this->stepsize_adaptation_.restart();
+
+          this->set_cross_chain_stepsize();          
+        }
+      } else {
+        bool update = this->covar_adaptation_.learn_covariance(this->z_.inv_e_metric_,
+							       this->z_.q);
+        if (update) {
+          this->init_stepsize(logger);
+
+          this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
+          this->stepsize_adaptation_.restart();
+        }
+      }
+    }
+    return s;
+  }
+
+  void disengage_adaptation() {
+    base_adapter::disengage_adaptation();
+    if (!this -> is_cross_chain_adapted()) {
+      this->stepsize_adaptation_.complete_adaptation(this->nom_epsilon_);
+    }
+  }
+};
+
+}  // mcmc
+}  // stan
+#endif

src/stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp

@@ -30,7 +30,6 @@ class adapt_diag_e_nuts : public diag_e_nuts<Model, BaseRNG>,
     if (this->adapt_flag_) {
       this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_,
                                                 s.accept_stat());
-
 
       if (this -> use_cross_chain_adapt()) {
         this -> add_cross_chain_sample(s.log_prob());

src/stan/mcmc/hmc/nuts/auto_e_nuts.hpp

@@ -0,0 +1,26 @@
+#ifndef STAN_MCMC_HMC_NUTS_AUTO_E_NUTS_HPP
+#define STAN_MCMC_HMC_NUTS_AUTO_E_NUTS_HPP
+
+#include <stan/mcmc/hmc/nuts/base_nuts.hpp>
+#include <stan/mcmc/hmc/hamiltonians/dense_e_point.hpp>
+#include <stan/mcmc/hmc/hamiltonians/auto_e_metric.hpp>
+#include <stan/mcmc/hmc/integrators/expl_leapfrog.hpp>
+
+namespace stan {
+  namespace mcmc {
+    /**
+     * The No-U-Turn sampler (NUTS) with multinomial sampling
+     * with a Gaussian-Euclidean disintegration and dense metric
+     */
+    template <class Model, class BaseRNG>
+    class auto_e_nuts : public base_nuts<Model, auto_e_metric,
+					 expl_leapfrog, BaseRNG> {
+    public:
+      auto_e_nuts(const Model& model, BaseRNG& rng)
+        : base_nuts<Model, auto_e_metric, expl_leapfrog,
+                    BaseRNG>(model, rng) { }
+    };
+
+  }  // mcmc
+}  // stan
+#endif

src/stan/mcmc/hmc/nuts/base_nuts.hpp

@@ -3,7 +3,7 @@
 
 #include <stan/callbacks/logger.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/fun.hpp>
 #include <stan/mcmc/hmc/base_hmc.hpp>
 #include <stan/mcmc/hmc/hamiltonians/ps_point.hpp>
 #include <algorithm>