Skip to content

Commit 616a597

Browse files
t4c1t4c1
committed
added OpenCL frechet (lc)cdf
1 parent 3ebfeb7 commit 616a597

10 files changed

Lines changed: 774 additions & 3 deletions

File tree

‎stan/math/opencl/prim.hpp‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -140,6 +140,9 @@
140140
#include <stan/math/opencl/prim/exponential_lccdf.hpp>
141141
#include <stan/math/opencl/prim/exponential_lcdf.hpp>
142142
#include <stan/math/opencl/prim/exponential_lpdf.hpp>
143+
#include <stan/math/opencl/prim/frechet_cdf.hpp>
144+
#include <stan/math/opencl/prim/frechet_lccdf.hpp>
145+
#include <stan/math/opencl/prim/frechet_lcdf.hpp>
143146
#include <stan/math/opencl/prim/frechet_lpdf.hpp>
144147
#include <stan/math/opencl/prim/gamma_lpdf.hpp>
145148
#include <stan/math/opencl/prim/gp_exp_quad_cov.hpp>

‎stan/math/opencl/prim/frechet_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_FRECHET_CDF_HPP
2+
#define STAN_MATH_OPENCL_PRIM_FRECHET_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 frechet cumulative distribution function for the given
18+
* location, and scale. If given containers of matching sizes
19+
* returns the product of probabilities.
20+
*
21+
* @tparam T_y_cl type of scalar outcome
22+
* @tparam T_shape_cl type of location
23+
* @tparam T_scale_cl type of scale
24+
* @param y (Sequence of) scalar(s).
25+
* @param alpha (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_shape_cl, typename T_scale_cl,
31+
require_all_prim_or_rev_kernel_expression_t<T_y_cl, T_shape_cl,
32+
T_scale_cl>* = nullptr,
33+
require_any_not_stan_scalar_t<T_y_cl, T_shape_cl, T_scale_cl>* = nullptr>
34+
return_type_t<T_y_cl, T_shape_cl, T_scale_cl> frechet_cdf(
35+
const T_y_cl& y, const T_shape_cl& alpha, const T_scale_cl& sigma) {
36+
static const char* function = "frechet_cdf(OpenCL)";
37+
using T_partials_return = partials_return_t<T_y_cl, T_shape_cl, T_scale_cl>;
38+
using std::isfinite;
39+
using std::isnan;
40+
41+
check_consistent_sizes(function, "Random variable", y, "Shape parameter",
42+
alpha, "Scale parameter", sigma);
43+
const size_t N = max_size(y, alpha, sigma);
44+
if (N == 0) {
45+
return 1.0;
46+
}
47+
48+
const auto& y_col = as_column_vector_or_scalar(y);
49+
const auto& alpha_col = as_column_vector_or_scalar(alpha);
50+
const auto& sigma_col = as_column_vector_or_scalar(sigma);
51+
52+
const auto& y_val = value_of(y_col);
53+
const auto& alpha_val = value_of(alpha_col);
54+
const auto& sigma_val = value_of(sigma_col);
55+
56+
auto check_y_positive
57+
= check_cl(function, "Random variable", y_val, "positive");
58+
auto y_positive = y_val>0;
59+
auto check_alpha_positive_finite
60+
= check_cl(function, "Shape parameter", alpha_val, "positive finite");
61+
auto alpha_positive_finite_expr = alpha_val>0 && isfinite(alpha_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 pow_n = pow(elt_divide(sigma_val, y_val), alpha_val);
67+
auto cdf_n = exp(-pow_n);
68+
auto cdf_expr = colwise_prod(cdf_n);
69+
70+
auto pow_n_alpha =elt_multiply(pow_n, alpha_val);
71+
auto y_deriv_tmp =elt_divide(pow_n_alpha, y_val);
72+
auto alpha_deriv_tmp = elt_multiply(pow_n, log(elt_divide(y_val, sigma_val)));
73+
auto sigma_deriv_tmp = elt_divide(pow_n_alpha, -sigma_val);
74+
75+
matrix_cl<double> cdf_cl;
76+
matrix_cl<double> y_deriv_cl;
77+
matrix_cl<double> alpha_deriv_cl;
78+
matrix_cl<double> sigma_deriv_cl;
79+
80+
results(check_y_positive, check_alpha_positive_finite, check_sigma_positive_finite,
81+
cdf_cl, y_deriv_cl, alpha_deriv_cl, sigma_deriv_cl)
82+
= expressions(
83+
y_positive, alpha_positive_finite_expr, sigma_positive_finite_expr,
84+
cdf_expr,
85+
calc_if<!is_constant<T_y_cl>::value>(y_deriv_tmp),
86+
calc_if<!is_constant<T_shape_cl>::value>(alpha_deriv_tmp),
87+
calc_if<!is_constant<T_scale_cl>::value>(sigma_deriv_tmp));
88+
89+
T_partials_return cdf = (from_matrix_cl(cdf_cl)).prod();
90+
91+
auto alpha_deriv = alpha_deriv_cl * cdf;
92+
auto y_deriv = y_deriv_cl * cdf;
93+
auto sigma_deriv = sigma_deriv_cl * cdf;
94+
95+
results(alpha_deriv_cl, y_deriv_cl, sigma_deriv_cl)
96+
= expressions(calc_if<!is_constant<T_shape_cl>::value>(alpha_deriv),
97+
calc_if<!is_constant<T_y_cl>::value>(y_deriv),
98+
calc_if<!is_constant<T_scale_cl>::value>(sigma_deriv));
99+
100+
operands_and_partials<decltype(y_col), decltype(alpha_col), decltype(sigma_col)>
101+
ops_partials(y_col, alpha_col, sigma_col);
102+
103+
if (!is_constant<T_y_cl>::value) {
104+
ops_partials.edge1_.partials_ = std::move(y_deriv_cl);
105+
}
106+
if (!is_constant<T_shape_cl>::value) {
107+
ops_partials.edge2_.partials_ = std::move(alpha_deriv_cl);
108+
}
109+
if (!is_constant<T_scale_cl>::value) {
110+
ops_partials.edge3_.partials_ = std::move(sigma_deriv_cl);
111+
}
112+
return ops_partials.build(cdf);
113+
}
114+
115+
} // namespace math
116+
} // namespace stan
117+
#endif
118+
#endif

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

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
#ifndef STAN_MATH_OPENCL_PRIM_FRECHET_LCCDF_HPP
2+
#define STAN_MATH_OPENCL_PRIM_FRECHET_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 frechet log complementary cumulative distribution function for
18+
* the given location, and scale. If given containers of matching sizes returns
19+
* the product of probabilities.
20+
*
21+
* @tparam T_y_cl type of scalar outcome
22+
* @tparam T_shape_cl type of location
23+
* @tparam T_scale_cl type of scale
24+
* @param y (Sequence of) scalar(s).
25+
* @param alpha (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_shape_cl, typename T_scale_cl,
31+
require_all_prim_or_rev_kernel_expression_t<T_y_cl, T_shape_cl,
32+
T_scale_cl>* = nullptr,
33+
require_any_not_stan_scalar_t<T_y_cl, T_shape_cl, T_scale_cl>* = nullptr>
34+
return_type_t<T_y_cl, T_shape_cl, T_scale_cl> frechet_lccdf(
35+
const T_y_cl& y, const T_shape_cl& alpha, const T_scale_cl& sigma) {
36+
static const char* function = "frechet_lccdf(OpenCL)";
37+
using T_partials_return = partials_return_t<T_y_cl, T_shape_cl, T_scale_cl>;
38+
using std::isfinite;
39+
using std::isnan;
40+
41+
check_consistent_sizes(function, "Random variable", y, "Shape parameter",
42+
alpha, "Scale parameter", sigma);
43+
const size_t N = max_size(y, alpha, sigma);
44+
if (N == 0) {
45+
return 1.0;
46+
}
47+
48+
const auto& y_col = as_column_vector_or_scalar(y);
49+
const auto& alpha_col = as_column_vector_or_scalar(alpha);
50+
const auto& sigma_col = as_column_vector_or_scalar(sigma);
51+
52+
const auto& y_val = value_of(y_col);
53+
const auto& alpha_val = value_of(alpha_col);
54+
const auto& sigma_val = value_of(sigma_col);
55+
56+
auto check_y_positive
57+
= check_cl(function, "Random variable", y_val, "positive");
58+
auto y_positive = y_val > 0;
59+
auto check_alpha_positive_finite
60+
= check_cl(function, "Shape parameter", alpha_val, "positive finite");
61+
auto alpha_positive_finite_expr = alpha_val > 0 && isfinite(alpha_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 pow_n = pow(elt_divide(sigma_val, y_val), alpha_val);
67+
auto exp_n = exp(-pow_n);
68+
auto lccdf_expr = colwise_sum(log1m(exp_n));
69+
70+
auto rep_deriv = elt_divide(pow_n, elt_divide(1.0, exp_n - 1.0));
71+
auto y_deriv = elt_multiply(elt_divide(-alpha_val, y_val), rep_deriv);
72+
auto alpha_deriv = elt_multiply(-log(elt_divide(y_val, sigma_val)), rep_deriv);
73+
auto sigma_deriv = elt_multiply(elt_divide(alpha_val, sigma_val), rep_deriv);
74+
75+
matrix_cl<double> lccdf_cl;
76+
matrix_cl<double> y_deriv_cl;
77+
matrix_cl<double> alpha_deriv_cl;
78+
matrix_cl<double> sigma_deriv_cl;
79+
80+
results(check_y_positive, check_alpha_positive_finite,
81+
check_sigma_positive_finite, lccdf_cl, y_deriv_cl, alpha_deriv_cl,
82+
sigma_deriv_cl)
83+
= expressions(y_positive, alpha_positive_finite_expr,
84+
sigma_positive_finite_expr, lccdf_expr,
85+
calc_if<!is_constant<T_y_cl>::value>(y_deriv),
86+
calc_if<!is_constant<T_shape_cl>::value>(alpha_deriv),
87+
calc_if<!is_constant<T_scale_cl>::value>(sigma_deriv));
88+
89+
T_partials_return lccdf = (from_matrix_cl(lccdf_cl)).sum();
90+
91+
operands_and_partials<decltype(y_col), decltype(alpha_col),
92+
decltype(sigma_col)>
93+
ops_partials(y_col, alpha_col, sigma_col);
94+
95+
if (!is_constant<T_y_cl>::value) {
96+
ops_partials.edge1_.partials_ = std::move(y_deriv_cl);
97+
}
98+
if (!is_constant<T_shape_cl>::value) {
99+
ops_partials.edge2_.partials_ = std::move(alpha_deriv_cl);
100+
}
101+
if (!is_constant<T_scale_cl>::value) {
102+
ops_partials.edge3_.partials_ = std::move(sigma_deriv_cl);
103+
}
104+
return ops_partials.build(lccdf);
105+
}
106+
107+
} // namespace math
108+
} // namespace stan
109+
#endif
110+
#endif

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

Lines changed: 108 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,108 @@
1+
#ifndef STAN_MATH_OPENCL_PRIM_FRECHET_LCDF_HPP
2+
#define STAN_MATH_OPENCL_PRIM_FRECHET_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 frechet log cumulative distribution function for
18+
* the given location, and scale. If given containers of matching sizes returns
19+
* the product of probabilities.
20+
*
21+
* @tparam T_y_cl type of scalar outcome
22+
* @tparam T_shape_cl type of location
23+
* @tparam T_scale_cl type of scale
24+
* @param y (Sequence of) scalar(s).
25+
* @param alpha (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_shape_cl, typename T_scale_cl,
31+
require_all_prim_or_rev_kernel_expression_t<T_y_cl, T_shape_cl,
32+
T_scale_cl>* = nullptr,
33+
require_any_not_stan_scalar_t<T_y_cl, T_shape_cl, T_scale_cl>* = nullptr>
34+
return_type_t<T_y_cl, T_shape_cl, T_scale_cl> frechet_lcdf(
35+
const T_y_cl& y, const T_shape_cl& alpha, const T_scale_cl& sigma) {
36+
static const char* function = "frechet_lcdf(OpenCL)";
37+
using T_partials_return = partials_return_t<T_y_cl, T_shape_cl, T_scale_cl>;
38+
using std::isfinite;
39+
using std::isnan;
40+
41+
check_consistent_sizes(function, "Random variable", y, "Shape parameter",
42+
alpha, "Scale parameter", sigma);
43+
const size_t N = max_size(y, alpha, sigma);
44+
if (N == 0) {
45+
return 1.0;
46+
}
47+
48+
const auto& y_col = as_column_vector_or_scalar(y);
49+
const auto& alpha_col = as_column_vector_or_scalar(alpha);
50+
const auto& sigma_col = as_column_vector_or_scalar(sigma);
51+
52+
const auto& y_val = value_of(y_col);
53+
const auto& alpha_val = value_of(alpha_col);
54+
const auto& sigma_val = value_of(sigma_col);
55+
56+
auto check_y_positive
57+
= check_cl(function, "Random variable", y_val, "positive");
58+
auto y_positive = y_val > 0;
59+
auto check_alpha_positive_finite
60+
= check_cl(function, "Shape parameter", alpha_val, "positive finite");
61+
auto alpha_positive_finite_expr = alpha_val > 0 && isfinite(alpha_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 pow_n = pow(elt_divide(sigma_val, y_val), alpha_val);
67+
auto lcdf_expr = colwise_sum(pow_n);
68+
69+
auto y_deriv = elt_multiply(elt_divide(alpha_val, y_val), pow_n);
70+
auto alpha_deriv = elt_multiply(log(elt_divide(y_val, sigma_val)), pow_n);
71+
auto sigma_deriv = elt_multiply(elt_divide(alpha_val, -sigma_val), pow_n);
72+
73+
matrix_cl<double> lcdf_cl;
74+
matrix_cl<double> y_deriv_cl;
75+
matrix_cl<double> alpha_deriv_cl;
76+
matrix_cl<double> sigma_deriv_cl;
77+
78+
results(check_y_positive, check_alpha_positive_finite,
79+
check_sigma_positive_finite, lcdf_cl, y_deriv_cl, alpha_deriv_cl,
80+
sigma_deriv_cl)
81+
= expressions(y_positive, alpha_positive_finite_expr,
82+
sigma_positive_finite_expr, lcdf_expr,
83+
calc_if<!is_constant<T_y_cl>::value>(y_deriv),
84+
calc_if<!is_constant<T_shape_cl>::value>(alpha_deriv),
85+
calc_if<!is_constant<T_scale_cl>::value>(sigma_deriv));
86+
87+
T_partials_return lcdf = -(from_matrix_cl(lcdf_cl)).sum();
88+
89+
operands_and_partials<decltype(y_col), decltype(alpha_col),
90+
decltype(sigma_col)>
91+
ops_partials(y_col, alpha_col, sigma_col);
92+
93+
if (!is_constant<T_y_cl>::value) {
94+
ops_partials.edge1_.partials_ = std::move(y_deriv_cl);
95+
}
96+
if (!is_constant<T_shape_cl>::value) {
97+
ops_partials.edge2_.partials_ = std::move(alpha_deriv_cl);
98+
}
99+
if (!is_constant<T_scale_cl>::value) {
100+
ops_partials.edge3_.partials_ = std::move(sigma_deriv_cl);
101+
}
102+
return ops_partials.build(lcdf);
103+
}
104+
105+
} // namespace math
106+
} // namespace stan
107+
#endif
108+
#endif

‎stan/math/prim/prob/frechet_cdf.hpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,9 @@
2222
namespace stan {
2323
namespace math {
2424

25-
template <typename T_y, typename T_shape, typename T_scale>
25+
template <typename T_y, typename T_shape, typename T_scale,
26+
require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
27+
T_y, T_shape, T_scale>* = nullptr>
2628
return_type_t<T_y, T_shape, T_scale> frechet_cdf(const T_y& y,
2729
const T_shape& alpha,
2830
const T_scale& sigma) {

‎stan/math/prim/prob/frechet_lccdf.hpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,9 @@
2222
namespace stan {
2323
namespace math {
2424

25-
template <typename T_y, typename T_shape, typename T_scale>
25+
template <typename T_y, typename T_shape, typename T_scale,
26+
require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
27+
T_y, T_shape, T_scale>* = nullptr>
2628
return_type_t<T_y, T_shape, T_scale> frechet_lccdf(const T_y& y,
2729
const T_shape& alpha,
2830
const T_scale& sigma) {

‎stan/math/prim/prob/frechet_lcdf.hpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,9 @@
2020
namespace stan {
2121
namespace math {
2222

23-
template <typename T_y, typename T_shape, typename T_scale>
23+
template <typename T_y, typename T_shape, typename T_scale,
24+
require_all_not_nonscalar_prim_or_rev_kernel_expression_t<
25+
T_y, T_shape, T_scale>* = nullptr>
2426
return_type_t<T_y, T_shape, T_scale> frechet_lcdf(const T_y& y,
2527
const T_shape& alpha,
2628
const T_scale& sigma) {

0 commit comments

Comments
 (0)