diff --git a/stan/math/opencl/kernels/poisson_log_glm_lpmf.hpp b/stan/math/opencl/kernels/poisson_log_glm_lpmf.hpp index 463b4510b2b..3aa5f8ee17a 100644 --- a/stan/math/opencl/kernels/poisson_log_glm_lpmf.hpp +++ b/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 + 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}}); diff --git a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp index e9be4b3b45f..c5b9ef6d455 100644 --- a/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp +++ b/stan/math/opencl/prim/poisson_log_glm_lpmf.hpp @@ -90,27 +90,23 @@ return_type_t poisson_log_glm_lpmf( matrix_cl theta_derivative_cl(N, 1); matrix_cl theta_derivative_sum_cl(wgs, 1); - const bool need_logp1 = include_summand::value; - const bool need_logp2 = include_summand::value; - matrix_cl logp_cl((need_logp1 || need_logp2) ? wgs : 0, 1); + const bool need_logp = include_summand::value; + matrix_cl 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 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)); diff --git a/stan/math/prim/prob/poisson_log_glm_lpmf.hpp b/stan/math/prim/prob/poisson_log_glm_lpmf.hpp index 0dda01658e8..d1e1586d960 100644 --- a/stan/math/prim/prob/poisson_log_glm_lpmf.hpp +++ b/stan/math/prim/prob/poisson_log_glm_lpmf.hpp @@ -114,10 +114,9 @@ return_type_t poisson_log_glm_lpmf( logp -= lgamma(forward_as(y_val) + 1); } } - if (include_summand::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, T_alpha, T_beta> diff --git a/test/unit/math/opencl/rev/poisson_log_glm_lpmf_test.cpp b/test/unit/math/opencl/rev/poisson_log_glm_lpmf_test.cpp index 2f24f72c282..20fca97a98e 100644 --- a/test/unit/math/opencl/rev/poisson_log_glm_lpmf_test.cpp +++ b/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 y{2, 0, 1, 2, 1, 0, 0, 1, 3, 0}; + Matrix 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 beta(3, 1); + beta << 1.17711, 1.24432, -0.596639; + var alpha = -1.04272; + var lp = stan::math::poisson_log_glm_lpmf( + stan::math::to_matrix_cl(y), stan::math::to_matrix_cl(x), alpha, beta); + var lp1 = stan::math::poisson_log_glm_lpmf(y, x, alpha, beta); + EXPECT_FLOAT_EQ(lp.val(), lp1.val()); +} + #endif diff --git a/test/unit/math/rev/prob/poisson_log_glm_lpmf_test.cpp b/test/unit/math/rev/prob/poisson_log_glm_lpmf_test.cpp index 35edba48a43..181ca1ee3f6 100644 --- a/test/unit/math/rev/prob/poisson_log_glm_lpmf_test.cpp +++ b/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 y{2, 0, 1, 2, 1, 0, 0, 1, 3, 0}; + Matrix 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 beta(3, 1); + beta << 1.17711, 1.24432, -0.596639; + var alpha = -1.04272; + Matrix alphavec = alpha * Matrix::Ones(); + Matrix theta(10, 1); + theta = x * beta + alphavec; + var lp = stan::math::poisson_log_lpmf(y, theta); + double lp_val = lp.val(); + var lp1 = stan::math::poisson_log_glm_lpmf(y, x, alpha, beta); + double lp1_val = lp1.val(); + EXPECT_FLOAT_EQ(lp_val, lp1_val); +}