Skip to content

Commit f390d82

Browse files
rok-cesnovarrok-cesnovar
authored
Merge pull request #2449 from bstatcomp/opencl_double_exponential_cdf
Opencl double_exponential and exp_mod_normal (lc)cdf
2 parents 1bf9657 + 1e4a550 commit f390d82

26 files changed

Lines changed: 2070 additions & 6 deletions

‎stan/math/opencl/prim.hpp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -134,7 +134,13 @@
134134
#include <stan/math/opencl/prim/divide_columns.hpp>
135135
#include <stan/math/opencl/prim/dot_product.hpp>
136136
#include <stan/math/opencl/prim/dot_self.hpp>
137+
#include <stan/math/opencl/prim/double_exponential_cdf.hpp>
138+
#include <stan/math/opencl/prim/double_exponential_lccdf.hpp>
139+
#include <stan/math/opencl/prim/double_exponential_lcdf.hpp>
137140
#include <stan/math/opencl/prim/double_exponential_lpdf.hpp>
141+
#include <stan/math/opencl/prim/exp_mod_normal_cdf.hpp>
142+
#include <stan/math/opencl/prim/exp_mod_normal_lccdf.hpp>
143+
#include <stan/math/opencl/prim/exp_mod_normal_lcdf.hpp>
138144
#include <stan/math/opencl/prim/exp_mod_normal_lpdf.hpp>
139145
#include <stan/math/opencl/prim/exponential_cdf.hpp>
140146
#include <stan/math/opencl/prim/exponential_lccdf.hpp>

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

Lines changed: 118 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,118 @@
1+
#ifndef STAN_MATH_OPENCL_PRIM_DOUBLE_EXPONENTIAL_CDF_HPP
2+
#define STAN_MATH_OPENCL_PRIM_DOUBLE_EXPONENTIAL_CDF_HPP
3+
#ifdef STAN_OPENCL
4+
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/operands_and_partials.hpp>
12+
13+
namespace stan {
14+
namespace math {
15+
16+
/** \ingroup opencl
17+
* Returns the double exponential cumulative density function. Given
18+
* containers of matching sizes, returns the product of probabilities.
19+
*
20+
* @tparam T_y_cl type of scalar outcome
21+
* @tparam T_loc_cl type of location
22+
* @tparam T_scale_cl type of scale
23+
* @param y (Sequence of) scalar(s).
24+
* @param mu (Sequence of) location(s).
25+
* @param sigma (Sequence of) scale(s).
26+
* @return The log of the product of densities.
27+
*/
28+
template <
29+
typename T_y_cl, typename T_loc_cl, typename T_scale_cl,
30+
require_all_prim_or_rev_kernel_expression_t<T_y_cl, T_loc_cl,
31+
T_scale_cl>* = nullptr,
32+
require_any_not_stan_scalar_t<T_y_cl, T_loc_cl, T_scale_cl>* = nullptr>
33+
return_type_t<T_y_cl, T_loc_cl, T_scale_cl> double_exponential_cdf(
34+
const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) {
35+
static const char* function = "double_exponential_cdf(OpenCL)";
36+
using T_partials_return = partials_return_t<T_y_cl, T_loc_cl, T_scale_cl>;
37+
using std::isfinite;
38+
using std::isnan;
39+
40+
check_consistent_sizes(function, "Random variable", y, "Location parameter",
41+
mu, "Scale parameter", sigma);
42+
const size_t N = max_size(y, mu, sigma);
43+
if (N == 0) {
44+
return 1.0;
45+
}
46+
47+
const auto& y_col = as_column_vector_or_scalar(y);
48+
const auto& mu_col = as_column_vector_or_scalar(mu);
49+
const auto& sigma_col = as_column_vector_or_scalar(sigma);
50+
51+
const auto& y_val = value_of(y_col);
52+
const auto& mu_val = value_of(mu_col);
53+
const auto& sigma_val = value_of(sigma_col);
54+
55+
auto check_y_not_nan
56+
= check_cl(function, "Random variable", y_val, "not NaN");
57+
auto y_not_nan_expr = !isnan(y_val);
58+
auto check_mu_finite
59+
= check_cl(function, "Location parameter", mu_val, "finite");
60+
auto mu_finite_expr = isfinite(mu_val);
61+
auto check_sigma_positive_finite
62+
= check_cl(function, "Scale parameter", sigma_val, "positive finite");
63+
auto sigma_positive_finite_expr = 0 < sigma_val && isfinite(sigma_val);
64+
65+
auto scaled_diff = elt_divide(y_val - mu_val, sigma_val);
66+
auto exp_scaled_diff = exp(scaled_diff);
67+
auto cond = y_val < mu_val;
68+
auto cdf_expr = colwise_prod(select(cond, 0.5 * exp_scaled_diff,
69+
1.0 - elt_divide(0.5, exp_scaled_diff)));
70+
71+
matrix_cl<double> cdf_cl;
72+
matrix_cl<double> mu_deriv_cl; // also temporarily stores exp_scaled_diff
73+
matrix_cl<double> y_deriv_cl;
74+
matrix_cl<double> sigma_deriv_cl; // also temporarily stores scaled_diff
75+
76+
results(check_y_not_nan, check_mu_finite, check_sigma_positive_finite, cdf_cl,
77+
mu_deriv_cl, sigma_deriv_cl)
78+
= expressions(
79+
y_not_nan_expr, mu_finite_expr, sigma_positive_finite_expr, cdf_expr,
80+
calc_if<!is_constant_all<T_y_cl, T_loc_cl, T_scale_cl>::value>(
81+
exp_scaled_diff),
82+
calc_if<!is_constant<T_scale_cl>::value>(scaled_diff));
83+
84+
T_partials_return cdf = (from_matrix_cl(cdf_cl)).prod();
85+
86+
operands_and_partials<decltype(y_col), decltype(mu_col), decltype(sigma_col)>
87+
ops_partials(y_col, mu_col, sigma_col);
88+
if (!is_constant_all<T_y_cl, T_loc_cl, T_scale_cl>::value) {
89+
auto cdf_div_sigma = elt_divide(cdf, sigma_val);
90+
auto y_deriv = select(cond, cdf_div_sigma,
91+
elt_divide(cdf_div_sigma, (2.0 * mu_deriv_cl - 1.0)));
92+
auto mu_deriv = -y_deriv;
93+
auto sigma_deriv
94+
= elt_multiply(mu_deriv, static_select<is_constant<T_scale_cl>::value>(
95+
mu_deriv_cl, sigma_deriv_cl));
96+
97+
results(mu_deriv_cl, y_deriv_cl, sigma_deriv_cl)
98+
= expressions(calc_if<!is_constant<T_loc_cl>::value>(mu_deriv),
99+
calc_if<!is_constant<T_y_cl>::value>(y_deriv),
100+
calc_if<!is_constant<T_scale_cl>::value>(sigma_deriv));
101+
102+
if (!is_constant<T_y_cl>::value) {
103+
ops_partials.edge1_.partials_ = std::move(y_deriv_cl);
104+
}
105+
if (!is_constant<T_loc_cl>::value) {
106+
ops_partials.edge2_.partials_ = std::move(mu_deriv_cl);
107+
}
108+
if (!is_constant<T_scale_cl>::value) {
109+
ops_partials.edge3_.partials_ = std::move(sigma_deriv_cl);
110+
}
111+
}
112+
return ops_partials.build(cdf);
113+
}
114+
115+
} // namespace math
116+
} // namespace stan
117+
#endif
118+
#endif

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

Lines changed: 109 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,109 @@
1+
#ifndef STAN_MATH_OPENCL_PRIM_DOUBLE_EXPONENTIAL_LCCDF_HPP
2+
#define STAN_MATH_OPENCL_PRIM_DOUBLE_EXPONENTIAL_LCCDF_HPP
3+
#ifdef STAN_OPENCL
4+
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/operands_and_partials.hpp>
12+
13+
namespace stan {
14+
namespace math {
15+
16+
/** \ingroup opencl
17+
* Returns the double exponential log complementary cumulative density
18+
* function. Given containers of matching sizes, returns the log sum of
19+
* probabilities.
20+
*
21+
* @tparam T_y_cl type of scalar outcome
22+
* @tparam T_loc_cl type of location
23+
* @tparam T_scale_cl type of scale
24+
* @param y (Sequence of) scalar(s).
25+
* @param mu (Sequence of) location(s).
26+
* @param sigma (Sequence of) scale(s).
27+
* @return The log of the product of densities.
28+
*/
29+
template <
30+
typename T_y_cl, typename T_loc_cl, typename T_scale_cl,
31+
require_all_prim_or_rev_kernel_expression_t<T_y_cl, T_loc_cl,
32+
T_scale_cl>* = nullptr,
33+
require_any_not_stan_scalar_t<T_y_cl, T_loc_cl, T_scale_cl>* = nullptr>
34+
return_type_t<T_y_cl, T_loc_cl, T_scale_cl> double_exponential_lccdf(
35+
const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) {
36+
static const char* function = "double_exponential_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_finite
63+
= check_cl(function, "Scale parameter", sigma_val, "positive finite");
64+
auto sigma_positive_finite_expr = 0 < sigma_val && isfinite(sigma_val);
65+
66+
auto sigma_inv = elt_divide(1.0, sigma_val);
67+
auto scaled_diff = elt_multiply(y_val - mu_val, sigma_inv);
68+
auto cond = y_val < mu_val;
69+
auto mu_deriv = select(
70+
cond, elt_divide(sigma_inv, 2.0 * exp(-scaled_diff) - 1.0), sigma_inv);
71+
auto ccdf_log_expr = colwise_sum(
72+
select(cond, log1m(0.5 * exp(scaled_diff)), LOG_HALF - scaled_diff));
73+
auto y_deriv = -mu_deriv;
74+
auto sigma_deriv = elt_multiply(mu_deriv, scaled_diff);
75+
76+
matrix_cl<double> ccdf_log_cl;
77+
matrix_cl<double> mu_deriv_cl;
78+
matrix_cl<double> y_deriv_cl;
79+
matrix_cl<double> sigma_deriv_cl;
80+
81+
results(check_y_not_nan, check_mu_finite, check_sigma_positive_finite,
82+
ccdf_log_cl, y_deriv_cl, mu_deriv_cl, sigma_deriv_cl)
83+
= expressions(y_not_nan_expr, mu_finite_expr, sigma_positive_finite_expr,
84+
ccdf_log_expr,
85+
calc_if<!is_constant<T_y_cl>::value>(y_deriv),
86+
calc_if<!is_constant<T_loc_cl>::value>(mu_deriv),
87+
calc_if<!is_constant<T_scale_cl>::value>(sigma_deriv));
88+
89+
T_partials_return ccdf_log = (from_matrix_cl(ccdf_log_cl)).sum();
90+
91+
operands_and_partials<decltype(y_col), decltype(mu_col), decltype(sigma_col)>
92+
ops_partials(y_col, mu_col, sigma_col);
93+
94+
if (!is_constant<T_y_cl>::value) {
95+
ops_partials.edge1_.partials_ = std::move(y_deriv_cl);
96+
}
97+
if (!is_constant<T_loc_cl>::value) {
98+
ops_partials.edge2_.partials_ = std::move(mu_deriv_cl);
99+
}
100+
if (!is_constant<T_scale_cl>::value) {
101+
ops_partials.edge3_.partials_ = std::move(sigma_deriv_cl);
102+
}
103+
return ops_partials.build(ccdf_log);
104+
}
105+
106+
} // namespace math
107+
} // namespace stan
108+
#endif
109+
#endif

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

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
#ifndef STAN_MATH_OPENCL_PRIM_DOUBLE_EXPONENTIAL_LCDF_HPP
2+
#define STAN_MATH_OPENCL_PRIM_DOUBLE_EXPONENTIAL_LCDF_HPP
3+
#ifdef STAN_OPENCL
4+
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/operands_and_partials.hpp>
12+
13+
namespace stan {
14+
namespace math {
15+
16+
/** \ingroup opencl
17+
* Returns the double exponential log cumulative density function. Given
18+
* containers of matching sizes, returns the log sum of probabilities.
19+
*
20+
* @tparam T_y_cl type of scalar outcome
21+
* @tparam T_loc_cl type of location
22+
* @tparam T_scale_cl type of scale
23+
* @param y (Sequence of) scalar(s).
24+
* @param mu (Sequence of) location(s).
25+
* @param sigma (Sequence of) scale(s).
26+
* @return The log of the product of densities.
27+
*/
28+
template <
29+
typename T_y_cl, typename T_loc_cl, typename T_scale_cl,
30+
require_all_prim_or_rev_kernel_expression_t<T_y_cl, T_loc_cl,
31+
T_scale_cl>* = nullptr,
32+
require_any_not_stan_scalar_t<T_y_cl, T_loc_cl, T_scale_cl>* = nullptr>
33+
return_type_t<T_y_cl, T_loc_cl, T_scale_cl> double_exponential_lcdf(
34+
const T_y_cl& y, const T_loc_cl& mu, const T_scale_cl& sigma) {
35+
static const char* function = "double_exponential_lcdf(OpenCL)";
36+
using T_partials_return = partials_return_t<T_y_cl, T_loc_cl, T_scale_cl>;
37+
using std::isfinite;
38+
using std::isnan;
39+
40+
check_consistent_sizes(function, "Random variable", y, "Location parameter",
41+
mu, "Scale parameter", sigma);
42+
const size_t N = max_size(y, mu, sigma);
43+
if (N == 0) {
44+
return 0.0;
45+
}
46+
47+
const auto& y_col = as_column_vector_or_scalar(y);
48+
const auto& mu_col = as_column_vector_or_scalar(mu);
49+
const auto& sigma_col = as_column_vector_or_scalar(sigma);
50+
51+
const auto& y_val = value_of(y_col);
52+
const auto& mu_val = value_of(mu_col);
53+
const auto& sigma_val = value_of(sigma_col);
54+
55+
auto check_y_not_nan
56+
= check_cl(function, "Random variable", y_val, "not NaN");
57+
auto y_not_nan_expr = !isnan(y_val);
58+
auto check_mu_finite
59+
= check_cl(function, "Location parameter", mu_val, "finite");
60+
auto mu_finite_expr = isfinite(mu_val);
61+
auto check_sigma_positive_finite
62+
= check_cl(function, "Scale parameter", sigma_val, "positive finite");
63+
auto sigma_positive_finite_expr = 0 < sigma_val && isfinite(sigma_val);
64+
65+
auto sigma_inv = elt_divide(1.0, sigma_val);
66+
auto scaled_diff = elt_multiply(y_val - mu_val, sigma_inv);
67+
auto cond = y_val < mu_val;
68+
auto y_deriv = select(cond, sigma_inv,
69+
elt_divide(sigma_inv, 2.0 * exp(scaled_diff) - 1.0));
70+
auto cdf_log_expr = colwise_sum(
71+
select(cond, LOG_HALF + scaled_diff, log1m(0.5 * exp(-scaled_diff))));
72+
auto mu_deriv = -y_deriv;
73+
auto sigma_deriv = elt_multiply(mu_deriv, scaled_diff);
74+
75+
matrix_cl<double> cdf_log_cl;
76+
matrix_cl<double> mu_deriv_cl;
77+
matrix_cl<double> y_deriv_cl;
78+
matrix_cl<double> sigma_deriv_cl;
79+
80+
results(check_y_not_nan, check_mu_finite, check_sigma_positive_finite,
81+
cdf_log_cl, y_deriv_cl, mu_deriv_cl, sigma_deriv_cl)
82+
= expressions(y_not_nan_expr, mu_finite_expr, sigma_positive_finite_expr,
83+
cdf_log_expr, calc_if<!is_constant<T_y_cl>::value>(y_deriv),
84+
calc_if<!is_constant<T_loc_cl>::value>(mu_deriv),
85+
calc_if<!is_constant<T_scale_cl>::value>(sigma_deriv));
86+
87+
T_partials_return cdf_log = (from_matrix_cl(cdf_log_cl)).sum();
88+
89+
operands_and_partials<decltype(y_col), decltype(mu_col), decltype(sigma_col)>
90+
ops_partials(y_col, mu_col, sigma_col);
91+
92+
if (!is_constant<T_y_cl>::value) {
93+
ops_partials.edge1_.partials_ = std::move(y_deriv_cl);
94+
}
95+
if (!is_constant<T_loc_cl>::value) {
96+
ops_partials.edge2_.partials_ = std::move(mu_deriv_cl);
97+
}
98+
if (!is_constant<T_scale_cl>::value) {
99+
ops_partials.edge3_.partials_ = std::move(sigma_deriv_cl);
100+
}
101+
return ops_partials.build(cdf_log);
102+
}
103+
104+
} // namespace math
105+
} // namespace stan
106+
#endif
107+
#endif

0 commit comments

Comments
 (0)