Add perfect forwarding and constexpr to reverse mode functions

make/compiler_flags

@@ -120,32 +120,8 @@ INC_GTEST ?= -I $(GTEST)/include -I $(GTEST)
 CPPFLAGS_BOOST ?= -DBOOST_DISABLE_ASSERTS
 CPPFLAGS_SUNDIALS ?= -DNO_FPRINTF_OUTPUT $(CPPFLAGS_OPTIM_SUNDIALS) $(CXXFLAGS_FLTO_SUNDIALS)
 #CPPFLAGS_GTEST ?=
-STAN_HAS_CXX17 ?= false
-ifeq ($(CXX_TYPE), gcc)
-  GCC_GE_73 := $(shell [ $(CXX_MAJOR) -gt 7 -o \( $(CXX_MAJOR) -eq 7 -a $(CXX_MINOR) -ge 1 \) ] && echo true)
-  ifeq ($(GCC_GE_73),true)
-    STAN_HAS_CXX17 := true
-  endif
-else ifeq ($(CXX_TYPE), clang)
-  CLANG_GE_5 := $(shell [ $(CXX_MAJOR) -gt 5 -o \( $(CXX_MAJOR) -eq 5 -a $(CXX_MINOR) -ge 0 \) ] && echo true)
-  ifeq ($(CLANG_GE_5),true)
-    STAN_HAS_CXX17 := true
-  endif
-else ifeq ($(CXX_TYPE), mingw32-gcc)
-  MINGW_GE_50 := $(shell [ $(CXX_MAJOR) -gt 5 -o \( $(CXX_MAJOR) -eq 5 -a $(CXX_MINOR) -ge 0 \) ] && echo true)
-  ifeq ($(MINGW_GE_50),true)
-    STAN_HAS_CXX17 := true
-  endif
-endif
-
-ifeq ($(STAN_HAS_CXX17), true)
-  CXXFLAGS_LANG ?= -std=c++17
-  CXXFLAGS_STANDARD ?= c++17
-else
-  $(warning "Stan cannot detect if your compiler has the C++17 standard. If it does, please set STAN_HAS_CXX17=true in your make/local file. C++17 support is mandatory in the next release of Stan. Defaulting to C++14")
-  CXXFLAGS_LANG ?= -std=c++1y
-  CXXFLAGS_STANDARD ?= c++1y
-endif
+CXXFLAGS_LANG ?= -std=c++17
+CXXFLAGS_STANDARD ?= c++17
 #CXXFLAGS_BOOST ?=
 CXXFLAGS_SUNDIALS ?= -pipe $(CXXFLAGS_OPTIM_SUNDIALS) $(CPPFLAGS_FLTO_SUNDIALS)
 #CXXFLAGS_GTEST

stan/math/prim/constraint/lb_constrain.hpp

@@ -149,11 +149,16 @@ inline auto lb_constrain(const T& x, const L& lb, return_type_t<T, L>& lp) {
  * @param[in] lb lower bound on output
  * @return lower-bound constrained value corresponding to inputs
  */
-template <typename T, typename L, require_not_std_vector_t<L>* = nullptr>
-inline auto lb_constrain(const std::vector<T>& x, const L& lb) {
+template <typename T, typename L, require_std_vector_t<T>* = nullptr,
+          require_not_std_vector_t<L>* = nullptr>
+inline auto lb_constrain(T&& x, L&& lb) {
   std::vector<plain_type_t<decltype(lb_constrain(x[0], lb))>> ret(x.size());
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lb_constrain(x[i], lb);
+    if constexpr (std::is_rvalue_reference_v<T&&>) {
+      ret[i] = lb_constrain(std::move(x[i]), lb);
+    } else {
+      ret[i] = lb_constrain(x[i], lb);
+    }
   }
   return ret;
 }
@@ -169,12 +174,16 @@ inline auto lb_constrain(const std::vector<T>& x, const L& lb) {
  * @param[in,out] lp reference to log probability to increment
  * @return lower-bound constrained value corresponding to inputs
  */
-template <typename T, typename L, require_not_std_vector_t<L>* = nullptr>
-inline auto lb_constrain(const std::vector<T>& x, const L& lb,
-                         return_type_t<T, L>& lp) {
+template <typename T, typename L, require_std_vector_t<T>* = nullptr,
+          require_not_std_vector_t<L>* = nullptr>
+inline auto lb_constrain(T&& x, L&& lb, return_type_t<T, L>& lp) {
   std::vector<plain_type_t<decltype(lb_constrain(x[0], lb))>> ret(x.size());
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lb_constrain(x[i], lb, lp);
+    if constexpr (std::is_rvalue_reference_v<T&&>) {
+      ret[i] = lb_constrain(std::move(x[i]), lb, lp);
+    } else {
+      ret[i] = lb_constrain(x[i], lb, lp);
+    }
   }
   return ret;
 }
@@ -189,12 +198,21 @@ inline auto lb_constrain(const std::vector<T>& x, const L& lb,
  * @param[in] lb lower bound on output
  * @return lower-bound constrained value corresponding to inputs
  */
-template <typename T, typename L>
-inline auto lb_constrain(const std::vector<T>& x, const std::vector<L>& lb) {
+template <typename T, typename L, require_all_std_vector_t<T, L>* = nullptr>
+inline auto lb_constrain(T&& x, L&& lb) {
   check_matching_dims("lb_constrain", "x", x, "lb", lb);
   std::vector<plain_type_t<decltype(lb_constrain(x[0], lb[0]))>> ret(x.size());
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lb_constrain(x[i], lb[i]);
+    if constexpr (std::is_rvalue_reference_v<
+                      T&&> && std::is_rvalue_reference_v<L&&>) {
+      ret[i] = lb_constrain(std::move(x[i]), std::move(lb[i]));
+    } else if constexpr (std::is_rvalue_reference_v<T&&>) {
+      ret[i] = lb_constrain(std::move(x[i]), lb[i]);
+    } else if constexpr (std::is_rvalue_reference_v<L&&>) {
+      ret[i] = lb_constrain(x[i], std::move(lb[i]));
+    } else {
+      ret[i] = lb_constrain(x[i], lb[i]);
+    }
   }
   return ret;
 }
@@ -210,13 +228,21 @@ inline auto lb_constrain(const std::vector<T>& x, const std::vector<L>& lb) {
  * @param[in,out] lp reference to log probability to increment
  * @return lower-bound constrained value corresponding to inputs
  */
-template <typename T, typename L>
-inline auto lb_constrain(const std::vector<T>& x, const std::vector<L>& lb,
-                         return_type_t<T, L>& lp) {
+template <typename T, typename L, require_all_std_vector_t<T, L>* = nullptr>
+inline auto lb_constrain(T&& x, L&& lb, return_type_t<T, L>& lp) {
   check_matching_dims("lb_constrain", "x", x, "lb", lb);
   std::vector<plain_type_t<decltype(lb_constrain(x[0], lb[0]))>> ret(x.size());
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = lb_constrain(x[i], lb[i], lp);
+    if constexpr (std::is_rvalue_reference_v<
+                      T&&> && std::is_rvalue_reference_v<L&&>) {
+      ret[i] = lb_constrain(std::move(x[i]), std::move(lb[i]), lp);
+    } else if constexpr (std::is_rvalue_reference_v<T&&>) {
+      ret[i] = lb_constrain(std::move(x[i]), lb[i], lp);
+    } else if constexpr (std::is_rvalue_reference_v<L&&>) {
+      ret[i] = lb_constrain(x[i], std::move(lb[i]), lp);
+    } else {
+      ret[i] = lb_constrain(x[i], lb[i], lp);
+    }
   }
   return ret;
 }
@@ -240,11 +266,11 @@ inline auto lb_constrain(const std::vector<T>& x, const std::vector<L>& lb,
  * @return lower-bound constrained value corresponding to inputs
  */
 template <bool Jacobian, typename T, typename L>
-inline auto lb_constrain(const T& x, const L& lb, return_type_t<T, L>& lp) {
-  if (Jacobian) {
-    return lb_constrain(x, lb, lp);
+inline auto lb_constrain(T&& x, L&& lb, return_type_t<T, L>& lp) {
+  if constexpr (Jacobian) {
+    return lb_constrain(std::forward<T>(x), std::forward<L>(lb), lp);
   } else {
-    return lb_constrain(x, lb);
+    return lb_constrain(std::forward<T>(x), std::forward<L>(lb));
   }
 }
 

stan/math/prim/constraint/ordered_constrain.hpp

@@ -51,12 +51,12 @@ inline plain_type_t<EigVec> ordered_constrain(const EigVec& x) {
  * @return Positive, increasing ordered vector.
  */
 template <typename EigVec, require_eigen_col_vector_t<EigVec>* = nullptr>
-inline auto ordered_constrain(const EigVec& x, value_type_t<EigVec>& lp) {
-  const auto& x_ref = to_ref(x);
-  if (likely(x.size() > 1)) {
-    lp += sum(x_ref.tail(x.size() - 1));
+inline auto ordered_constrain(EigVec&& x, value_type_t<EigVec>& lp) {
+  auto&& x_ref = to_ref(std::forward<EigVec>(x));
+  if (likely(x_ref.size() > 1)) {
+    lp += sum(x_ref.tail(x_ref.size() - 1));
   }
-  return ordered_constrain(x_ref);
+  return ordered_constrain(std::forward<decltype(x_ref)>(x_ref));
 }
 
 /**
@@ -78,11 +78,11 @@ inline auto ordered_constrain(const EigVec& x, value_type_t<EigVec>& lp) {
  * @return Positive, increasing ordered vector.
  */
 template <bool Jacobian, typename T, require_not_std_vector_t<T>* = nullptr>
-inline auto ordered_constrain(const T& x, return_type_t<T>& lp) {
-  if (Jacobian) {
-    return ordered_constrain(x, lp);
+inline auto ordered_constrain(T&& x, return_type_t<T>& lp) {
+  if constexpr (Jacobian) {
+    return ordered_constrain(std::forward<T>(x), lp);
   } else {
-    return ordered_constrain(x);
+    return ordered_constrain(std::forward<T>(x));
   }
 }
 
@@ -105,9 +105,10 @@ inline auto ordered_constrain(const T& x, return_type_t<T>& lp) {
  * @return Positive, increasing ordered vector.
  */
 template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr>
-inline auto ordered_constrain(const T& x, return_type_t<T>& lp) {
-  return apply_vector_unary<T>::apply(
-      x, [&lp](auto&& v) { return ordered_constrain<Jacobian>(v, lp); });
+inline auto ordered_constrain(T&& x, return_type_t<T>& lp) {
+  return apply_vector_unary<T>::apply(std::forward<T>(x), [&lp](auto&& v) {
+    return ordered_constrain<Jacobian>(std::forward<decltype(v)>(v), lp);
+  });
 }
 
 }  // namespace math

stan/math/prim/constraint/simplex_constrain.hpp

@@ -99,12 +99,11 @@ inline plain_type_t<Vec> simplex_constrain(const Vec& y,
  * @return simplex of dimensionality one greater than `y`
  */
 template <bool Jacobian, typename Vec, require_not_std_vector_t<Vec>* = nullptr>
-inline plain_type_t<Vec> simplex_constrain(const Vec& y,
-                                           return_type_t<Vec>& lp) {
-  if (Jacobian) {
-    return simplex_constrain(y, lp);
+inline auto simplex_constrain(Vec&& y, return_type_t<Vec>& lp) {
+  if constexpr (Jacobian) {
+    return simplex_constrain(std::forward<Vec>(y), lp);
   } else {
-    return simplex_constrain(y);
+    return simplex_constrain(std::forward<Vec>(y));
   }
 }
 
@@ -125,9 +124,10 @@ inline plain_type_t<Vec> simplex_constrain(const Vec& y,
  * @return simplex of dimensionality one greater than `y`
  */
 template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr>
-inline auto simplex_constrain(const T& y, return_type_t<T>& lp) {
-  return apply_vector_unary<T>::apply(
-      y, [&lp](auto&& v) { return simplex_constrain<Jacobian>(v, lp); });
+inline auto simplex_constrain(T&& y, return_type_t<T>& lp) {
+  return apply_vector_unary<T>::apply(std::forward<T>(y), [&lp](auto&& v) {
+    return simplex_constrain<Jacobian>(std::forward<decltype(v)>(v), lp);
+  });
 }
 
 }  // namespace math

stan/math/prim/constraint/ub_constrain.hpp

@@ -159,11 +159,16 @@ inline auto ub_constrain(const T& x, const U& ub,
  * @param[in] ub upper bound on output
  * @return lower-bound constrained value corresponding to inputs
  */
-template <typename T, typename U, require_not_std_vector_t<U>* = nullptr>
-inline auto ub_constrain(const std::vector<T>& x, const U& ub) {
+template <typename T, typename U, require_std_vector_t<T>* = nullptr,
+          require_not_std_vector_t<U>* = nullptr>
+inline auto ub_constrain(T&& x, const U& ub) {
   std::vector<plain_type_t<decltype(ub_constrain(x[0], ub))>> ret(x.size());
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = ub_constrain(x[i], ub);
+    if constexpr (std::is_rvalue_reference_v<T&&>) {
+      ret[i] = ub_constrain(std::move(x[i]), ub);
+    } else {
+      ret[i] = ub_constrain(x[i], ub);
+    }
   }
   return ret;
 }
@@ -179,12 +184,16 @@ inline auto ub_constrain(const std::vector<T>& x, const U& ub) {
  * @param[in,out] lp reference to log probability to increment
  * @return lower-bound constrained value corresponding to inputs
  */
-template <typename T, typename U, require_not_std_vector_t<U>* = nullptr>
-inline auto ub_constrain(const std::vector<T>& x, const U& ub,
-                         return_type_t<T, U>& lp) {
+template <typename T, typename U, require_std_vector_t<T>* = nullptr,
+          require_not_std_vector_t<U>* = nullptr>
+inline auto ub_constrain(T&& x, const U& ub, return_type_t<T, U>& lp) {
   std::vector<plain_type_t<decltype(ub_constrain(x[0], ub))>> ret(x.size());
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = ub_constrain(x[i], ub, lp);
+    if constexpr (std::is_rvalue_reference_v<T&&>) {
+      ret[i] = ub_constrain(std::move(x[i]), ub, lp);
+    } else {
+      ret[i] = ub_constrain(x[i], ub, lp);
+    }
   }
   return ret;
 }
@@ -199,12 +208,21 @@ inline auto ub_constrain(const std::vector<T>& x, const U& ub,
  * @param[in] ub upper bound on output
  * @return lower-bound constrained value corresponding to inputs
  */
-template <typename T, typename U>
-inline auto ub_constrain(const std::vector<T>& x, const std::vector<U>& ub) {
+template <typename T, typename U, require_all_std_vector_t<T, U>* = nullptr>
+inline auto ub_constrain(T&& x, U&& ub) {
   check_matching_dims("ub_constrain", "x", x, "ub", ub);
   std::vector<plain_type_t<decltype(ub_constrain(x[0], ub[0]))>> ret(x.size());
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = ub_constrain(x[i], ub[i]);
+    if constexpr (std::is_rvalue_reference_v<
+                      T&&> && std::is_rvalue_reference_v<U&&>) {
+      ret[i] = ub_constrain(std::move(x[i]), std::move(ub[i]));
+    } else if constexpr (std::is_rvalue_reference_v<T&&>) {
+      ret[i] = ub_constrain(std::move(x[i]), ub[i]);
+    } else if constexpr (std::is_rvalue_reference_v<U&&>) {
+      ret[i] = ub_constrain(x[i], std::move(ub[i]));
+    } else {
+      ret[i] = ub_constrain(x[i], ub[i]);
+    }
   }
   return ret;
 }
@@ -220,13 +238,21 @@ inline auto ub_constrain(const std::vector<T>& x, const std::vector<U>& ub) {
  * @param[in,out] lp reference to log probability to increment
  * @return lower-bound constrained value corresponding to inputs
  */
-template <typename T, typename U>
-inline auto ub_constrain(const std::vector<T>& x, const std::vector<U>& ub,
-                         return_type_t<T, U>& lp) {
+template <typename T, typename U, require_all_std_vector_t<T, U>* = nullptr>
+inline auto ub_constrain(T&& x, U&& ub, return_type_t<T, U>& lp) {
   check_matching_dims("ub_constrain", "x", x, "ub", ub);
   std::vector<plain_type_t<decltype(ub_constrain(x[0], ub[0]))>> ret(x.size());
   for (size_t i = 0; i < x.size(); ++i) {
-    ret[i] = ub_constrain(x[i], ub[i], lp);
+    if constexpr (std::is_rvalue_reference_v<
+                      T&&> && std::is_rvalue_reference_v<U&&>) {
+      ret[i] = ub_constrain(std::move(x[i]), std::move(ub[i]), lp);
+    } else if constexpr (std::is_rvalue_reference_v<T&&>) {
+      ret[i] = ub_constrain(std::move(x[i]), ub[i], lp);
+    } else if constexpr (std::is_rvalue_reference_v<U&&>) {
+      ret[i] = ub_constrain(x[i], std::move(ub[i]), lp);
+    } else {
+      ret[i] = ub_constrain(x[i], ub[i], lp);
+    }
   }
   return ret;
 }
@@ -250,11 +276,11 @@ inline auto ub_constrain(const std::vector<T>& x, const std::vector<U>& ub,
  * @return lower-bound constrained value corresponding to inputs
  */
 template <bool Jacobian, typename T, typename U>
-inline auto ub_constrain(const T& x, const U& ub, return_type_t<T, U>& lp) {
-  if (Jacobian) {
-    return ub_constrain(x, ub, lp);
+inline auto ub_constrain(T&& x, U&& ub, return_type_t<T, U>& lp) {
+  if constexpr (Jacobian) {
+    return ub_constrain(std::forward<T>(x), std::forward<U>(ub), lp);
   } else {
-    return ub_constrain(x, ub);
+    return ub_constrain(std::forward<T>(x), std::forward<U>(ub));
   }
 }
 

stan/math/prim/constraint/unit_vector_constrain.hpp

@@ -73,11 +73,11 @@ inline plain_type_t<T1> unit_vector_constrain(const T1& y, T2& lp) {
  * @return Unit length vector of dimension K
  */
 template <bool Jacobian, typename T, require_not_std_vector_t<T>* = nullptr>
-inline auto unit_vector_constrain(const T& y, return_type_t<T>& lp) {
-  if (Jacobian) {
-    return unit_vector_constrain(y, lp);
+inline auto unit_vector_constrain(T&& y, return_type_t<T>& lp) {
+  if constexpr (Jacobian) {
+    return unit_vector_constrain(std::forward<T>(y), lp);
   } else {
-    return unit_vector_constrain(y);
+    return unit_vector_constrain(std::forward<T>(y));
   }
 }
 
@@ -98,9 +98,10 @@ inline auto unit_vector_constrain(const T& y, return_type_t<T>& lp) {
  * @return Unit length vector of dimension K
  */
 template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr>
-inline auto unit_vector_constrain(const T& y, return_type_t<T>& lp) {
-  return apply_vector_unary<T>::apply(
-      y, [&lp](auto&& v) { return unit_vector_constrain<Jacobian>(v, lp); });
+inline auto unit_vector_constrain(T&& y, return_type_t<T>& lp) {
+  return apply_vector_unary<T>::apply(std::forward<T>(y), [&lp](auto&& v) {
+    return unit_vector_constrain<Jacobian>(std::forward<decltype(v)>(v), lp);
+  });
 }
 
 }  // namespace math

stan/math/prim/fun/grad_reg_inc_gamma.hpp

@@ -49,9 +49,9 @@ namespace math {
  (a-1)_k\right) \frac{1}{z^k} \end{array} \f]
  */
 template <typename T1, typename T2>
-return_type_t<T1, T2> grad_reg_inc_gamma(T1 a, T2 z, T1 g, T1 dig,
-                                         double precision = 1e-6,
-                                         int max_steps = 1e5) {
+inline return_type_t<T1, T2> grad_reg_inc_gamma(T1 a, T2 z, T1 g, T1 dig,
+                                                double precision = 1e-6,
+                                                int max_steps = 1e5) {
   using std::exp;
   using std::fabs;
   using std::log;

stan/math/prim/meta/is_autodiff.hpp

@@ -19,6 +19,14 @@ struct is_autodiff
     : bool_constant<math::disjunction<is_var<std::decay_t<T>>,
                                       is_fvar<std::decay_t<T>>>::value> {};
 
+template <typename... Types>
+inline constexpr bool is_autodiff_v
+    = math::conjunction<is_autodiff<Types>...>::value;
+
+template <typename... Types>
+inline constexpr bool is_autodiffable_v
+    = math::conjunction<is_autodiff<scalar_type_t<Types>>...>::value;
+
 /*! \ingroup require_stan_scalar_real */
 /*! \defgroup autodiff_types autodiff  */
 /*! \addtogroup autodiff_types */