Skip to content

Commit dbfa52c

Browse files
GidonFrischkornGidonFrischkorn
committed
Merge branch 'develop' of https://github.com/stan-dev/math into feature/issue-2788-inv-gaussian
2 parents 5e6f41b + d20b238 commit dbfa52c

39 files changed

Lines changed: 2365 additions & 1270 deletions

‎Jenkinsfile‎

100755100644
Lines changed: 261 additions & 590 deletions
Large diffs are not rendered by default.

‎runTests.py‎

Lines changed: 26 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -116,6 +116,13 @@ def processCLIArgs():
116116
action="store_true",
117117
help="Don't run tests, just try to make them.",
118118
)
119+
parser.add_argument(
120+
"-t",
121+
"--test-only",
122+
dest="test_only",
123+
action="store_true",
124+
help="Don't make tests, just try to run them (must already be made).",
125+
)
119126
parser.add_argument(
120127
"--run-all",
121128
dest="run_all",
@@ -405,20 +412,21 @@ def main():
405412
if inputs.changed:
406413
tests = findChangedTests(inputs.debug)
407414
else:
408-
# pass 0: generate all auto-generated tests
409-
if any("test/prob" in arg for arg in inputs.tests):
410-
generateTests(inputs.j)
411-
412-
if inputs.do_jumbo:
413-
jumboFiles = generateJumboTests(inputs.tests)
414-
if inputs.e == -1:
415-
if inputs.j == 1:
416-
num_expr_test_files = 1
415+
if not inputs.test_only:
416+
# pass 0: generate all auto-generated tests
417+
if any("test/prob" in arg for arg in inputs.tests):
418+
generateTests(inputs.j)
419+
420+
if inputs.do_jumbo:
421+
jumboFiles = generateJumboTests(inputs.tests)
422+
if inputs.e == -1:
423+
if inputs.j == 1:
424+
num_expr_test_files = 1
425+
else:
426+
num_expr_test_files = inputs.j * 4
417427
else:
418-
num_expr_test_files = inputs.j * 4
419-
else:
420-
num_expr_test_files = inputs.e
421-
handleExpressionTests(inputs.tests, inputs.only_functions, num_expr_test_files)
428+
num_expr_test_files = inputs.e
429+
handleExpressionTests(inputs.tests, inputs.only_functions, num_expr_test_files)
422430

423431
tests = findTests(inputs.tests, inputs.f, inputs.do_jumbo)
424432

@@ -431,10 +439,11 @@ def main():
431439

432440
try:
433441
# pass 1: make test executables
434-
for batch in batched(tests):
435-
if inputs.debug:
436-
print("Test batch: ", batch)
437-
makeTest(" ".join(batch), inputs.j)
442+
if not inputs.test_only:
443+
for batch in batched(tests):
444+
if inputs.debug:
445+
print("Test batch: ", batch)
446+
makeTest(" ".join(batch), inputs.j)
438447
if not inputs.make_only:
439448
# pass 2: run test targets
440449
for t in tests:

‎stan/math/opencl/kernels/device_functions/std_normal_lcdf.hpp‎

Lines changed: 19 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ static constexpr const char* std_normal_lcdf_device_function
2929
if (isnan(lcdf_n)) {
3030
lcdf_n = 0;
3131
}
32-
} else if (scaled_y > -20.0) {
32+
} else if (scaled_y > -4.0) {
3333
// CDF(x) = 1/2 - 1/2 erf(-x) = 1/2 erfc(-x)
3434
lcdf_n = log(erfc(-scaled_y)) - M_LN2;
3535
} else if (10.0 * log(fabs(scaled_y)) < log(DBL_MAX)) {
@@ -76,10 +76,15 @@ static constexpr const char* std_normal_lcdf_device_function
7676
t = 1.0 / (1.0 + 0.3275911 * scaled_y);
7777
t2 = t * t;
7878
t4 = pow(t, 4);
79-
dnlcdf = 0.5 * M_2_SQRTPI
80-
/ (exp(x2) - 0.254829592 + 0.284496736 * t
81-
- 1.421413741 * t2 + 1.453152027 * t2 * t
82-
- 1.061405429 * t4);
79+
// A&S 7.1.26 keeps exp(-x2) in the numerator, as R's pnorm
80+
// does; refs in stan/math/prim/prob/std_normal_lcdf.hpp
81+
const double exp_m_x2 = exp(-x2);
82+
dnlcdf
83+
= 0.5 * M_2_SQRTPI * exp_m_x2
84+
/ (1.0
85+
- exp_m_x2
86+
* (0.254829592 - 0.284496736 * t + 1.421413741 * t2
87+
- 1.453152027 * t2 * t + 1.061405429 * t4));
8388
} else if (scaled_y > 2.5) {
8489
t = scaled_y - 2.7;
8590
t2 = t * t;
@@ -116,6 +121,14 @@ static constexpr const char* std_normal_lcdf_device_function
116121
dnlcdf = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2
117122
+ 0.04700850676 * t2 * t - 0.03478651979 * t4
118123
- 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6);
124+
} else if (scaled_y < -29.0) {
125+
// asymptotic Mills ratio, DLMF 7.12.1; grows linearly as
126+
// -2*scaled_y, same 1/x^2 series shape as R's pnorm uses
127+
const double inv_x2 = 1.0 / x2;
128+
dnlcdf
129+
= -2.0 * scaled_y
130+
/ (1.0
131+
+ inv_x2 * (-0.5 + inv_x2 * (0.75 + inv_x2 * -1.875)));
119132
} else if (10.0 * log(fabs(scaled_y)) < log(DBL_MAX)) {
120133
t = 1.0 / (1.0 - 0.3275911 * scaled_y);
121134
t2 = t * t;
@@ -124,10 +137,7 @@ static constexpr const char* std_normal_lcdf_device_function
124137
= M_2_SQRTPI
125138
/ (0.254829592 * t - 0.284496736 * t2 + 1.421413741 * t2 * t
126139
- 1.453152027 * t4 + 1.061405429 * t4 * t);
127-
if (scaled_y < -29.0) {
128-
dnlcdf += 0.0015065154280332 * x2
129-
- 0.3993154819705530 * scaled_y - 4.2919418242931700;
130-
} else if (scaled_y < -17.0) {
140+
if (scaled_y < -17.0) {
131141
dnlcdf += 0.0001263257217272 * x2 * scaled_y
132142
+ 0.0123586859488623 * x2
133143
- 0.0860505264736028 * scaled_y - 1.252783383752970;

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

Lines changed: 5 additions & 78 deletions
Original file line numberDiff line numberDiff line change
@@ -2,16 +2,13 @@
22
#define STAN_MATH_OPENCL_PRIM_NORMAL_LCCDF_HPP
33
#ifdef STAN_OPENCL
44

5-
#include <stan/math/prim/meta.hpp>
6-
#include <stan/math/prim/err.hpp>
7-
#include <stan/math/prim/fun/constants.hpp>
8-
#include <stan/math/prim/fun/elt_divide.hpp>
9-
#include <stan/math/prim/fun/elt_multiply.hpp>
10-
#include <stan/math/opencl/kernel_generator.hpp>
11-
#include <stan/math/prim/functor/partials_propagator.hpp>
5+
#include <stan/math/opencl/prim/normal_lcdf.hpp>
126

137
namespace stan {
148
namespace math {
9+
namespace internal {
10+
constexpr char normal_lccdf_opencl_func[] = "normal_lccdf(OpenCL)";
11+
} // namespace internal
1512

1613
/** \ingroup opencl
1714
* Returns the normal log complementary cumulative distribution function
@@ -33,77 +30,7 @@ template <
3330
require_any_not_stan_scalar_t<T_y_cl, T_loc_cl, T_scale_cl>* = nullptr>
3431
inline return_type_t<T_y_cl, T_loc_cl, T_scale_cl> normal_lccdf(
3532
const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) {
36-
static constexpr const char* function = "normal_lccdf(OpenCL)";
37-
using T_partials_return = partials_return_t<T_y_cl, T_loc_cl, T_scale_cl>;
38-
using std::isfinite;
39-
using std::isnan;
40-
41-
check_consistent_sizes(function, "Random variable", y, "Location parameter",
42-
mu, "Scale parameter", sigma);
43-
const size_t N = max_size(y, mu, sigma);
44-
if (N == 0) {
45-
return 0.0;
46-
}
47-
48-
const auto& y_col = as_column_vector_or_scalar(y);
49-
const auto& mu_col = as_column_vector_or_scalar(mu);
50-
const auto& sigma_col = as_column_vector_or_scalar(sigma);
51-
52-
const auto& y_val = value_of(y_col);
53-
const auto& mu_val = value_of(mu_col);
54-
const auto& sigma_val = value_of(sigma_col);
55-
56-
auto check_y_not_nan
57-
= check_cl(function, "Random variable", y_val, "not NaN");
58-
auto y_not_nan_expr = !isnan(y_val);
59-
auto check_mu_finite
60-
= check_cl(function, "Location parameter", mu_val, "finite");
61-
auto mu_finite_expr = isfinite(mu_val);
62-
auto check_sigma_positive
63-
= check_cl(function, "Scale parameter", sigma_val, "positive");
64-
auto sigma_positive_expr = 0 < sigma_val;
65-
66-
auto scaled_diff = elt_divide(y_val - mu_val, sigma_val * SQRT_TWO);
67-
matrix_cl<double> one_m_erf = select(
68-
scaled_diff < -37.5 * INV_SQRT_TWO, 2.0,
69-
select(scaled_diff < -5.0 * INV_SQRT_TWO, 2.0 - erfc(-scaled_diff),
70-
select(scaled_diff > 8.25 * INV_SQRT_TWO, 0.0,
71-
1.0 - erf(scaled_diff))));
72-
auto lccdf_expr = log(one_m_erf);
73-
auto mu_deriv = select(scaled_diff > 8.25 * INV_SQRT_TWO, INFTY,
74-
SQRT_TWO_OVER_SQRT_PI
75-
* elt_divide(exp(-square(scaled_diff)),
76-
elt_multiply(one_m_erf, sigma_val)));
77-
auto y_deriv = -mu_deriv;
78-
auto sigma_deriv = elt_multiply(mu_deriv, scaled_diff) * SQRT_TWO;
79-
80-
matrix_cl<double> lccdf_cl;
81-
matrix_cl<double> y_deriv_cl;
82-
matrix_cl<double> mu_deriv_cl;
83-
matrix_cl<double> sigma_deriv_cl;
84-
85-
results(check_y_not_nan, check_mu_finite, check_sigma_positive)
86-
= expressions(y_not_nan_expr, mu_finite_expr, sigma_positive_expr);
87-
results(lccdf_cl, y_deriv_cl, mu_deriv_cl, sigma_deriv_cl)
88-
= expressions(lccdf_expr, calc_if<is_autodiff_v<T_y_cl>>(y_deriv),
89-
calc_if<is_autodiff_v<T_loc_cl>>(mu_deriv),
90-
calc_if<is_autodiff_v<T_scale_cl>>(sigma_deriv));
91-
92-
T_partials_return lccdf
93-
= LOG_HALF * lccdf_cl.size() + sum(from_matrix_cl(lccdf_cl));
94-
95-
auto ops_partials = make_partials_propagator(y_col, mu_col, sigma_col);
96-
97-
if constexpr (is_autodiff_v<T_y_cl>) {
98-
partials<0>(ops_partials) = std::move(y_deriv_cl);
99-
}
100-
if constexpr (is_autodiff_v<T_loc_cl>) {
101-
partials<1>(ops_partials) = std::move(mu_deriv_cl);
102-
}
103-
if constexpr (is_autodiff_v<T_scale_cl>) {
104-
partials<2>(ops_partials) = std::move(sigma_deriv_cl);
105-
}
106-
return ops_partials.build(lccdf);
33+
return normal_lcdf<internal::normal_lccdf_opencl_func>(-y, -mu, sigma);
10734
}
10835

10936
} // namespace math

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

Lines changed: 27 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -13,6 +13,7 @@
1313
namespace stan {
1414
namespace math {
1515
namespace internal {
16+
constexpr char normal_lcdf_opencl_func[] = "normal_lcdf(OpenCL)";
1617
const char opencl_normal_lcdf_impl[] = STRINGIFY(
1718
double x2 = normal_lcdf_scaled_diff * normal_lcdf_scaled_diff;
1819
double normal_lcdf_n = 0;
@@ -29,7 +30,7 @@ const char opencl_normal_lcdf_impl[] = STRINGIFY(
2930
if (isnan(normal_lcdf_n)) {
3031
normal_lcdf_n = 0;
3132
}
32-
} else if (normal_lcdf_scaled_diff > -20.0) {
33+
} else if (normal_lcdf_scaled_diff > -4.0) {
3334
// CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x)
3435
normal_lcdf_n = log(erfc(-normal_lcdf_scaled_diff)) - M_LN2;
3536
} else if (10.0 * log(fabs(normal_lcdf_scaled_diff)) < log(DBL_MAX)) {
@@ -58,7 +59,8 @@ const char opencl_normal_lcdf_impl[] = STRINGIFY(
5859
// NOLINTBEGIN
5960
const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY(
6061
double normal_ldncdf = 0.0; double t = 0.0; double t2 = 0.0;
61-
double t4 = 0.0;
62+
double t4 = 0.0; double normal_lcdf_exp_m_x2 = 0.0;
63+
double normal_lcdf_inv_x2 = 0.0;
6264

6365
// calculate using piecewise function
6466
// (due to instability / inaccuracy in the various approximations)
@@ -67,10 +69,15 @@ const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY(
6769
t = 1.0 / (1.0 + 0.3275911 * normal_lcdf_deriv_scaled_diff);
6870
t2 = t * t;
6971
t4 = pow(t, 4);
72+
// A&S 7.1.26 keeps exp(-x2) in the numerator, as R's pnorm do_del
73+
// does; refs in stan/math/prim/prob/normal_lcdf.hpp
74+
normal_lcdf_exp_m_x2 = exp(-x2);
7075
normal_ldncdf
71-
= 0.5 * M_2_SQRTPI
72-
/ (exp(x2) - 0.254829592 + 0.284496736 * t - 1.421413741 * t2
73-
+ 1.453152027 * t2 * t - 1.061405429 * t4);
76+
= 0.5 * M_2_SQRTPI * normal_lcdf_exp_m_x2
77+
/ (1.0
78+
- normal_lcdf_exp_m_x2
79+
* (0.254829592 - 0.284496736 * t + 1.421413741 * t2
80+
- 1.453152027 * t2 * t + 1.061405429 * t4));
7481
} else if (normal_lcdf_deriv_scaled_diff > 2.5) {
7582
// in the trouble area where all of the standard numerical
7683
// approximations are unstable - bridge the gap using Taylor
@@ -114,6 +121,17 @@ const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY(
114121
normal_ldncdf = 0.6245634904 - 0.9521866949 * t + 0.3986215682 * t2
115122
+ 0.04700850676 * t2 * t - 0.03478651979 * t4
116123
- 0.01772675404 * t4 * t + 0.0006577254811 * pow(t, 6);
124+
} else if (normal_lcdf_deriv_scaled_diff < -29.0) {
125+
// asymptotic Mills ratio, DLMF 7.12.1; grows linearly as -2*scaled_diff,
126+
// so no quadratic fit can track it. Same 1/x^2 series shape as R's pnorm
127+
normal_lcdf_inv_x2 = 1.0 / x2;
128+
normal_ldncdf
129+
= -2.0 * normal_lcdf_deriv_scaled_diff
130+
/ (1.0
131+
+ normal_lcdf_inv_x2
132+
* (-0.5
133+
+ normal_lcdf_inv_x2
134+
* (0.75 + normal_lcdf_inv_x2 * -1.875)));
117135
} else if (10.0 * log(fabs(normal_lcdf_deriv_scaled_diff)) < log(DBL_MAX)) {
118136
// approximation derived from Abramowitz and Stegun (1964) 7.1.26
119137
// use fact that erf(x)=-erf(-x)
@@ -128,11 +146,7 @@ const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY(
128146
- 1.453152027 * t4 + 1.061405429 * t4 * t);
129147
// check if we need to add a correction term
130148
// (from cubic fit of residuals)
131-
if (normal_lcdf_deriv_scaled_diff < -29.0) {
132-
normal_ldncdf += 0.0015065154280332 * x2
133-
- 0.3993154819705530 * normal_lcdf_deriv_scaled_diff
134-
- 4.2919418242931700;
135-
} else if (normal_lcdf_deriv_scaled_diff < -17.0) {
149+
if (normal_lcdf_deriv_scaled_diff < -17.0) {
136150
normal_ldncdf += 0.0001263257217272 * x2 * normal_lcdf_deriv_scaled_diff
137151
+ 0.0123586859488623 * x2
138152
- 0.0860505264736028 * normal_lcdf_deriv_scaled_diff
@@ -174,13 +188,14 @@ const char opencl_normal_lcdf_ldncdf_impl[] = STRINGIFY(
174188
* @return The log of the product of densities.
175189
*/
176190
template <
177-
typename T_y_cl, typename T_loc_cl, typename T_scale_cl,
191+
const char* func = internal::normal_lcdf_opencl_func, typename T_y_cl,
192+
typename T_loc_cl, typename T_scale_cl,
178193
require_all_prim_or_rev_kernel_expression_t<T_y_cl, T_loc_cl,
179194
T_scale_cl>* = nullptr,
180195
require_any_not_stan_scalar_t<T_y_cl, T_loc_cl, T_scale_cl>* = nullptr>
181196
inline return_type_t<T_y_cl, T_loc_cl, T_scale_cl> normal_lcdf(
182197
const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) {
183-
static constexpr const char* function = "normal_lcdf(OpenCL)";
198+
static constexpr const char* function = func;
184199
using std::isfinite;
185200
using std::isnan;
186201

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

Lines changed: 5 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -2,16 +2,13 @@
22
#define STAN_MATH_OPENCL_PRIM_STD_NORMAL_LCCDF_HPP
33
#ifdef STAN_OPENCL
44

5-
#include <stan/math/prim/meta.hpp>
6-
#include <stan/math/prim/err.hpp>
7-
#include <stan/math/prim/fun/constants.hpp>
8-
#include <stan/math/prim/fun/elt_divide.hpp>
9-
#include <stan/math/prim/fun/elt_multiply.hpp>
10-
#include <stan/math/opencl/kernel_generator.hpp>
11-
#include <stan/math/prim/functor/partials_propagator.hpp>
5+
#include <stan/math/opencl/prim/std_normal_lcdf.hpp>
126

137
namespace stan {
148
namespace math {
9+
namespace internal {
10+
constexpr char std_normal_lccdf_opencl_func[] = "std_normal_lccdf(OpenCL)";
11+
} // namespace internal
1512

1613
/** \ingroup opencl
1714
* Returns the log standard normal complementary cumulative distribution
@@ -25,47 +22,7 @@ template <typename T_y_cl,
2522
require_all_prim_or_rev_kernel_expression_t<T_y_cl>* = nullptr,
2623
require_any_not_stan_scalar_t<T_y_cl>* = nullptr>
2724
inline return_type_t<T_y_cl> std_normal_lccdf(const T_y_cl& y) {
28-
static constexpr const char* function = "std_normal_lccdf(OpenCL)";
29-
using T_partials_return = partials_return_t<T_y_cl>;
30-
using std::isfinite;
31-
using std::isnan;
32-
33-
const size_t N = math::size(y);
34-
if (N == 0) {
35-
return 1.0;
36-
}
37-
38-
const auto& y_col = as_column_vector_or_scalar(y);
39-
const auto& y_val = value_of(y_col);
40-
41-
auto check_y_not_nan
42-
= check_cl(function, "Random variable", y_val, "not NaN");
43-
auto y_not_nan_expr = !isnan(y_val);
44-
45-
auto scaled_y = y_val * INV_SQRT_TWO;
46-
auto one_m_erf
47-
= select(y_val < -37.5, 2.0,
48-
select(y_val < -5.0, 2.0 - erfc(-scaled_y),
49-
select(y_val > 8.25, 0.0, 1.0 - erf(scaled_y))));
50-
auto lccdf_expr = colwise_sum(log(one_m_erf));
51-
auto y_deriv = -select(
52-
y_val > 8.25, INFTY,
53-
SQRT_TWO_OVER_SQRT_PI * elt_divide(exp(-square(scaled_y)), one_m_erf));
54-
55-
matrix_cl<double> lccdf_cl;
56-
matrix_cl<double> y_deriv_cl;
57-
58-
results(check_y_not_nan, lccdf_cl, y_deriv_cl) = expressions(
59-
y_not_nan_expr, lccdf_expr, calc_if<is_autodiff_v<T_y_cl>>(y_deriv));
60-
61-
T_partials_return lccdf = from_matrix_cl(lccdf_cl).sum() + LOG_HALF * N;
62-
63-
auto ops_partials = make_partials_propagator(y_col);
64-
65-
if constexpr (is_autodiff_v<T_y_cl>) {
66-
partials<0>(ops_partials) = std::move(y_deriv_cl);
67-
}
68-
return ops_partials.build(lccdf);
25+
return std_normal_lcdf<internal::std_normal_lccdf_opencl_func>(-y);
6926
}
7027

7128
} // namespace math

0 commit comments

Comments
 (0)