Merge branch 'develop' of https://github.com/stan-dev/stan into develop

.gitignore

@@ -41,10 +41,13 @@ src/test/gm/model_specs/printer.cpp
 src/test/mcmc/models/improper
 src/test/mcmc/models/improper.exe
 src/test/mcmc/models/improper.cpp
+src/test/mcmc/hmc/hamiltonians/models/funnel.cpp
+src/test/mcmc/hmc/integrators/models/gauss.cpp
 src/test/gm/arguments/test_model
 src/test/gm/arguments/test_model.exe
 src/test/gm/arguments/test_model.cpp
 
+
 src/test/mcmc/models/*
 !src/test/mcmc/models/*.stan
 

make/tests

@@ -156,6 +156,22 @@ ifeq ($(findstring test-,$(MAKECMDGOALS))$(findstring runtest/,$(MAKECMDGOALS)),
 	$@ --gtest_output="xml:$(basename $@).xml"
 endif
 
+src/test/mcmc/hmc/hamiltonians/models/%.cpp : src/test/mcmc/hmc/hamiltonians/models/%.stan bin/stanc$(EXE)
+	bin$(PATH_SEPARATOR)stanc$(EXE) --name=$(basename $(@F)) --o=$@ --no_main $<
+
+src/test/mcmc/hmc/integrators/models/%.cpp : src/test/mcmc/hmc/integrators/models/%.stan bin/stanc$(EXE)
+	bin$(PATH_SEPARATOR)stanc$(EXE) --name=$(basename $(@F)) --o=$@ --no_main $<
+
+.PRECIOUS: src/test/mcmc/hmc/hamiltonians/models/funnel.cpp
+test/mcmc/hmc/hamiltonians/unit_e_metric.o test/mcmc/hmc/hamiltonians/diag_e_metric.o test/mcmc/hmc/hamiltonians/dense_e_metric.o: test/%.o : src/test/%_test.cpp $(LIBGTEST) src/test/mcmc/hmc/hamiltonians/models/funnel.cpp
+	@mkdir -p $(dir $@)
+	$(COMPILE.c) -O$O $(CFLAGS_GTEST) $< $(OUTPUT_OPTION)
+
+.PRECIOUS: src/test/mcmc/hmc/integrators/models/gauss.cpp
+test/mcmc/hmc/integrators/expl_leapfrog.o: test/%.o : src/test/%_test.cpp $(LIBGTEST) src/test/mcmc/hmc/integrators/models/gauss.cpp
+	@mkdir -p $(dir $@)
+	$(COMPILE.c) -O$O $(CFLAGS_GTEST) $< $(OUTPUT_OPTION)
+
 test/gm/arguments/argument_configuration$(EXE) : test/gm/arguments/argument_configuration.o bin/libstan.a bin/libstanc.a src/test/gm/arguments/test_model$(EXE)
 	@mkdir -p $(dir $@)
 	$(LINK.c) -O$O $(GTEST_MAIN) $< $(CFLAGS_GTEST) $(OUTPUT_OPTION) $(LIBGTEST) $(LDLIBS) $(LDLIBS_STANC)
@@ -222,6 +238,7 @@ endif
 # tests.
 ## 
 .PHONY: runtest/%
+.PRECIOUS: test/%$(EXE)
 .SECONDARY:
 runtest/% : %$(EXE)
 	@echo '------------------------------------------------------------'

src/stan/agrad/partials_vari.hpp

@@ -79,12 +79,12 @@ namespace stan {
       agrad::vari** all_varis;
       double* all_partials;
 
-      VectorView<double*, is_vector<T1>::value> d_x1;
-      VectorView<double*, is_vector<T2>::value> d_x2;
-      VectorView<double*, is_vector<T3>::value> d_x3;
-      VectorView<double*, is_vector<T4>::value> d_x4;
-      VectorView<double*, is_vector<T5>::value> d_x5;
-      VectorView<double*, is_vector<T6>::value> d_x6;
+      VectorView<double*, is_vector<T1>::value, is_constant_struct<T1>::value> d_x1;
+      VectorView<double*, is_vector<T2>::value, is_constant_struct<T2>::value> d_x2;
+      VectorView<double*, is_vector<T3>::value, is_constant_struct<T3>::value> d_x3;
+      VectorView<double*, is_vector<T4>::value, is_constant_struct<T4>::value> d_x4;
+      VectorView<double*, is_vector<T5>::value, is_constant_struct<T5>::value> d_x5;
+      VectorView<double*, is_vector<T6>::value, is_constant_struct<T6>::value> d_x6;
 
       OperandsAndPartials(const T1& x1=0, const T2& x2=0, const T3& x3=0, 
                           const T4& x4=0, const T5& x5=0, const T6& x6=0)

src/stan/gm/command.hpp

@@ -371,18 +371,18 @@ namespace stan {
                                                        cont_params, disc_params, init_grad,
                                                        &std::cout);
           } catch (std::domain_error e) {
-            std::cout << "Rejecting inititialization at zero because of log_prob_grad failure." << std::endl;
+            std::cout << "Rejecting initialization at zero because of log_prob_grad failure." << std::endl;
             return error_codes::OK;
           }
 
           if (!boost::math::isfinite(init_log_prob)) {
-            std::cout << "Rejecting inititialization at zero because of vanishing density." << std::endl;
+            std::cout << "Rejecting initialization at zero because of vanishing density." << std::endl;
             return 0;
           }
 
           for (size_t i = 0; i < init_grad.size(); ++i) {
             if (!boost::math::isfinite(init_grad[i])) {
-              std::cout << "Rejecting inititialization at zero because of divergent gradient." << std::endl;
+              std::cout << "Rejecting initialization at zero because of divergent gradient." << std::endl;
               return 0;
             }
           }
@@ -477,18 +477,18 @@ namespace stan {
                                                      &std::cout);
 
         } catch (std::domain_error e) {
-          std::cout << "Rejecting user-specified inititialization because of log_prob_grad failure." << std::endl;
+          std::cout << "Rejecting user-specified initialization because of log_prob_grad failure." << std::endl;
           return 0;
         }
 
         if (!boost::math::isfinite(init_log_prob)) {
-          std::cout << "Rejecting user-specified inititialization because of vanishing density." << std::endl;
+          std::cout << "Rejecting user-specified initialization because of vanishing density." << std::endl;
           return 0;
         }
 
         for (size_t i = 0; i < init_grad.size(); ++i) {
           if (!boost::math::isfinite(init_grad[i])) {
-            std::cout << "Rejecting user-specified inititialization because of divergent gradient." << std::endl;
+            std::cout << "Rejecting user-specified initialization because of divergent gradient." << std::endl;
             return 0;
           }
         }

src/stan/io/dump.hpp

@@ -435,6 +435,7 @@ namespace stan {
 
       bool scan_single_char(char c_expected) {
         int c = in_.peek();
+        if (in_.fail()) return false;
         if (c != c_expected)
           return false;
         char c_skip;
@@ -454,6 +455,7 @@ namespace stan {
       bool scan_char(char c_expected) {
         char c;
         in_ >> c;
+        if (in_.fail()) return false;
         if (c != c_expected) {
           in_.putback(c);
           return false;
@@ -464,6 +466,7 @@ namespace stan {
       bool scan_name_unquoted() {
         char c;
         in_ >> c; // 
+        if (in_.fail()) return false;
         if (!std::isalpha(c)) return false;
         name_.push_back(c); 
         while (in_.get(c)) { // get turns off auto space skip

src/stan/io/writer.hpp

@@ -130,9 +130,7 @@ namespace stan {
        * @throw std::runtime_error if y is lower than the lower bound provided.
        */
       void scalar_lb_unconstrain(double lb, T& y) {
-        if (y < lb)
-          BOOST_THROW_EXCEPTION(std::runtime_error ("y is lower than the lower bound"));
-        data_r_.push_back(log(y - lb));
+        data_r_.push_back(stan::prob::lb_free(y,lb));
       }
 
       /**
@@ -146,9 +144,7 @@ namespace stan {
        * @throw std::runtime_error if y is higher than the upper bound provided.
        */
       void scalar_ub_unconstrain(double ub, T& y) {
-        if (y > ub)
-          BOOST_THROW_EXCEPTION(std::runtime_error ("y is higher than the lower bound"));
-        data_r_.push_back(log(ub - y));
+        data_r_.push_back(stan::prob::ub_free(y,ub));
       }
 
       /**
@@ -164,9 +160,7 @@ namespace stan {
        * @throw std::runtime_error if y is not between the lower and upper bounds
        */
       void scalar_lub_unconstrain(double lb, double ub, T& y) {
-        if (y < lb || y > ub)
-          BOOST_THROW_EXCEPTION(std::runtime_error ("y is not between the lower and upper bounds"));
-        data_r_.push_back(stan::math::logit((y - lb) / (ub - lb)));
+        data_r_.push_back(stan::prob::lub_free(y,lb,ub));
       }
 
       /**
@@ -180,9 +174,7 @@ namespace stan {
        * @throw std::runtime_error if y is not between -1.0 and 1.0
        */
       void corr_unconstrain(T& y) {
-        if (y > 1.0 || y < -1.0)
-          BOOST_THROW_EXCEPTION(std::runtime_error ("y is not between -1.0 and 1.0"));
-        data_r_.push_back(atanh(y));
+        data_r_.push_back(stan::prob::corr_free(y));
       }
 
       /**
@@ -197,9 +189,7 @@ namespace stan {
        * @throw std::runtime_error if y is not between 0.0 and 1.0
         */
       void prob_unconstrain(T& y) {
-        if (y > 1.0 || y < 0.0)
-          BOOST_THROW_EXCEPTION(std::runtime_error ("y is not between 0.0 and 1.0"));
-        data_r_.push_back(stan::math::logit(y));
+        data_r_.push_back(stan::prob::prob_free(y));
       }
 
       /**
@@ -242,6 +232,7 @@ namespace stan {
        * @throw std::runtime_error if vector is not in ascending order.
        */
       void positive_ordered_unconstrain(vector_t& y) {
+        // reimplements pos_ordered_free in prob to avoid malloc
         if (y.size() == 0) return;
         stan::math::check_positive_ordered("stan::io::positive_ordered_unconstrain(%1%)", y, "Vector");
         data_r_.push_back(log(y[0]));

src/stan/meta/traits.hpp

@@ -222,8 +222,34 @@ namespace stan {
   };
 
 
+  /**
+   *  VectorView is a template metaprogram that takes its argument and
+   *  allows it to be used like a vector. There are three template parameters
+   *  - T: Type of the thing to be wrapped. For example, double, var, vector<double>, etc.
+   *  - is_array: Boolean variable indicating whether the underlying type is an array.
+   *  - throw_if_accessed: Boolean variable indicating whether this instance should
+   *       not be used and should throw if operator[] is used.
+   *
+   *  For a scalar value, it broadcasts the single value when using
+   *  operator[].
+   *
+   *  For a vector, operator[] looks into the value passed in.
+   *  Note: this is not safe. It is possible to read past the size of
+   *  an array.
+   *
+   *  Uses: 
+   *    Read arguments to prob functions as vectors, even if scalars, so
+   *    they can be read by common code (and scalars automatically
+   *    broadcast up to behave like vectors) : VectorView of immutable
+   *    const array of double* (no allocation)
+   *
+   *    Build up derivatives into common storage : VectorView of
+   *    mutable shared array (no allocation because allocated on
+   *    auto-diff arena memory)
+   */
   template <typename T,
-            bool is_array = stan::is_vector_like<T>::value>
+            bool is_array = stan::is_vector_like<T>::value,
+            bool throw_if_accessed = false>
   class VectorView {
   public: 
     typedef typename scalar_type<T>::type scalar_t;
@@ -238,15 +264,23 @@ namespace stan {
     VectorView(scalar_t* x) : x_(x) { }
 
     scalar_t& operator[](int i) {
-      if (is_array) return x_[i];
-      else return x_[0];
+      if (throw_if_accessed) 
+        throw std::out_of_range("VectorView: this cannot be accessed");
+      if (is_array) 
+        return x_[i];
+      else 
+        return x_[0];
     }
   private:
     scalar_t* x_;
   };
 
-  template <typename T, bool is_array>
-  class VectorView<const T, is_array> {
+  /**
+   *
+   *  VectorView that has const correctness.
+   */
+  template <typename T, bool is_array, bool throw_if_accessed>
+  class VectorView<const T, is_array, throw_if_accessed> {
   public:
     typedef typename scalar_type<T>::type scalar_t;
 
@@ -259,17 +293,21 @@ namespace stan {
     template <int R, int C>
     VectorView(const Eigen::Matrix<scalar_t,R,C>& m) : x_(&m(0)) { }
 
-    const scalar_t operator[](int i) const {
-      if (is_array) return x_[i];
-      else return x_[0];
+    const scalar_t& operator[](int i) const {
+      if (throw_if_accessed) 
+        throw std::out_of_range("VectorView: this cannot be accessed");
+      if (is_array)
+        return x_[i];
+      else 
+        return x_[0];
     }
   private:
     const scalar_t* x_;
   };
 
   // simplify to hold value in common case where it's more efficient
   template <>
-  class VectorView<const double, false> {
+  class VectorView<const double, false, false> {
   public:
     VectorView(double x) : x_(x) { }
     double operator[](int /* i */)  const {
@@ -279,6 +317,21 @@ namespace stan {
     const double x_;
   };
 
+
+  /**
+   *
+   *  DoubleVectorView allocates double values to be used as
+   *  intermediate values. There are 2 template parameters:
+   *  - used: boolean variable indicating whether this instance
+   *      is used. If this is false, there is no storage allocated
+   *      and operator[] throws.
+   *  - is_vec: boolean variable indicating whether this instance
+   *      should allocate a vector, if it is used. If this is false,
+   *      the instance will only allocate a single double value.
+   *      If this is true, it will allocate the number requested.
+   *
+   *  These values are mutable.
+   */
   template<bool used, bool is_vec>
   class DoubleVectorView {
   public:

src/stan/prob/distributions/univariate/continuous/beta.hpp

@@ -270,10 +270,6 @@ namespace stan {
       agrad::OperandsAndPartials<T_y, T_scale_succ, T_scale_fail> 
         operands_and_partials(y, alpha, beta);
 
-      std::fill(operands_and_partials.all_partials,
-                operands_and_partials.all_partials 
-                + operands_and_partials.nvaris, 0.0);
-
       // Explicit return for extreme values
       // The gradients are technically ill-defined, but treated as zero
       for (size_t i = 0; i < stan::length(y); i++) {
@@ -430,10 +426,6 @@ namespace stan {
       agrad::OperandsAndPartials<T_y, T_scale_succ, T_scale_fail> 
         operands_and_partials(y, alpha, beta);
 
-      std::fill(operands_and_partials.all_partials,
-                operands_and_partials.all_partials 
-                + operands_and_partials.nvaris, 0.0);
-
       // Compute CDF and its gradients
       using boost::math::ibeta;
       using boost::math::ibeta_derivative;
@@ -563,10 +555,6 @@ namespace stan {
       agrad::OperandsAndPartials<T_y, T_scale_succ, T_scale_fail> 
         operands_and_partials(y, alpha, beta);
 
-      std::fill(operands_and_partials.all_partials,
-                operands_and_partials.all_partials
-                + operands_and_partials.nvaris, 0.0);
-
       // Compute CDF and its gradients
       using boost::math::ibeta;
       using boost::math::ibeta_derivative;

src/stan/prob/distributions/univariate/continuous/cauchy.hpp

@@ -199,10 +199,6 @@ namespace stan {
       agrad::OperandsAndPartials<T_y, T_loc, T_scale>
         operands_and_partials(y, mu, sigma);
 
-      std::fill(operands_and_partials.all_partials,
-                operands_and_partials.all_partials
-                + operands_and_partials.nvaris, 0.0);
-
       // Explicit return for extreme values
       // The gradients are technically ill-defined, but treated as zero
       for (size_t i = 0; i < stan::length(y); i++) {
@@ -305,10 +301,6 @@ namespace stan {
       agrad::OperandsAndPartials<T_y, T_loc, T_scale> 
         operands_and_partials(y, mu, sigma);
 
-      std::fill(operands_and_partials.all_partials,
-                operands_and_partials.all_partials 
-                + operands_and_partials.nvaris, 0.0);
-
       // Compute CDFLog and its gradients
       using std::atan;
       using stan::math::pi;
@@ -382,10 +374,6 @@ namespace stan {
       agrad::OperandsAndPartials<T_y, T_loc, T_scale> 
         operands_and_partials(y, mu, sigma);
 
-      std::fill(operands_and_partials.all_partials,
-                operands_and_partials.all_partials 
-                + operands_and_partials.nvaris, 0.0);
-
       // Compute CDFLog and its gradients
       using std::atan;
       using stan::math::pi;