1 #ifndef STAN_MATH_PRIM_FUN_LUB_CONSTRAIN_HPP
2 #define STAN_MATH_PRIM_FUN_LUB_CONSTRAIN_HPP
3
4 #include <stan/math/prim/fun/Eigen.hpp>
5 #include <stan/math/prim/meta.hpp>
6 #include <stan/math/prim/err.hpp>
7 #include <stan/math/prim/fun/add.hpp>
8 #include <stan/math/prim/fun/exp.hpp>
9 #include <stan/math/prim/fun/elt_multiply.hpp>
10 #include <stan/math/prim/fun/identity_constrain.hpp>
11 #include <stan/math/prim/fun/inv_logit.hpp>
12 #include <stan/math/prim/fun/log.hpp>
13 #include <stan/math/prim/fun/log1p.hpp>
14 #include <stan/math/prim/fun/lb_constrain.hpp>
15 #include <stan/math/prim/fun/multiply.hpp>
16 #include <stan/math/prim/fun/subtract.hpp>
17 #include <stan/math/prim/fun/sum.hpp>
18 #include <stan/math/prim/fun/ub_constrain.hpp>
19 #include <cmath>
20
21 namespace stan {
22 namespace math {
23
24 /**
25 * Return the lower and upper-bounded scalar derived by
26 * transforming the specified free scalar given the specified
27 * lower and upper bounds.
28 *
29 * <p>The transform is the transformed and scaled inverse logit,
30 *
31 * <p>\f$f(x) = L + (U - L) \mbox{logit}^{-1}(x)\f$
32 *
33 * @tparam T Scalar.
34 * @tparam L Scalar.
35 * @tparam U Scalar.
36 * @param[in] x Free scalar to transform.
37 * @param[in] lb Lower bound.
38 * @param[in] ub Upper bound.
39 * @return Lower- and upper-bounded scalar derived from transforming
40 * the free scalar.
41 * @throw std::domain_error if ub <= lb
42 */
43 template <typename T, typename L, typename U,
44 require_all_stan_scalar_t<T, L, U>* = nullptr,
45 require_not_var_t<return_type_t<T, L, U>>* = nullptr>
lub_constrain(T && x,L && lb,U && ub)46 inline auto lub_constrain(T&& x, L&& lb, U&& ub) {
47 const bool is_lb_inf = value_of(lb) == NEGATIVE_INFTY;
48 const bool is_ub_inf = value_of(ub) == INFTY;
49 if (unlikely(is_ub_inf && is_lb_inf)) {
50 return identity_constrain(x, lb, ub);
51 } else if (unlikely(is_ub_inf)) {
52 return lb_constrain(identity_constrain(x, ub), lb);
53 } else if (unlikely(is_lb_inf)) {
54 return ub_constrain(identity_constrain(x, lb), ub);
55 } else {
56 check_less("lub_constrain", "lb", value_of(lb), value_of(ub));
57 return (ub - lb) * inv_logit(x) + lb;
58 }
59 }
60
61 /**
62 * Return the lower- and upper-bounded scalar derived by
63 * transforming the specified free scalar given the specified
64 * lower and upper bounds and increment the specified log
65 * density with the log absolute Jacobian determinant.
66 *
67 * <p>The transform is as defined in
68 * <code>lub_constrain(T, double, double)</code>. The log absolute
69 * Jacobian determinant is given by
70 *
71 * <p>\f$\log \left| \frac{d}{dx} \left(
72 * L + (U-L) \mbox{logit}^{-1}(x) \right)
73 * \right|\f$
74 *
75 * <p>\f$ {} = \log |
76 * (U-L)
77 * \, (\mbox{logit}^{-1}(x))
78 * \, (1 - \mbox{logit}^{-1}(x)) |\f$
79 *
80 * <p>\f$ {} = \log (U - L) + \log (\mbox{logit}^{-1}(x))
81 * + \log (1 - \mbox{logit}^{-1}(x))\f$
82 *
83 * @tparam T Scalar.
84 * @tparam L Scalar.
85 * @tparam U Scalar.
86 * @param[in] x Free scalar to transform.
87 * @param[in] lb Lower bound.
88 * @param[in] ub Upper bound.
89 * @param[in,out] lp Log probability scalar reference.
90 * @return Lower- and upper-bounded scalar derived from transforming
91 * the free scalar.
92 * @throw std::domain_error if ub <= lb
93 */
94 template <typename T, typename L, typename U,
95 require_all_stan_scalar_t<T, L, U>* = nullptr,
96 require_not_var_t<return_type_t<T, L, U>>* = nullptr>
lub_constrain(T && x,L && lb,U && ub,return_type_t<T,L,U> & lp)97 inline auto lub_constrain(T&& x, L&& lb, U&& ub, return_type_t<T, L, U>& lp) {
98 const bool is_lb_inf = value_of(lb) == NEGATIVE_INFTY;
99 const bool is_ub_inf = value_of(ub) == INFTY;
100 if (unlikely(is_ub_inf && is_lb_inf)) {
101 return identity_constrain(x, ub, lb);
102 } else if (unlikely(is_ub_inf)) {
103 return lb_constrain(identity_constrain(x, ub), lb, lp);
104 } else if (unlikely(is_lb_inf)) {
105 return ub_constrain(identity_constrain(x, lb), ub, lp);
106 } else {
107 check_less("lub_constrain", "lb", value_of(lb), value_of(ub));
108 const auto diff = ub - lb;
109 lp += add(log(diff), subtract(-abs(x), multiply(2.0, log1p_exp(-abs(x)))));
110 return diff * inv_logit(x) + lb;
111 }
112 }
113
114 /**
115 * Overload for Eigen matrix and scalar bounds.
116 */
117 template <typename T, typename L, typename U, require_eigen_t<T>* = nullptr,
118 require_all_stan_scalar_t<L, U>* = nullptr,
119 require_not_var_t<return_type_t<T, L, U>>* = nullptr>
lub_constrain(const T & x,const L & lb,const U & ub)120 inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
121 return eval(
122 x.unaryExpr([ub, lb](auto&& xx) { return lub_constrain(xx, lb, ub); }));
123 }
124
125 /**
126 * Overload for Eigen matrix and scalar bounds plus lp.
127 */
128 template <typename T, typename L, typename U, require_eigen_t<T>* = nullptr,
129 require_all_stan_scalar_t<L, U>* = nullptr,
130 require_not_var_t<return_type_t<T, L, U>>* = nullptr>
lub_constrain(const T & x,const L & lb,const U & ub,return_type_t<T,L,U> & lp)131 inline auto lub_constrain(const T& x, const L& lb, const U& ub,
132 return_type_t<T, L, U>& lp) {
133 return eval(x.unaryExpr(
134 [lb, ub, &lp](auto&& xx) { return lub_constrain(xx, lb, ub, lp); }));
135 }
136
137 /**
138 * Overload for Eigen matrix with matrix lower bound and scalar upper
139 * bound.
140 */
141 template <typename T, typename L, typename U,
142 require_all_eigen_t<T, L>* = nullptr,
143 require_stan_scalar_t<U>* = nullptr,
144 require_not_var_t<return_type_t<T, L, U>>* = nullptr>
lub_constrain(const T & x,const L & lb,const U & ub)145 inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
146 check_matching_dims("lub_constrain", "x", x, "lb", lb);
147 return eval(x.binaryExpr(
148 lb, [ub](auto&& x, auto&& lb) { return lub_constrain(x, lb, ub); }));
149 }
150
151 /**
152 * Overload for Eigen matrix with matrix lower bound and scalar upper
153 * bound plus lp.
154 */
155 template <typename T, typename L, typename U,
156 require_all_eigen_t<T, L>* = nullptr,
157 require_stan_scalar_t<U>* = nullptr,
158 require_not_var_t<return_type_t<T, L, U>>* = nullptr>
lub_constrain(const T & x,const L & lb,const U & ub,return_type_t<T,L,U> & lp)159 inline auto lub_constrain(const T& x, const L& lb, const U& ub,
160 return_type_t<T, L, U>& lp) {
161 check_matching_dims("lub_constrain", "x", x, "lb", lb);
162 return eval(x.binaryExpr(lb, [ub, &lp](auto&& x, auto&& lb) {
163 return lub_constrain(x, lb, ub, lp);
164 }));
165 }
166
167 /**
168 * Overload for Eigen matrix with scalar lower bound and matrix upper
169 * bound.
170 */
171 template <typename T, typename L, typename U,
172 require_all_eigen_t<T, U>* = nullptr,
173 require_stan_scalar_t<L>* = nullptr,
174 require_not_var_t<return_type_t<T, L, U>>* = nullptr>
lub_constrain(const T & x,const L & lb,const U & ub)175 inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
176 check_matching_dims("lub_constrain", "x", x, "ub", ub);
177 return eval(x.binaryExpr(
178 ub, [lb](auto&& x, auto&& ub) { return lub_constrain(x, lb, ub); }));
179 }
180
181 /**
182 * Overload for Eigen matrix with scalar lower bound and matrix upper
183 * bound plus lp.
184 */
185 template <typename T, typename L, typename U,
186 require_all_eigen_t<T, U>* = nullptr,
187 require_stan_scalar_t<L>* = nullptr,
188 require_not_var_t<return_type_t<T, L, U>>* = nullptr>
lub_constrain(const T & x,const L & lb,const U & ub,return_type_t<T,L,U> & lp)189 inline auto lub_constrain(const T& x, const L& lb, const U& ub,
190 return_type_t<T, L, U>& lp) {
191 check_matching_dims("lub_constrain", "x", x, "ub", ub);
192 return eval(x.binaryExpr(ub, [lb, &lp](auto&& x, auto&& ub) {
193 return lub_constrain(x, lb, ub, lp);
194 }));
195 }
196
197 /**
198 * Overload for Eigen matrix and matrix bounds.
199 */
200 template <typename T, typename L, typename U,
201 require_all_eigen_t<T, L, U>* = nullptr,
202 require_not_var_t<return_type_t<T, L, U>>* = nullptr>
lub_constrain(const T & x,const L & lb,const U & ub)203 inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
204 check_matching_dims("lub_constrain", "x", x, "lb", lb);
205 check_matching_dims("lub_constrain", "x", x, "ub", ub);
206 auto x_ref = to_ref(x);
207 auto lb_ref = to_ref(lb);
208 auto ub_ref = to_ref(ub);
209 promote_scalar_t<return_type_t<T, L, U>, T> x_ret(x.rows(), x.cols());
210 for (Eigen::Index j = 0; j < x_ref.cols(); ++j) {
211 for (Eigen::Index i = 0; i < x_ref.rows(); ++i) {
212 x_ret.coeffRef(i, j) = lub_constrain(
213 x_ref.coeff(i, j), lb_ref.coeff(i, j), ub_ref.coeff(i, j));
214 }
215 }
216 return x_ret;
217 }
218
219 /**
220 * Overload for Eigen matrix and matrix bounds plus lp.
221 */
222 template <typename T, typename L, typename U,
223 require_all_eigen_t<T, L, U>* = nullptr,
224 require_not_var_t<return_type_t<T, L, U>>* = nullptr>
lub_constrain(const T & x,const L & lb,const U & ub,return_type_t<T,L,U> & lp)225 inline auto lub_constrain(const T& x, const L& lb, const U& ub,
226 return_type_t<T, L, U>& lp) {
227 check_matching_dims("lub_constrain", "x", x, "lb", lb);
228 check_matching_dims("lub_constrain", "x", x, "ub", ub);
229 auto x_ref = to_ref(x);
230 auto lb_ref = to_ref(lb);
231 auto ub_ref = to_ref(ub);
232 promote_scalar_t<return_type_t<T, L, U>, T> x_ret(x.rows(), x.cols());
233 for (Eigen::Index j = 0; j < x_ref.cols(); ++j) {
234 for (Eigen::Index i = 0; i < x_ref.rows(); ++i) {
235 x_ret.coeffRef(i, j) = lub_constrain(
236 x_ref.coeff(i, j), lb_ref.coeff(i, j), ub_ref.coeff(i, j), lp);
237 }
238 }
239 return x_ret;
240 }
241
242 /**
243 * Overload for array of x and non-array lb and ub
244 */
245 template <typename T, typename L, typename U,
246 require_all_not_std_vector_t<L, U>* = nullptr>
lub_constrain(const std::vector<T> & x,const L & lb,const U & ub)247 inline auto lub_constrain(const std::vector<T>& x, const L& lb, const U& ub) {
248 std::vector<plain_type_t<decltype(lub_constrain(x[0], lb, ub))>> ret(
249 x.size());
250 for (size_t i = 0; i < x.size(); ++i) {
251 ret[i] = lub_constrain(x[i], lb, ub);
252 }
253 return ret;
254 }
255
256 /**
257 * Overload for array of x and non-array lb and ub with lp
258 */
259 template <typename T, typename L, typename U,
260 require_all_not_std_vector_t<L, U>* = nullptr>
lub_constrain(const std::vector<T> & x,const L & lb,const U & ub,return_type_t<T,L,U> & lp)261 inline auto lub_constrain(const std::vector<T>& x, const L& lb, const U& ub,
262 return_type_t<T, L, U>& lp) {
263 std::vector<plain_type_t<decltype(lub_constrain(x[0], lb, ub))>> ret(
264 x.size());
265 for (size_t i = 0; i < x.size(); ++i) {
266 ret[i] = lub_constrain(x[i], lb, ub, lp);
267 }
268 return ret;
269 }
270
271 /**
272 * Overload for array of x and ub and non-array lb
273 */
274 template <typename T, typename L, typename U,
275 require_not_std_vector_t<L>* = nullptr>
lub_constrain(const std::vector<T> & x,const L & lb,const std::vector<U> & ub)276 inline auto lub_constrain(const std::vector<T>& x, const L& lb,
277 const std::vector<U>& ub) {
278 check_matching_dims("lub_constrain", "x", x, "ub", ub);
279 std::vector<plain_type_t<decltype(lub_constrain(x[0], lb, ub[0]))>> ret(
280 x.size());
281 for (size_t i = 0; i < x.size(); ++i) {
282 ret[i] = lub_constrain(x[i], lb, ub[i]);
283 }
284 return ret;
285 }
286
287 /**
288 * Overload for array of x and ub and non-array lb with lp
289 */
290 template <typename T, typename L, typename U,
291 require_not_std_vector_t<L>* = nullptr>
lub_constrain(const std::vector<T> & x,const L & lb,const std::vector<U> & ub,return_type_t<T,L,U> & lp)292 inline auto lub_constrain(const std::vector<T>& x, const L& lb,
293 const std::vector<U>& ub,
294 return_type_t<T, L, U>& lp) {
295 check_matching_dims("lub_constrain", "x", x, "ub", ub);
296 std::vector<plain_type_t<decltype(lub_constrain(x[0], lb, ub[0]))>> ret(
297 x.size());
298 for (size_t i = 0; i < x.size(); ++i) {
299 ret[i] = lub_constrain(x[i], lb, ub[i], lp);
300 }
301 return ret;
302 }
303
304 /**
305 * Overload for array of x and lb and non-array ub
306 */
307 template <typename T, typename L, typename U,
308 require_not_std_vector_t<U>* = nullptr>
lub_constrain(const std::vector<T> & x,const std::vector<L> & lb,const U & ub)309 inline auto lub_constrain(const std::vector<T>& x, const std::vector<L>& lb,
310 const U& ub) {
311 check_matching_dims("lub_constrain", "x", x, "lb", lb);
312 std::vector<plain_type_t<decltype(lub_constrain(x[0], lb[0], ub))>> ret(
313 x.size());
314 for (size_t i = 0; i < x.size(); ++i) {
315 ret[i] = lub_constrain(x[i], lb[i], ub);
316 }
317 return ret;
318 }
319
320 /**
321 * Overload for array of x and lb and non-array ub with lp
322 */
323 template <typename T, typename L, typename U,
324 require_not_std_vector_t<U>* = nullptr>
lub_constrain(const std::vector<T> & x,const std::vector<L> & lb,const U & ub,return_type_t<T,L,U> & lp)325 inline auto lub_constrain(const std::vector<T>& x, const std::vector<L>& lb,
326 const U& ub, return_type_t<T, L, U>& lp) {
327 check_matching_dims("lub_constrain", "x", x, "lb", lb);
328 std::vector<plain_type_t<decltype(lub_constrain(x[0], lb[0], ub))>> ret(
329 x.size());
330 for (size_t i = 0; i < x.size(); ++i) {
331 ret[i] = lub_constrain(x[i], lb[i], ub, lp);
332 }
333 return ret;
334 }
335
336 /**
337 * Overload for array of x, lb, and ub with lp
338 */
339 template <typename T, typename L, typename U>
lub_constrain(const std::vector<T> & x,const std::vector<L> & lb,const std::vector<U> & ub)340 inline auto lub_constrain(const std::vector<T>& x, const std::vector<L>& lb,
341 const std::vector<U>& ub) {
342 check_matching_dims("lub_constrain", "x", x, "lb", lb);
343 check_matching_dims("lub_constrain", "x", x, "ub", ub);
344 std::vector<plain_type_t<decltype(lub_constrain(x[0], lb[0], ub[0]))>> ret(
345 x.size());
346 for (size_t i = 0; i < x.size(); ++i) {
347 ret[i] = lub_constrain(x[i], lb[i], ub[i]);
348 }
349 return ret;
350 }
351
352 /**
353 * Overload for array of x, lb, and ub
354 */
355 template <typename T, typename L, typename U>
lub_constrain(const std::vector<T> & x,const std::vector<L> & lb,const std::vector<U> & ub,return_type_t<T,L,U> & lp)356 inline auto lub_constrain(const std::vector<T>& x, const std::vector<L>& lb,
357 const std::vector<U>& ub,
358 return_type_t<T, L, U>& lp) {
359 check_matching_dims("lub_constrain", "x", x, "lb", lb);
360 check_matching_dims("lub_constrain", "x", x, "ub", ub);
361 std::vector<plain_type_t<decltype(lub_constrain(x[0], lb[0], ub[0]))>> ret(
362 x.size());
363 for (size_t i = 0; i < x.size(); ++i) {
364 ret[i] = lub_constrain(x[i], lb[i], ub[i], lp);
365 }
366 return ret;
367 }
368
369 /**
370 * Return the lower and upper-bounded scalar derived by
371 * transforming the specified free scalar given the specified
372 * lower and upper bounds. If the `Jacobian` parameter is `true`, the log
373 * density accumulator is incremented with the log absolute Jacobian determinant
374 * of the transform. All of the transforms are specified with their Jacobians
375 * in the *Stan Reference Manual* chapter Constraint Transforms.
376 *
377 * @tparam Jacobian if `true`, increment log density accumulator with log
378 * absolute Jacobian determinant of constraining transform
379 * @tparam T A type inheriting from `Eigen::EigenBase`, a `var_value` with inner
380 * type inheriting from `Eigen::EigenBase`, a standard vector, or a scalar
381 * @tparam L A type inheriting from `Eigen::EigenBase`, a `var_value` with inner
382 * type inheriting from `Eigen::EigenBase`, a standard vector, or a scalar
383 * @tparam U A type inheriting from `Eigen::EigenBase`, a `var_value` with inner
384 * type inheriting from `Eigen::EigenBase`, a standard vector, or a scalar
385 * @param[in] x Free scalar to transform
386 * @param[in] lb Lower bound
387 * @param[in] ub Upper bound
388 * @param[in, out] lp log density accumulator
389 * @return Lower- and upper-bounded scalar derived from transforming the free
390 * scalar
391 * @throw std::domain_error if `ub <= lb`
392 */
393 template <bool Jacobian, typename T, typename L, typename U>
lub_constrain(const T & x,const L & lb,const U & ub,return_type_t<T,L,U> & lp)394 inline auto lub_constrain(const T& x, const L& lb, const U& ub,
395 return_type_t<T, L, U>& lp) {
396 if (Jacobian) {
397 return lub_constrain(x, lb, ub, lp);
398 } else {
399 return lub_constrain(x, lb, ub);
400 }
401 }
402
403 } // namespace math
404 } // namespace stan
405
406 #endif
407