1 #ifndef STAN_MATH_PRIM_FUN_UB_CONSTRAIN_HPP
2 #define STAN_MATH_PRIM_FUN_UB_CONSTRAIN_HPP
3 
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/constants.hpp>
7 #include <stan/math/prim/fun/exp.hpp>
8 #include <stan/math/prim/fun/identity_constrain.hpp>
9 #include <stan/math/prim/fun/identity_free.hpp>
10 #include <stan/math/prim/fun/subtract.hpp>
11 #include <stan/math/prim/fun/sum.hpp>
12 #include <stan/math/prim/fun/value_of.hpp>
13 #include <cmath>
14 
15 namespace stan {
16 namespace math {
17 
18 /**
19  * Return the upper-bounded value for the specified unconstrained
20  * matrix and upper bound.
21  *
22  * <p>The transform is
23  *
24  * <p>\f$f(x) = U - \exp(x)\f$
25  *
26  * <p>where \f$U\f$ is the upper bound.
27  *
28  * @tparam T type of Matrix
29  * @tparam U type of upper bound
30  * @param[in] x free Matrix.
31  * @param[in] ub upper bound
32  * @return matrix constrained to have upper bound
33  */
34 template <typename T, typename U, require_all_stan_scalar_t<T, U>* = nullptr,
35           require_all_not_st_var<T, U>* = nullptr>
ub_constrain(const T & x,const U & ub)36 inline auto ub_constrain(const T& x, const U& ub) {
37   if (value_of_rec(ub) == INFTY) {
38     return identity_constrain(x, ub);
39   } else {
40     return subtract(ub, exp(x));
41   }
42 }
43 
44 /**
45  * Return the upper-bounded value for the specified unconstrained
46  * scalar and upper bound and increment the specified log
47  * probability reference with the log absolute Jacobian
48  * determinant of the transform.
49  *
50  * <p>The transform is as specified for
51  * <code>ub_constrain(T, double)</code>.  The log absolute Jacobian
52  * determinant is
53  *
54  * <p>\f$ \log | \frac{d}{dx} -\mbox{exp}(x) + U |
55  *     = \log | -\mbox{exp}(x) + 0 | = x\f$.
56  *
57  * @tparam T type of scalar
58  * @tparam U type of upper bound
59  * @tparam S type of log probability
60  * @param[in] x free scalar
61  * @param[in] ub upper bound
62  * @param[in,out] lp log density
63  * @return scalar constrained to have upper bound
64  */
65 template <typename T, typename U, require_all_stan_scalar_t<T, U>* = nullptr,
66           require_all_not_st_var<T, U>* = nullptr>
ub_constrain(const T & x,const U & ub,std::decay_t<return_type_t<T,U>> & lp)67 inline auto ub_constrain(const T& x, const U& ub,
68                          std::decay_t<return_type_t<T, U>>& lp) {
69   if (value_of_rec(ub) == INFTY) {
70     return identity_constrain(x, ub);
71   } else {
72     lp += x;
73     return subtract(ub, exp(x));
74   }
75 }
76 
77 /**
78  * Specialization of `ub_constrain` to apply a scalar upper bound elementwise
79  *  to each input.
80  *
81  * @tparam T A type inheriting from `EigenBase`.
82  * @tparam U Scalar.
83  * @param[in] x unconstrained input
84  * @param[in] ub upper bound on output
85  * @return upper-bound constrained value corresponding to inputs
86  */
87 template <typename T, typename U, require_eigen_t<T>* = nullptr,
88           require_stan_scalar_t<U>* = nullptr,
89           require_all_not_st_var<T, U>* = nullptr>
ub_constrain(const T & x,const U & ub)90 inline auto ub_constrain(const T& x, const U& ub) {
91   return eval(x.unaryExpr([ub](auto&& xx) { return ub_constrain(xx, ub); }));
92 }
93 
94 /**
95  * Specialization of `ub_constrain` to apply a scalar upper bound elementwise
96  *  to each input.
97  *
98  * @tparam T A type inheriting from `EigenBase`.
99  * @tparam U Scalar.
100  * @param[in] x unconstrained input
101  * @param[in] ub upper bound on output
102  * @param[in,out] lp reference to log probability to increment
103  * @return upper-bound constrained value corresponding to inputs
104  */
105 template <typename T, typename U, require_eigen_t<T>* = nullptr,
106           require_stan_scalar_t<U>* = nullptr,
107           require_all_not_st_var<T, U>* = nullptr>
ub_constrain(const T & x,const U & ub,std::decay_t<return_type_t<T,U>> & lp)108 inline auto ub_constrain(const T& x, const U& ub,
109                          std::decay_t<return_type_t<T, U>>& lp) {
110   return eval(
111       x.unaryExpr([ub, &lp](auto&& xx) { return ub_constrain(xx, ub, lp); }));
112 }
113 
114 /**
115  * Specialization of `ub_constrain` to apply a matrix of upper bounds
116  * elementwise to each input element.
117  *
118  * @tparam T A type inheriting from `EigenBase`.
119  * @tparam U A type inheriting from `EigenBase`.
120  * @param[in] x unconstrained input
121  * @param[in] ub upper bound on output
122  * @return upper-bound constrained value corresponding to inputs
123  */
124 template <typename T, typename U, require_all_eigen_t<T, U>* = nullptr,
125           require_all_not_st_var<T, U>* = nullptr>
ub_constrain(const T & x,const U & ub)126 inline auto ub_constrain(const T& x, const U& ub) {
127   check_matching_dims("ub_constrain", "x", x, "ub", ub);
128   return eval(x.binaryExpr(
129       ub, [](auto&& xx, auto&& ubb) { return ub_constrain(xx, ubb); }));
130 }
131 
132 /**
133  * Specialization of `ub_constrain` to apply a matrix of upper bounds
134  * elementwise to each input element.
135  *
136  * @tparam T A type inheriting from `EigenBase`.
137  * @tparam U A type inheriting from `EigenBase`.
138  * @param[in] x unconstrained input
139  * @param[in] ub upper bound on output
140  * @param[in,out] lp reference to log probability to increment
141  * @return upper-bound constrained value corresponding to inputs
142  */
143 template <typename T, typename U, require_all_eigen_t<T, U>* = nullptr,
144           require_all_not_st_var<T, U>* = nullptr>
ub_constrain(const T & x,const U & ub,std::decay_t<return_type_t<T,U>> & lp)145 inline auto ub_constrain(const T& x, const U& ub,
146                          std::decay_t<return_type_t<T, U>>& lp) {
147   check_matching_dims("ub_constrain", "x", x, "ub", ub);
148   return eval(x.binaryExpr(
149       ub, [&lp](auto&& xx, auto&& ubb) { return ub_constrain(xx, ubb, lp); }));
150 }
151 
152 /**
153  * Specialization of `ub_constrain` to apply a scalar upper bound elementwise
154  *  to each input element.
155  *
156  * @tparam T A Any type with a Scalar `scalar_type`.
157  * @tparam U Scalar.
158  * @param[in] x unconstrained input
159  * @param[in] ub upper bound on output
160  * @return lower-bound constrained value corresponding to inputs
161  */
162 template <typename T, typename U, require_not_std_vector_t<U>* = nullptr>
ub_constrain(const std::vector<T> & x,const U & ub)163 inline auto ub_constrain(const std::vector<T>& x, const U& ub) {
164   std::vector<plain_type_t<decltype(ub_constrain(x[0], ub))>> ret(x.size());
165   for (size_t i = 0; i < x.size(); ++i) {
166     ret[i] = ub_constrain(x[i], ub);
167   }
168   return ret;
169 }
170 
171 /**
172  * Specialization of `ub_constrain` to apply a scalar upper bound elementwise
173  *  to each input element.
174  *
175  * @tparam T A Any type with a Scalar `scalar_type`.
176  * @tparam U Scalar.
177  * @param[in] x unconstrained input
178  * @param[in] ub upper bound on output
179  * @param[in,out] lp reference to log probability to increment
180  * @return lower-bound constrained value corresponding to inputs
181  */
182 template <typename T, typename U, require_not_std_vector_t<U>* = nullptr>
ub_constrain(const std::vector<T> & x,const U & ub,return_type_t<T,U> & lp)183 inline auto ub_constrain(const std::vector<T>& x, const U& ub,
184                          return_type_t<T, U>& lp) {
185   std::vector<plain_type_t<decltype(ub_constrain(x[0], ub))>> ret(x.size());
186   for (size_t i = 0; i < x.size(); ++i) {
187     ret[i] = ub_constrain(x[i], ub, lp);
188   }
189   return ret;
190 }
191 
192 /**
193  * Specialization of `ub_constrain` to apply a container of upper bounds
194  * elementwise to each input element.
195  *
196  * @tparam T A Any type with a Scalar `scalar_type`.
197  * @tparam U A type inheriting from `EigenBase` or a standard vector.
198  * @param[in] x unconstrained input
199  * @param[in] ub upper bound on output
200  * @return lower-bound constrained value corresponding to inputs
201  */
202 template <typename T, typename U>
ub_constrain(const std::vector<T> & x,const std::vector<U> & ub)203 inline auto ub_constrain(const std::vector<T>& x, const std::vector<U>& ub) {
204   check_matching_dims("ub_constrain", "x", x, "ub", ub);
205   std::vector<plain_type_t<decltype(ub_constrain(x[0], ub[0]))>> ret(x.size());
206   for (size_t i = 0; i < x.size(); ++i) {
207     ret[i] = ub_constrain(x[i], ub[i]);
208   }
209   return ret;
210 }
211 
212 /**
213  * Specialization of `ub_constrain` to apply a container of upper bounds
214  * elementwise to each input element.
215  *
216  * @tparam T A Any type with a Scalar `scalar_type`.
217  * @tparam U A type inheriting from `EigenBase` or a standard vector.
218  * @param[in] x unconstrained input
219  * @param[in] ub upper bound on output
220  * @param[in,out] lp reference to log probability to increment
221  * @return lower-bound constrained value corresponding to inputs
222  */
223 template <typename T, typename U>
ub_constrain(const std::vector<T> & x,const std::vector<U> & ub,return_type_t<T,U> & lp)224 inline auto ub_constrain(const std::vector<T>& x, const std::vector<U>& ub,
225                          return_type_t<T, U>& lp) {
226   check_matching_dims("ub_constrain", "x", x, "ub", ub);
227   std::vector<plain_type_t<decltype(ub_constrain(x[0], ub[0]))>> ret(x.size());
228   for (size_t i = 0; i < x.size(); ++i) {
229     ret[i] = ub_constrain(x[i], ub[i], lp);
230   }
231   return ret;
232 }
233 
234 /**
235  * Specialization of `ub_constrain` to apply a container of upper bounds
236  * elementwise to each input element. If the `Jacobian` parameter is `true`, the
237  * log density accumulator is incremented with the log absolute Jacobian
238  * determinant of the transform.  All of the transforms are specified with their
239  * Jacobians in the *Stan Reference Manual* chapter Constraint Transforms.
240  *
241  * @tparam Jacobian if `true`, increment log density accumulator with log
242  * absolute Jacobian determinant of constraining transform
243  * @tparam T A type inheriting from `Eigen::EigenBase`, a `var_value` with inner
244  * type inheriting from `Eigen::EigenBase`, a standard vector, or a scalar
245  * @tparam U A type inheriting from `Eigen::EigenBase`, a `var_value` with inner
246  * type inheriting from `Eigen::EigenBase`, a standard vector, or a scalar
247  * @param[in] x unconstrained input
248  * @param[in] ub upper bound on output
249  * @param[in, out] lp log density accumulator
250  * @return lower-bound constrained value corresponding to inputs
251  */
252 template <bool Jacobian, typename T, typename U>
ub_constrain(const T & x,const U & ub,return_type_t<T,U> & lp)253 inline auto ub_constrain(const T& x, const U& ub, return_type_t<T, U>& lp) {
254   if (Jacobian) {
255     return ub_constrain(x, ub, lp);
256   } else {
257     return ub_constrain(x, ub);
258   }
259 }
260 
261 }  // namespace math
262 }  // namespace stan
263 
264 #endif
265