Added variadic implementation of rev parallel_sum

stan/math/prim/scal.hpp

@@ -129,6 +129,7 @@
 #include <stan/math/prim/scal/fun/value_of.hpp>
 #include <stan/math/prim/scal/fun/value_of_rec.hpp>
 
+#include <stan/math/prim/scal/functor/apply.hpp>
 #include <stan/math/prim/scal/functor/reduce_sum.hpp>
 
 #include <stan/math/prim/prob/bernoulli_ccdf_log.hpp>

stan/math/prim/scal/functor/apply.hpp

@@ -0,0 +1,50 @@
+#ifndef STAN_MATH_PRIM_SCAL_FUNCTOR_APPLY_HPP
+#define STAN_MATH_PRIM_SCAL_FUNCTOR_APPLY_HPP
+
+#include <functional>
+#include <tuple>
+#include <utility>
+
+namespace stan {
+namespace math {
+/*
+ * Invoke the functor f with arguments given in t and indexed in the index
+ * sequence I
+ *
+ * @tparam F Type of functor
+ * @tparam Tuple Type of tuple containing arguments
+ * @tparam I Index sequence going from 0 to std::tuple_size<T>::value - 1
+ * inclusive
+ * @param f functor callable
+ * @param t tuple of arguments
+ * @param i placeholder variable for index sequence
+ */
+template <class F, class Tuple, std::size_t... I>
+constexpr decltype(auto) apply_impl(F&& f, Tuple&& t,
+                                    std::index_sequence<I...> i) {
+  return f(std::forward<decltype(std::get<I>(t))>(std::get<I>(t))...);
+}
+
+/*
+ * Call the functor f with the tuple of arguments t, like:
+ *
+ * f(std::get<0>(t), std::get<1>(t), ...)
+ *
+ * TODO: replace this with implementation in C++ std when C++17 is available
+ *
+ * @tparam F Type of functor
+ * @tparam Tuple Type of tuple containing arguments
+ * @param f functor callable
+ * @param t tuple of arguments
+ */
+template <class F, class Tuple>
+constexpr decltype(auto) apply(F&& f, Tuple&& t) {
+  return apply_impl(std::forward<F>(f), std::forward<Tuple>(t),
+                    std::make_index_sequence<
+                        std::tuple_size<std::remove_reference_t<Tuple>>{}>{});
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

stan/math/prim/scal/functor/reduce_sum.hpp

@@ -19,59 +19,66 @@ namespace math {
 
 namespace internal {
 
-
-// base definition => compile error
-template <typename ReduceFunction, typename Enable = void, typename M, typename T, typename... Args>
-struct reduce_sum_impl {};
-
-// todo, double check if I need enable if here
-template <typename ReduceFunction, typename M, typename T, typename... Args>
-struct reduce_sum_impl<ReduceFunction, require_arithmetic_t<T>, M, T, Args...> {
-  using vmapped_t = std::vector<M>;
-  // todo(Steve): Need helper funcs for manipulating / iterating tuples
-  std::tuple<std::vector<Args>...> arg_;
-  const vmapped_t& vmapped_;
-  T terms_sum_;
-
-  // todo(Steve): put back recursive reducer
-    reduce_sum_impl(const vmapped_t& vmapped, const T& init, Args&&... args)
-        : vmapped_(vmapped),
-          arg_(std::make_tuple(std::forward<Args>(args)...))
-          terms_sum_(value_of(init)) {}
-
-    reduce_sum_impl(reduce_sum_impl& other, tbb::split)
-        : vmapped_(other.vmapped_),
-        arg_(other.arg_),
-        terms_sum_(0.0) {}
+template <typename ReduceFunction, typename M, typename... Args>
+struct reduce_sum_impl {
+  struct recursive_reducer {
+    using vmapped_t = std::vector<M>;
+    std::tuple<const Args&...> args_tuple_;
+    const vmapped_t& vmapped_;
+    double sum_;
+
+    recursive_reducer(const vmapped_t& vmapped, const Args&... args)
+      : vmapped_(vmapped),
+	args_tuple_(args...),
+	sum_(0.0) {}
+
+    recursive_reducer(recursive_reducer& other, tbb::split)
+      : vmapped_(other.vmapped_),
+	args_tuple_(other.args_tuple_),
+	sum_(other.sum_) {}
 
     void operator()(const tbb::blocked_range<size_t>& r) {
       if (r.empty()) {
-        return;
+	return;
       }
-
+    
       auto start = vmapped_.begin();
       std::advance(start, r.begin());
       auto end = vmapped_.begin();
       std::advance(end, r.end());
-
+    
       const vmapped_t sub_slice(start, end);
-
-      terms_sum_ += ReduceFunction()(r.begin(), r.end() - 1, sub_slice, unpack_tuple_func(args));
+
+      sum_ += apply([&](auto&&... args) {
+	  return ReduceFunction()(r.begin(), r.end() - 1, sub_slice, args...);
+	}, args_tuple_);
     }
-
-    void join(const reduce_sum_impl& child) {
-      terms_sum_ += child.terms_sum_;
+  
+    void join(const recursive_reducer& child) {
+      sum_ += child.sum_;
     }
+  };
 
-  // Todo: Better name for this than operator(), collect()?
-  template <typename... OtherArgs>
-  T operator()(const vmapped_t& vmapped, T init, std::size_t grainsize,
-               Args&&... args) const {
+  double operator()(const std::vector<M>& vmapped, std::size_t grainsize,
+	       const Args&... args) const {
     const std::size_t num_jobs = vmapped.size();
-    reduce_sum_impl<ReduceFunction, M, T, OtherArgs...> worker(vmapped, init, args...);
+
+    if (num_jobs == 0)
+      return 0.0;
+
+    recursive_reducer worker(vmapped, args...);
+
+#ifdef STAN_DETERMINISTIC
+    tbb::static_partitioner partitioner;
+    tbb::parallel_deterministic_reduce(
+        tbb::blocked_range<std::size_t>(0, num_jobs, grainsize), worker,
+        partitioner);
+#else
     tbb::parallel_reduce(
         tbb::blocked_range<std::size_t>(0, num_jobs, grainsize), worker);
-    return std::move(worker.terms_sum_);
+#endif
+
+    return worker.sum_;
   }
 };
 
@@ -82,17 +89,13 @@ struct reduce_sum_impl<ReduceFunction, require_arithmetic_t<T>, M, T, Args...> {
  * that any internal state of the functor is causing trouble. Thus,
  * the functor must be default constructible without any arguments.
  */
-template <typename ReduceFunction, typename M, typename T, typename... Args>
-constexpr T reduce_sum(const std::vector<M>& vmapped, T init,
-                       std::size_t grainsize, Args&&... args) {
-  typedef T return_base_t;
-  // void here but need to figure out enable_if stuff
-  // We do this somewhere in the opencl code
-  return internal::reduce_sum_impl<ReduceFunction, M, T, void, Args...>()(
-    vmapped, init, grainsize, std::forward<Args>(args)...);
+template <typename ReduceFunction, typename M, typename... Args>
+constexpr double reduce_sum(const std::vector<M>& vmapped,
+			    std::size_t grainsize, const Args&... args) {
+  return internal::reduce_sum_impl<ReduceFunction, M, Args...>()
+    (vmapped, grainsize, args...);
 }
 
-
 }  // namespace math
 }  // namespace stan
 

stan/math/rev/functor/adj_jac_apply.hpp

@@ -13,40 +13,6 @@ namespace stan {
 namespace math {
 
 namespace internal {
-/**
- * Invoke the functor f with arguments given in t and indexed in the index
- * sequence I
- *
- * @tparam F Type of functor
- * @tparam Tuple Type of tuple containing arguments
- * @tparam I Index sequence going from 0 to std::tuple_size<T>::value - 1
- * inclusive
- * @param f functor callable
- * @param t tuple of arguments
- * @param i placeholder variable for index sequence
- */
-template <class F, class Tuple, std::size_t... I>
-constexpr auto apply_impl(const F& f, const Tuple& t,
-                          std::index_sequence<I...> i) {
-  return f(std::get<I>(t)...);
-}
-
-/**
- * Call the functor f with the tuple of arguments t, like:
- *
- * f(std::get<0>(t), std::get<1>(t), ...)
- *
- * TODO: replace this with implementation in C++ std when C++17 is available
- *
- * @tparam F Type of functor
- * @tparam Tuple Type of tuple containing arguments
- * @param f functor callable
- * @param t tuple of arguments
- */
-template <class F, class Tuple>
-constexpr auto apply(const F& f, const Tuple& t) {
-  return apply_impl(f, t, std::make_index_sequence<std::tuple_size<Tuple>{}>{});
-}
 
 /**
  * Store the adjoint in y_vi[0] in y_adj
@@ -533,7 +499,7 @@ struct adj_jac_vari : public vari {
     internal::build_y_adj(y_vi_, M_, y_adj);
     auto y_adj_jacs = f_.multiply_adjoint_jacobian(is_var_, y_adj);
 
-    internal::apply(
+    apply(
         [this](auto&&... args) { this->accumulate_adjoints(args...); },
         y_adj_jacs);
   }