Skip to content

Commit 4e1f542

Browse files
SteveBronderSteveBronder
committed
Fix normal_lccdf for opencl so that LOG_HALF is added N times to the result
1 parent 400c94d commit 4e1f542

2 files changed

Lines changed: 175 additions & 24 deletions

File tree

‎stan/math/opencl/prim/normal_lccdf.hpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -64,12 +64,12 @@ return_type_t<T_y_cl, T_loc_cl, T_scale_cl> normal_lccdf(
6464
auto sigma_positive_expr = 0 < sigma_val;
6565

6666
auto scaled_diff = elt_divide(y_val - mu_val, sigma_val * SQRT_TWO);
67-
auto one_m_erf = select(
67+
matrix_cl<double> one_m_erf = select(
6868
scaled_diff < -37.5 * INV_SQRT_TWO, 2.0,
6969
select(scaled_diff < -5.0 * INV_SQRT_TWO, 2.0 - erfc(-scaled_diff),
7070
select(scaled_diff > 8.25 * INV_SQRT_TWO, 0.0,
7171
1.0 - erf(scaled_diff))));
72-
auto lccdf_expr = colwise_sum(log(one_m_erf));
72+
auto lccdf_expr = log(one_m_erf);
7373
auto mu_deriv = select(scaled_diff > 8.25 * INV_SQRT_TWO, INFTY,
7474
SQRT_TWO_OVER_SQRT_PI
7575
* elt_divide(exp(-square(scaled_diff)),
@@ -89,7 +89,7 @@ return_type_t<T_y_cl, T_loc_cl, T_scale_cl> normal_lccdf(
8989
calc_if<!is_constant<T_loc_cl>::value>(mu_deriv),
9090
calc_if<!is_constant<T_scale_cl>::value>(sigma_deriv));
9191

92-
T_partials_return lccdf = LOG_HALF + sum(from_matrix_cl(lccdf_cl));
92+
T_partials_return lccdf = LOG_HALF * lccdf_cl.size() + sum(from_matrix_cl(lccdf_cl));
9393

9494
auto ops_partials = make_partials_propagator(y_col, mu_col, sigma_col);
9595

‎test/unit/math/opencl/rev/normal_lccdf_test.cpp‎

Lines changed: 172 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -5,35 +5,186 @@
55
#include <test/unit/math/opencl/util.hpp>
66
#include <vector>
77

8+
namespace exp_mod_normal_lccdf_test {
89

9-
auto normal_lccdf_functor
10-
= [](const auto& y, const auto& mu, const auto& sigma) {
11-
return stan::math::normal_lccdf(y, mu, sigma);
10+
TEST(ProbDistributionsDoubleExpModNormalLccdf, error_checking) {
11+
int N = 3;
12+
13+
Eigen::VectorXd y(N);
14+
y << 0.3, 0.8, 1.0;
15+
Eigen::VectorXd y_size(N - 1);
16+
y_size << 0.3, 0.8;
17+
Eigen::VectorXd y_value(N);
18+
y_value << 0.3, NAN, 0.5;
19+
20+
Eigen::VectorXd mu(N);
21+
mu << 0.3, 0.8, 1.0;
22+
Eigen::VectorXd mu_size(N - 1);
23+
mu_size << 0.3, 0.8;
24+
Eigen::VectorXd mu_value(N);
25+
mu_value << 0.3, -INFINITY, 0.5;
26+
27+
Eigen::VectorXd sigma(N);
28+
sigma << 0.3, 0.8, 1.0;
29+
Eigen::VectorXd sigma_size(N - 1);
30+
sigma_size << 0.3, 0.8;
31+
Eigen::VectorXd sigma_value(N);
32+
sigma_value << 0.3, 0, 0.5;
33+
34+
Eigen::VectorXd lambda(N);
35+
lambda << 0.4, 0.4, 1.4;
36+
Eigen::VectorXd lambda_size(N - 1);
37+
lambda_size << 0.3, 0.8;
38+
Eigen::VectorXd lambda_value(N);
39+
lambda_value << 0.3, 0, 0.5;
40+
41+
stan::math::matrix_cl<double> y_cl(y);
42+
stan::math::matrix_cl<double> y_size_cl(y_size);
43+
stan::math::matrix_cl<double> y_value_cl(y_value);
44+
stan::math::matrix_cl<double> mu_cl(mu);
45+
stan::math::matrix_cl<double> mu_size_cl(mu_size);
46+
stan::math::matrix_cl<double> mu_value_cl(mu_value);
47+
stan::math::matrix_cl<double> sigma_cl(sigma);
48+
stan::math::matrix_cl<double> sigma_size_cl(sigma_size);
49+
stan::math::matrix_cl<double> sigma_value_cl(sigma_value);
50+
stan::math::matrix_cl<double> lambda_cl(lambda);
51+
stan::math::matrix_cl<double> lambda_size_cl(lambda_size);
52+
stan::math::matrix_cl<double> lambda_value_cl(lambda_value);
53+
54+
EXPECT_NO_THROW(
55+
stan::math::exp_mod_normal_lccdf(y_cl, mu_cl, sigma_cl, lambda_cl));
56+
57+
EXPECT_THROW(
58+
stan::math::exp_mod_normal_lccdf(y_size_cl, mu_cl, sigma_cl, lambda_cl),
59+
std::invalid_argument);
60+
EXPECT_THROW(
61+
stan::math::exp_mod_normal_lccdf(y_cl, mu_size_cl, sigma_cl, lambda_cl),
62+
std::invalid_argument);
63+
EXPECT_THROW(
64+
stan::math::exp_mod_normal_lccdf(y_cl, mu_cl, sigma_size_cl, lambda_cl),
65+
std::invalid_argument);
66+
EXPECT_THROW(
67+
stan::math::exp_mod_normal_lccdf(y_cl, mu_cl, sigma_cl, lambda_size_cl),
68+
std::invalid_argument);
69+
70+
EXPECT_THROW(
71+
stan::math::exp_mod_normal_lccdf(y_value_cl, mu_cl, sigma_cl, lambda_cl),
72+
std::domain_error);
73+
EXPECT_THROW(
74+
stan::math::exp_mod_normal_lccdf(y_cl, mu_value_cl, sigma_cl, lambda_cl),
75+
std::domain_error);
76+
EXPECT_THROW(
77+
stan::math::exp_mod_normal_lccdf(y_cl, mu_cl, sigma_value_cl, lambda_cl),
78+
std::domain_error);
79+
EXPECT_THROW(
80+
stan::math::exp_mod_normal_lccdf(y_cl, mu_cl, sigma_cl, lambda_value_cl),
81+
std::domain_error);
82+
}
83+
84+
auto exp_mod_normal_lccdf_functor
85+
= [](const auto& y, const auto& mu, const auto& sigma, const auto& lambda) {
86+
return stan::math::exp_mod_normal_lccdf(y, mu, sigma, lambda);
1287
};
1388

89+
TEST(ProbDistributionsDoubleExpModNormalLccdf, opencl_matches_cpu_small) {
90+
int N = 3;
91+
int M = 2;
92+
93+
Eigen::VectorXd y(N);
94+
y << -0.3, 1.8, 1.4;
95+
Eigen::VectorXd mu(N);
96+
mu << 0.3, 0.8, 1.0;
97+
Eigen::VectorXd sigma(N);
98+
sigma << 0.3, 0.8, 1.0;
99+
Eigen::VectorXd lambda(N);
100+
lambda << 0.3, 0.4, 1.1;
101+
102+
stan::math::test::compare_cpu_opencl_prim_rev(exp_mod_normal_lccdf_functor, y,
103+
mu, sigma, lambda);
104+
stan::math::test::compare_cpu_opencl_prim_rev(
105+
exp_mod_normal_lccdf_functor, y.transpose().eval(), mu.transpose().eval(),
106+
sigma.transpose().eval(), lambda.transpose().eval());
107+
}
108+
109+
TEST(ProbDistributionsDoubleExpModNormalLccdf,
110+
opencl_matches_cpu_small_y_pos_inf) {
111+
int N = 3;
112+
int M = 2;
14113

114+
Eigen::VectorXd y(N);
115+
y << -0.3, 1.8, INFINITY;
116+
Eigen::VectorXd mu(N);
117+
mu << 0.3, 0.8, 1.0;
118+
Eigen::VectorXd sigma(N);
119+
sigma << 0.3, 0.8, 1.0;
120+
Eigen::VectorXd lambda(N);
121+
lambda << 0.3, 0.4, 1.1;
15122

16-
TEST(ProbDistributionsNormalLccdf, opencl_matches_cpu_big) {
123+
stan::math::test::compare_cpu_opencl_prim_rev(exp_mod_normal_lccdf_functor, y,
124+
mu, sigma, lambda);
125+
stan::math::test::compare_cpu_opencl_prim_rev(
126+
exp_mod_normal_lccdf_functor, y.transpose().eval(), mu.transpose().eval(),
127+
sigma.transpose().eval(), lambda.transpose().eval());
128+
}
129+
130+
TEST(ProbDistributionsDoubleExpModNormalLccdf,
131+
opencl_matches_cpu_small_y_neg_inf) {
132+
int N = 3;
133+
int M = 2;
134+
135+
Eigen::VectorXd y(N);
136+
y << -0.3, 1.8, -INFINITY;
137+
Eigen::VectorXd mu(N);
138+
mu << 0.3, 0.8, 1.0;
139+
Eigen::VectorXd sigma(N);
140+
sigma << 0.3, 0.8, 1.0;
141+
Eigen::VectorXd lambda(N);
142+
lambda << 0.3, 0.4, 1.1;
143+
144+
stan::math::test::compare_cpu_opencl_prim_rev(exp_mod_normal_lccdf_functor, y,
145+
mu, sigma, lambda);
146+
stan::math::test::compare_cpu_opencl_prim_rev(
147+
exp_mod_normal_lccdf_functor, y.transpose().eval(), mu.transpose().eval(),
148+
sigma.transpose().eval(), lambda.transpose().eval());
149+
}
150+
151+
TEST(ProbDistributionsDoubleExpModNormalLccdf, opencl_broadcast_y) {
152+
int N = 3;
153+
154+
double y_scal = 12.3;
155+
Eigen::VectorXd mu(N);
156+
mu << 0.5, 1.2, 1.0;
157+
Eigen::VectorXd sigma(N);
158+
sigma << 0.3, 0.8, 1.0;
159+
Eigen::VectorXd lambda(N);
160+
lambda << 0.3, 0.4, 1.1;
161+
162+
stan::math::test::test_opencl_broadcasting_prim_rev<0>(
163+
exp_mod_normal_lccdf_functor, y_scal, mu, sigma, lambda);
164+
stan::math::test::test_opencl_broadcasting_prim_rev<0>(
165+
exp_mod_normal_lccdf_functor, y_scal, mu.transpose().eval(), sigma,
166+
lambda.transpose().eval());
167+
}
168+
169+
TEST(ProbDistributionsDoubleExpModNormalLccdf, opencl_matches_cpu_big) {
17170
int N = 153;
18171

19-
std::srand(123);
20-
for (int i = 0; i < 10; ++i) {
172+
Eigen::Matrix<double, Eigen::Dynamic, 1> y
173+
= Eigen::Array<double, Eigen::Dynamic, 1>::Random(N, 1).abs();
21174
Eigen::Matrix<double, Eigen::Dynamic, 1> mu
22-
= Eigen::Array<double, Eigen::Dynamic, 1>::Random(N, 1) + 1.0;
175+
= Eigen::Array<double, Eigen::Dynamic, 1>::Random(N, 1).abs();
23176
Eigen::Matrix<double, Eigen::Dynamic, 1> sigma
24-
= Eigen::Array<double, Eigen::Dynamic, 1>::Random(N, 1).abs() + 0.01;
25-
Eigen::Matrix<double, Eigen::Dynamic, 1> y = (mu.array() * sigma.array()).matrix();
26-
std::cout << "Iter: " << i << " mu, sigma, y" << std::endl;
27-
for (int j = 0; j < N; j++) {
28-
std::cout << mu(j) << ", " << sigma(j) << ", " << y(j) << std::endl;
29-
}
30-
std::cout << "-----------compare_cpu_opencl_prim_rev" << std::endl;
31-
stan::math::test::compare_cpu_opencl_prim_rev(normal_lccdf_functor, y, mu,
32-
sigma);
33-
std::cout << "-----------compare_cpu_opencl_prim_rev transpose" << std::endl;
177+
= Eigen::Array<double, Eigen::Dynamic, 1>::Random(N, 1).abs().array()
178+
+ 0.1;
179+
Eigen::Matrix<double, Eigen::Dynamic, 1> lambda
180+
= Eigen::Array<double, Eigen::Dynamic, 1>::Random(N, 1).abs();
181+
182+
stan::math::test::compare_cpu_opencl_prim_rev(exp_mod_normal_lccdf_functor, y,
183+
mu, sigma, lambda);
34184
stan::math::test::compare_cpu_opencl_prim_rev(
35-
normal_lccdf_functor, y.transpose().eval(), mu.transpose().eval(),
36-
sigma.transpose().eval());
37-
}
185+
exp_mod_normal_lccdf_functor, y.transpose().eval(), mu.transpose().eval(),
186+
sigma.transpose().eval(), lambda.transpose().eval());
38187
}
39-
#endif
188+
} // namespace exp_mod_normal_lccdf_test
189+
190+
#endif

0 commit comments

Comments
 (0)