1 #ifndef STAN_MATH_PRIM_FUN_LUB_FREE_HPP
2 #define STAN_MATH_PRIM_FUN_LUB_FREE_HPP
3 
4 #include <stan/math/prim/meta.hpp>
5 #include <stan/math/prim/err.hpp>
6 #include <stan/math/prim/fun/divide.hpp>
7 #include <stan/math/prim/fun/eval.hpp>
8 #include <stan/math/prim/fun/logit.hpp>
9 #include <stan/math/prim/fun/lb_free.hpp>
10 #include <stan/math/prim/fun/ub_free.hpp>
11 #include <stan/math/prim/fun/subtract.hpp>
12 
13 namespace stan {
14 namespace math {
15 
16 /** lub_free
17  * Return the unconstrained variable that transforms to the
18  * y given the specified bounds.
19  *
20  * <p>The transform in `lub_constrain`,
21  * is reversed by a transformed and scaled logit,
22  *
23  * <p>\f$f^{-1}(y) = \mbox{logit}(\frac{y - L}{U - L})\f$
24  *
25  * where \f$U\f$ and \f$L\f$ are the lower and upper bounds.
26  *
27  * @tparam T type of bounded object
28  * @tparam L type of lower bound
29  * @tparam U type of upper bound
30  * @param y constrained value
31  * @param lb lower bound
32  * @param ub upper bound
33  * @return the free object that transforms to the input
34  *   given the bounds
35  * @throw std::invalid_argument if the lower bound is greater than
36  *   the upper bound, y is less than the lower bound, or y is
37  *   greater than the upper bound
38  */
39 ///@{
40 /**
41  * Overload for all scalar arguments
42  */
43 template <typename T, typename L, typename U,
44           require_not_std_vector_t<T>* = nullptr,
45           require_all_stan_scalar_t<L, U>* = nullptr>
lub_free(T && y,L && lb,U && ub)46 inline auto lub_free(T&& y, 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_free(y, lb, ub);
51   } else if (unlikely(is_ub_inf)) {
52     return lb_free(identity_free(y, ub), lb);
53   } else if (unlikely(is_lb_inf)) {
54     return ub_free(identity_free(y, lb), ub);
55   } else {
56     auto&& y_ref = to_ref(std::forward<T>(y));
57     check_bounded("lub_free", "Bounded variable", value_of(y_ref), value_of(lb),
58                   value_of(ub));
59     return eval(logit(divide(subtract(std::forward<decltype(y_ref)>(y_ref), lb),
60                              subtract(ub, lb))));
61   }
62 }
63 
64 /**
65  * Overload for matrix constrained variable, matrix lower bound, scalar upper
66  * bound
67  */
68 template <typename T, typename L, typename U,
69           require_all_eigen_t<T, L>* = nullptr,
70           require_stan_scalar_t<U>* = nullptr>
lub_free(T && y,L && lb,U && ub)71 inline auto lub_free(T&& y, L&& lb, U&& ub) {
72   check_matching_dims("lub_free", "y", y, "lb", lb);
73   auto&& y_ref = to_ref(std::forward<T>(y));
74   auto&& lb_ref = to_ref(std::forward<L>(lb));
75   promote_scalar_t<return_type_t<T, L, U>, T> ret(y.rows(), y.cols());
76   for (Eigen::Index j = 0; j < y.cols(); ++j) {
77     for (Eigen::Index i = 0; i < y.rows(); ++i) {
78       ret.coeffRef(i, j) = lub_free(y_ref.coeff(i, j), lb_ref.coeff(i, j), ub);
79     }
80   }
81   return ret;
82 }
83 
84 /**
85  * Overload for matrix constrained variable, matrix upper bound, scalar lower
86  * bound
87  */
88 template <typename T, typename L, typename U,
89           require_all_eigen_t<T, U>* = nullptr,
90           require_stan_scalar_t<L>* = nullptr>
lub_free(T && y,L && lb,U && ub)91 inline auto lub_free(T&& y, L&& lb, U&& ub) {
92   check_matching_dims("lub_free", "y", y, "ub", ub);
93   auto&& y_ref = to_ref(std::forward<T>(y));
94   auto&& ub_ref = to_ref(std::forward<U>(ub));
95   promote_scalar_t<return_type_t<T, L, U>, T> ret(y.rows(), y.cols());
96   for (Eigen::Index j = 0; j < y.cols(); ++j) {
97     for (Eigen::Index i = 0; i < y.rows(); ++i) {
98       ret.coeffRef(i, j) = lub_free(y_ref.coeff(i, j), lb, ub_ref.coeff(i, j));
99     }
100   }
101   return ret;
102 }
103 
104 /**
105  * Overload for matrix constrained variable, matrix upper bound, matrix lower
106  * bound
107  */
108 template <typename T, typename L, typename U,
109           require_all_eigen_t<T, L, U>* = nullptr>
lub_free(T && y,L && lb,U && ub)110 inline auto lub_free(T&& y, L&& lb, U&& ub) {
111   check_matching_dims("lub_free", "y", y, "lb", lb);
112   check_matching_dims("lub_free", "y", y, "ub", ub);
113   auto&& y_ref = to_ref(std::forward<T>(y));
114   auto&& lb_ref = to_ref(std::forward<L>(lb));
115   auto&& ub_ref = to_ref(std::forward<U>(ub));
116   promote_scalar_t<return_type_t<T, L, U>, T> ret(y.rows(), y.cols());
117   for (Eigen::Index j = 0; j < y.cols(); ++j) {
118     for (Eigen::Index i = 0; i < y.rows(); ++i) {
119       ret.coeffRef(i, j)
120           = lub_free(y_ref.coeff(i, j), lb_ref.coeff(i, j), ub_ref.coeff(i, j));
121     }
122   }
123   return ret;
124 }
125 
126 /**
127  * Overload for `std::vector` constrained variable
128  */
129 template <typename T, typename L, typename U,
130           require_all_not_std_vector_t<L, U>* = nullptr>
lub_free(const std::vector<T> y,const L & lb,const U & ub)131 inline auto lub_free(const std::vector<T> y, const L& lb, const U& ub) {
132   std::vector<decltype(lub_free(y[0], lb, ub))> ret(y.size());
133   for (Eigen::Index i = 0; i < y.size(); ++i) {
134     ret[i] = lub_free(y[i], lb, ub);
135   }
136   return ret;
137 }
138 
139 /**
140  * Overload for `std::vector` constrained variable and `std::vector` upper bound
141  */
142 template <typename T, typename L, typename U,
143           require_not_std_vector_t<L>* = nullptr>
lub_free(const std::vector<T> y,const L & lb,const std::vector<U> & ub)144 inline auto lub_free(const std::vector<T> y, const L& lb,
145                      const std::vector<U>& ub) {
146   check_matching_dims("lub_free", "y", y, "ub", ub);
147   std::vector<decltype(lub_free(y[0], lb, ub[0]))> ret(y.size());
148   for (Eigen::Index i = 0; i < y.size(); ++i) {
149     ret[i] = lub_free(y[i], lb, ub[i]);
150   }
151   return ret;
152 }
153 
154 /**
155  * Overload for `std::vector` constrained variable and `std::vector` lower bound
156  */
157 template <typename T, typename L, typename U,
158           require_not_std_vector_t<U>* = nullptr>
lub_free(const std::vector<T> y,const std::vector<L> & lb,const U & ub)159 inline auto lub_free(const std::vector<T> y, const std::vector<L>& lb,
160                      const U& ub) {
161   check_matching_dims("lub_free", "y", y, "lb", lb);
162   std::vector<decltype(lub_free(y[0], lb[0], ub))> ret(y.size());
163   for (Eigen::Index i = 0; i < y.size(); ++i) {
164     ret[i] = lub_free(y[i], lb[i], ub);
165   }
166   return ret;
167 }
168 
169 /**
170  * Overload for `std::vector` constrained variable and `std::vector` constraints
171  */
172 template <typename T, typename L, typename U>
lub_free(const std::vector<T> y,const std::vector<L> & lb,const std::vector<U> & ub)173 inline auto lub_free(const std::vector<T> y, const std::vector<L>& lb,
174                      const std::vector<U>& ub) {
175   check_matching_dims("lub_free", "y", y, "lb", lb);
176   check_matching_dims("lub_free", "y", y, "ub", ub);
177   std::vector<decltype(lub_free(y[0], lb[0], ub[0]))> ret(y.size());
178   for (Eigen::Index i = 0; i < y.size(); ++i) {
179     ret[i] = lub_free(y[i], lb[i], ub[i]);
180   }
181   return ret;
182 }
183 ///@}
184 
185 }  // namespace math
186 }  // namespace stan
187 #endif
188