get vector working for prim and reverse mode

Author: SteveBronder

stan/math/prim/prob/normal_log.hpp

@@ -29,7 +29,7 @@ namespace math {
  * @tparam T_loc Type of location parameter.
  */
 template <bool propto, typename T_y, typename T_loc, typename T_scale>
-inline return_type_t<T_y, T_loc, T_scale> normal_log(const T_y& y,
+inline auto normal_log(const T_y& y,
                                                      const T_loc& mu,
                                                      const T_scale& sigma) {
   return normal_lpdf<propto>(y, mu, sigma);
@@ -39,7 +39,7 @@ inline return_type_t<T_y, T_loc, T_scale> normal_log(const T_y& y,
  * @deprecated use <code>normal_lpdf</code>
  */
 template <typename T_y, typename T_loc, typename T_scale>
-inline return_type_t<T_y, T_loc, T_scale> normal_log(const T_y& y,
+inline auto normal_log(const T_y& y,
                                                      const T_loc& mu,
                                                      const T_scale& sigma) {
   return normal_lpdf<false>(y, mu, sigma);

stan/math/prim/prob/normal_lpdf.hpp

@@ -69,11 +69,12 @@ struct prob_broadcaster<T, require_stan_scalar_t<T>> {
     return *this;
   }
   inline auto ret() noexcept {
+    static_assert(!is_var<T>::value, "NOOO");
     return ret_;
   }
-  template <typename T1>
-  static auto zero(T1&& /* */) {
-    return T(0);
+  template <typename... Types>
+  static auto zero(int /* */) {
+    return return_type_t<Types...>(0);
   }
 
 };
@@ -86,7 +87,7 @@ struct prob_broadcaster<T, require_eigen_t<T>> {
 
   template <typename EigArr, require_eigen_t<EigArr>* = nullptr>
   inline auto operator=(EigArr&& x) {
-    ret_ = sum(x);
+    ret_ = x;
     return *this;
   }
 
@@ -124,9 +125,9 @@ struct prob_broadcaster<T, require_eigen_t<T>> {
     return std::move(ret_);
   }
 
-  template <typename T1>
-  static auto zero(T1&& size) {
-    return Eigen::Array<value_type_t<T>, -1, 1>::Constant(0, size).eval();
+  template <typename... Types>
+  static auto zero(int size) {
+    return Eigen::Array<return_type_t<Types...>, -1, 1>::Constant(0, size).eval();
   }
 
 };
@@ -182,10 +183,10 @@ inline auto normal_lpdf(const T_y& y,
   using ret_t = prob_return_t<RetType, T_partials_return>;
   const size_t N = max_size(y, mu, sigma);
   if (size_zero(y, mu, sigma)) {
-    return ret_t::zero(N);
+    return ret_t::template zero<T_y, T_loc, T_scale>(N);
   }
   if (!include_summand<propto, T_y, T_loc, T_scale>::value) {
-    return ret_t::zero(N);
+    return ret_t::template zero<T_y, T_loc, T_scale>(N);
   }
 
   operands_and_partials<T_y_ref, T_mu_ref, T_sigma_ref> ops_partials(

stan/math/rev/functor/operands_and_partials.hpp

@@ -87,7 +87,7 @@ template <typename Scalar1, typename Scalar2, typename T3, require_var_t<Scalar1
 inline void update_adjoints(Scalar1 x, Scalar2 y, const T3& z) noexcept {
   x.adj() += sum(z.adj() * y);
 }
-template <typename Matrix1, typename Matrix2, typename T3
+template <typename Matrix1, typename Matrix2, typename T3,
           require_rev_matrix_t<Matrix1>* = nullptr,
           require_st_arithmetic<Matrix2>* = nullptr,
           require_eigen_t<T3>* = nullptr>
@@ -214,19 +214,19 @@ class operands_and_partials<Op1, Op2, Op3, Op4, Op5, var> {
             operand5 = edge5_.operand(),
             partial5 = edge5_.partial()]() mutable {
       if (!is_constant<Op1>::value) {
-        internal::update_adjoints(operand1, partial1, vi);
+        internal::update_adjoints(operand1, partial1, ret);
       }
       if (!is_constant<Op2>::value) {
-        internal::update_adjoints(operand2, partial2, vi);
+        internal::update_adjoints(operand2, partial2, ret);
       }
       if (!is_constant<Op3>::value) {
-        internal::update_adjoints(operand3, partial3, vi);
+        internal::update_adjoints(operand3, partial3, ret);
       }
       if (!is_constant<Op4>::value) {
-        internal::update_adjoints(operand4, partial4, vi);
+        internal::update_adjoints(operand4, partial4, ret);
       }
       if (!is_constant<Op5>::value) {
-        internal::update_adjoints(operand5, partial5, vi);
+        internal::update_adjoints(operand5, partial5, ret);
       }
     });
     return plain_type_t<promote_scalar_t<var, T>>(ret);

test/unit/math/rev/prob/normal_log_test.cpp

@@ -21,3 +21,13 @@ TEST(ProbDistributionsNormal, intVsDouble) {
     EXPECT_FLOAT_EQ(lp1adj, lp2adj);
   }
 }
+
+
+TEST(ProbNormal, test_vlpdf) {
+  using stan::math::var;
+  Eigen::Matrix<var, -1, 1> Y = Eigen::Matrix<double, -1, 1>::Random(5);
+  Eigen::Matrix<var, -1, 1> Mu = Eigen::Matrix<double, -1, 1>::Random(5);
+  Eigen::Matrix<var, -1, 1> Sigma = stan::math::abs(Eigen::Matrix<double, -1, 1>::Random(5));
+  Eigen::Matrix<var, -1, 1> A = stan::math::normal_lpdf<false, stan::math::ProbReturnType::Vector>(Y, Mu, Sigma);
+
+}