Bugfix/poisson glm

stan/math/opencl/kernels/poisson_log_glm_lpmf.hpp

@@ -32,9 +32,7 @@ static const char* poisson_log_glm_kernel_code = STRINGIFY(
      * it is a scalar)
      * @param is_alpha_vector 0 or 1 - whether alpha is a vector (alternatively
      * it is a scalar)
-     * @param need_logp1 interpreted as boolean - whether first part of
-     * logp_global needs to be computed
-     * @param need_logp2 interpreted as boolean - whether second part of
+     * @param need_logp interpreted as boolean - whether first part of
      * logp_global needs to be computed
      */
     __kernel void poisson_log_glm(
@@ -43,7 +41,7 @@ static const char* poisson_log_glm_kernel_code = STRINGIFY(
         const __global int* y_global, const __global double* x,
         const __global double* alpha, const __global double* beta, const int N,
         const int M, const int is_y_vector, const int is_alpha_vector,
-        const int need_logp1, const int need_logp2) {
+        const int need_logp) {
       const int gid = get_global_id(0);
       const int lid = get_local_id(0);
       const int lsize = get_local_size(0);
@@ -68,12 +66,10 @@ static const char* poisson_log_glm_kernel_code = STRINGIFY(
           // this signals that an exception must be raised
           theta_derivative = NAN;
         }
-        if (need_logp1) {
+        if (need_logp) {
           logp = -lgamma(y + 1);
         }
-        if (need_logp2) {
-          logp += y * theta - exp_theta;
-        }
+        logp += y * theta - exp_theta;
         theta_derivative_global[gid] = theta_derivative;
       }
       // Sum theta_derivative, calculated by different threads.
@@ -94,23 +90,21 @@ static const char* poisson_log_glm_kernel_code = STRINGIFY(
         theta_derivative_sum[wg_id] = local_storage[0];
       }
 
-      if (need_logp1 || need_logp2) {
-        // Sum logp, calculated by different threads.
-        barrier(CLK_LOCAL_MEM_FENCE);
-        local_storage[lid] = logp;
-        barrier(CLK_LOCAL_MEM_FENCE);
-        for (int step = lsize / REDUCTION_STEP_SIZE; step > 0;
-             step /= REDUCTION_STEP_SIZE) {
-          if (lid < step) {
-            for (int i = 1; i < REDUCTION_STEP_SIZE; i++) {
-              local_storage[lid] += local_storage[lid + step * i];
-            }
+      // Sum logp, calculated by different threads.
+      barrier(CLK_LOCAL_MEM_FENCE);
+      local_storage[lid] = logp;
+      barrier(CLK_LOCAL_MEM_FENCE);
+      for (int step = lsize / REDUCTION_STEP_SIZE; step > 0;
+           step /= REDUCTION_STEP_SIZE) {
+        if (lid < step) {
+          for (int i = 1; i < REDUCTION_STEP_SIZE; i++) {
+            local_storage[lid] += local_storage[lid + step * i];
           }
-          barrier(CLK_LOCAL_MEM_FENCE);
-        }
-        if (lid == 0) {
-          logp_global[wg_id] = local_storage[0];
         }
+        barrier(CLK_LOCAL_MEM_FENCE);
+      }
+      if (lid == 0) {
+        logp_global[wg_id] = local_storage[0];
       }
     }
     // \cond
@@ -122,7 +116,7 @@ static const char* poisson_log_glm_kernel_code = STRINGIFY(
  * poisson_log_glm_lpmf() \endlink
  */
 const kernel_cl<out_buffer, out_buffer, out_buffer, in_buffer, in_buffer,
-                in_buffer, in_buffer, int, int, int, int, int, int>
+                in_buffer, in_buffer, int, int, int, int, int>
     poisson_log_glm("poisson_log_glm", {poisson_log_glm_kernel_code},
                     {{"REDUCTION_STEP_SIZE", 4}, {"LOCAL_SIZE_", 64}});
 

stan/math/opencl/prim/poisson_log_glm_lpmf.hpp

@@ -90,27 +90,23 @@ return_type_t<T_alpha, T_beta> poisson_log_glm_lpmf(
 
   matrix_cl<double> theta_derivative_cl(N, 1);
   matrix_cl<double> theta_derivative_sum_cl(wgs, 1);
-  const bool need_logp1 = include_summand<propto>::value;
-  const bool need_logp2 = include_summand<propto, T_partials_return>::value;
-  matrix_cl<double> logp_cl((need_logp1 || need_logp2) ? wgs : 0, 1);
+  const bool need_logp = include_summand<propto>::value;
+  matrix_cl<double> logp_cl(wgs, 1);
 
   try {
     opencl_kernels::poisson_log_glm(
         cl::NDRange(local_size * wgs), cl::NDRange(local_size),
         theta_derivative_cl, theta_derivative_sum_cl, logp_cl, y_cl, x_cl,
-        alpha_cl, beta_cl, N, M, y_cl.size() != 1, size(alpha) != 1, need_logp1,
-        need_logp2);
+        alpha_cl, beta_cl, N, M, y_cl.size() != 1, size(alpha) != 1, need_logp);
   } catch (const cl::Error& e) {
     check_opencl_error(function, e);
   }
   Matrix<T_partials_return, Dynamic, 1> theta_derivative_partial_sum(wgs);
   theta_derivative_partial_sum = from_matrix_cl(theta_derivative_sum_cl);
   double theta_derivative_sum = sum(theta_derivative_partial_sum);
-  if (need_logp1 || need_logp2) {
-    Eigen::VectorXd logp_partial_sum(wgs);
-    logp_partial_sum = from_matrix_cl(logp_cl);
-    logp += sum(logp_partial_sum);
-  }
+  Eigen::VectorXd logp_partial_sum(wgs);
+  logp_partial_sum = from_matrix_cl(logp_cl);
+  logp += sum(logp_partial_sum);
   if (!std::isfinite(theta_derivative_sum)) {
     check_nonnegative(function, "Vector of dependent variables",
                       from_matrix_cl(y_cl));

stan/math/prim/prob/poisson_log_glm_lpmf.hpp

@@ -114,10 +114,9 @@ return_type_t<T_x_scalar, T_alpha, T_beta> poisson_log_glm_lpmf(
       logp -= lgamma(forward_as<double>(y_val) + 1);
     }
   }
-  if (include_summand<propto, T_partials_return>::value) {
-    logp += sum(as_array_or_scalar(y_val_vec) * theta.array()
-                - exp(theta.array()));
-  }
+
+  logp += sum(as_array_or_scalar(y_val_vec) * theta.array()
+              - exp(theta.array()));
 
   operands_and_partials<Eigen::Matrix<T_x_scalar, T_x_rows, Eigen::Dynamic>,
                         T_alpha, T_beta>

test/unit/math/opencl/rev/poisson_log_glm_lpmf_test.cpp

@@ -330,4 +330,20 @@ TEST(ProbDistributionsPoissonLogGLM, gpu_matches_cpu_big) {
                   alpha_var2.adj().eval());
 }
 
+TEST(ProbDistributionsPoissonLogGLM, gpu_matches_cpu_poisson_log_vars_propto) {
+  vector<int> y{2, 0, 1, 2, 1, 0, 0, 1, 3, 0};
+  Matrix<double, Dynamic, Dynamic> x(10, 3);
+  x << -1.87936, 0.55093, -2.50689, 4.78584, 0.988523, -2.46141, 1.46229,
+      2.21497, 1.72734, -0.916165, -0.563808, 0.165279, -0.752066, 1.43575,
+      0.296557, -0.738422, 0.438686, 0.664492, 0.518203, -0.27485, 2, 0, 1, 2,
+      1, 0, 0, 1, 3, 0;
+  Matrix<var, Dynamic, 1> beta(3, 1);
+  beta << 1.17711, 1.24432, -0.596639;
+  var alpha = -1.04272;
+  var lp = stan::math::poisson_log_glm_lpmf<true>(
+      stan::math::to_matrix_cl(y), stan::math::to_matrix_cl(x), alpha, beta);
+  var lp1 = stan::math::poisson_log_glm_lpmf<true>(y, x, alpha, beta);
+  EXPECT_FLOAT_EQ(lp.val(), lp1.val());
+}
+
 #endif

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

@@ -470,3 +470,23 @@ TEST(ProbDistributionsPoissonLogGLM, glm_matches_poisson_log_error_checking) {
   EXPECT_THROW(stan::math::poisson_log_glm_lpmf(y, x, alpha, betaw2),
                std::domain_error);
 }
+
+TEST(ProbDistributionsPoissonLogGLM, glm_matches_poisson_log_vars_propto) {
+  vector<int> y{2, 0, 1, 2, 1, 0, 0, 1, 3, 0};
+  Matrix<double, Dynamic, Dynamic> x(10, 3);
+  x << -1.87936, 0.55093, -2.50689, 4.78584, 0.988523, -2.46141, 1.46229,
+      2.21497, 1.72734, -0.916165, -0.563808, 0.165279, -0.752066, 1.43575,
+      0.296557, -0.738422, 0.438686, 0.664492, 0.518203, -0.27485, 2, 0, 1, 2,
+      1, 0, 0, 1, 3, 0;
+  Matrix<var, Dynamic, 1> beta(3, 1);
+  beta << 1.17711, 1.24432, -0.596639;
+  var alpha = -1.04272;
+  Matrix<var, Dynamic, 1> alphavec = alpha * Matrix<var, 10, 1>::Ones();
+  Matrix<var, Dynamic, 1> theta(10, 1);
+  theta = x * beta + alphavec;
+  var lp = stan::math::poisson_log_lpmf<true>(y, theta);
+  double lp_val = lp.val();
+  var lp1 = stan::math::poisson_log_glm_lpmf<true>(y, x, alpha, beta);
+  double lp1_val = lp1.val();
+  EXPECT_FLOAT_EQ(lp_val, lp1_val);
+}