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