Automatic Differentiation
 
Loading...
Searching...
No Matches
lub_constrain.hpp
Go to the documentation of this file.
1#ifndef STAN_MATH_REV_CONSTRAINT_LUB_CONSTRAIN_HPP
2#define STAN_MATH_REV_CONSTRAINT_LUB_CONSTRAIN_HPP
3
4#include <stan/math/rev/core.hpp>
5#include <stan/math/prim/meta.hpp>
6#include <stan/math/rev/constraint/lb_constrain.hpp>
7#include <stan/math/prim/constraint/lub_constrain.hpp>
8#include <stan/math/rev/constraint/ub_constrain.hpp>
9#include <cmath>
10
11namespace stan {
12namespace math {
13
33template <typename T, typename L, typename U,
34 require_all_stan_scalar_t<T, L, U>* = nullptr,
35 require_var_t<return_type_t<T, L, U>>* = nullptr>
36inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
37 using std::exp;
38 const auto lb_val = value_of(lb);
39 const auto ub_val = value_of(ub);
40 const bool is_lb_inf = lb_val == NEGATIVE_INFTY;
41 const bool is_ub_inf = ub_val == INFTY;
42 if (unlikely(is_ub_inf && is_lb_inf)) {
43 return identity_constrain(x, ub, lb);
44 } else if (unlikely(is_ub_inf)) {
45 return lb_constrain(identity_constrain(x, ub), lb);
46 } else if (unlikely(is_lb_inf)) {
47 return ub_constrain(identity_constrain(x, lb), ub);
48 } else {
49 check_less("lub_constrain", "lb", lb_val, ub_val);
50 auto diff = ub_val - lb_val;
51 double inv_logit_x = inv_logit(value_of(x));
52 return make_callback_var(
53 diff * inv_logit_x + lb_val,
54 [x, ub, lb, diff, inv_logit_x](auto& vi) mutable {
55 if (!is_constant<T>::value) {
56 forward_as<var>(x).adj()
57 += vi.adj() * diff * inv_logit_x * (1.0 - inv_logit_x);
58 }
59 if (!is_constant<L>::value) {
60 forward_as<var>(lb).adj() += vi.adj() * (1.0 - inv_logit_x);
61 }
62 if (!is_constant<U>::value) {
63 forward_as<var>(ub).adj() += vi.adj() * inv_logit_x;
64 }
65 });
66 }
67}
68
102template <typename T, typename L, typename U,
103 require_all_stan_scalar_t<T, L, U>* = nullptr,
104 require_var_t<return_type_t<T, L, U>>* = nullptr>
105inline auto lub_constrain(const T& x, const L& lb, const U& ub,
106 return_type_t<T, L, U>& lp) {
107 using std::exp;
108 const auto lb_val = value_of(lb);
109 const auto ub_val = value_of(ub);
110 const bool is_lb_inf = lb_val == NEGATIVE_INFTY;
111 const bool is_ub_inf = ub_val == INFTY;
112 if (unlikely(is_ub_inf && is_lb_inf)) {
113 return identity_constrain(x, ub, lb);
114 } else if (unlikely(is_ub_inf)) {
115 return lb_constrain(identity_constrain(x, ub), lb, lp);
116 } else if (unlikely(is_lb_inf)) {
117 return ub_constrain(identity_constrain(x, lb), ub, lp);
118 } else {
119 check_less("lub_constrain", "lb", lb_val, ub_val);
120 auto neg_abs_x = -abs(value_of(x));
121 auto diff = ub_val - lb_val;
122 double inv_logit_x = inv_logit(value_of(x));
123 lp += (log(diff) + (neg_abs_x - (2.0 * log1p_exp(neg_abs_x))));
124 return make_callback_var(
125 diff * inv_logit_x + lb_val,
126 [x, ub, lb, diff, lp, inv_logit_x](auto& vi) mutable {
127 if (!is_constant<T>::value) {
128 forward_as<var>(x).adj()
129 += vi.adj() * diff * inv_logit_x * (1.0 - inv_logit_x)
130 + lp.adj() * (1.0 - 2.0 * inv_logit_x);
131 }
132 if (!is_constant<L>::value && !is_constant<U>::value) {
133 const auto one_over_diff = 1.0 / diff;
134 forward_as<var>(lb).adj()
135 += vi.adj() * (1.0 - inv_logit_x) + -one_over_diff * lp.adj();
136 forward_as<var>(ub).adj()
137 += vi.adj() * inv_logit_x + one_over_diff * lp.adj();
138 } else if (!is_constant<L>::value) {
139 forward_as<var>(lb).adj()
140 += vi.adj() * (1.0 - inv_logit_x) + (-1.0 / diff) * lp.adj();
141 } else if (!is_constant<U>::value) {
142 forward_as<var>(ub).adj()
143 += vi.adj() * inv_logit_x + (1.0 / diff) * lp.adj();
144 }
145 });
146 }
147}
148
152template <typename T, typename L, typename U, require_matrix_t<T>* = nullptr,
153 require_all_stan_scalar_t<L, U>* = nullptr,
154 require_var_t<return_type_t<T, L, U>>* = nullptr>
155inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
156 using std::exp;
157 using ret_type = return_var_matrix_t<T, T, L, U>;
158 const auto lb_val = value_of(lb);
159 const auto ub_val = value_of(ub);
160 const bool is_lb_inf = lb_val == NEGATIVE_INFTY;
161 const bool is_ub_inf = ub_val == INFTY;
162 if (unlikely(is_ub_inf && is_lb_inf)) {
163 return ret_type(identity_constrain(x, ub, lb));
164 } else if (unlikely(is_ub_inf)) {
165 return ret_type(lb_constrain(identity_constrain(x, ub), lb));
166 } else if (unlikely(is_lb_inf)) {
167 return ret_type(ub_constrain(identity_constrain(x, lb), ub));
168 } else {
169 arena_t<T> arena_x = x;
170 check_less("lub_constrain", "lb", lb_val, ub_val);
171 const auto diff = ub_val - lb_val;
172 auto inv_logit_x = to_arena(inv_logit(arena_x.val().array()));
173 arena_t<ret_type> ret = diff * inv_logit_x + lb_val;
174 reverse_pass_callback([arena_x, ub, lb, ret, diff, inv_logit_x]() mutable {
175 if (!is_constant<T>::value) {
176 using T_var = arena_t<promote_scalar_t<var, T>>;
177 forward_as<T_var>(arena_x).adj().array()
178 += ret.adj().array() * diff * inv_logit_x * (1.0 - inv_logit_x);
179 }
180 if (!is_constant<L>::value) {
181 forward_as<var>(lb).adj()
182 += (ret.adj().array() * (1.0 - inv_logit_x)).sum();
183 }
184 if (!is_constant<U>::value) {
185 forward_as<var>(ub).adj() += (ret.adj().array() * inv_logit_x).sum();
186 }
187 });
188 return ret_type(ret);
189 }
190}
191
195template <typename T, typename L, typename U, require_matrix_t<T>* = nullptr,
196 require_all_stan_scalar_t<L, U>* = nullptr,
197 require_var_t<return_type_t<T, L, U>>* = nullptr>
198inline auto lub_constrain(const T& x, const L& lb, const U& ub,
199 return_type_t<T, L, U>& lp) {
200 using std::exp;
201 using ret_type = return_var_matrix_t<T, T, L, U>;
202 const auto lb_val = value_of(lb);
203 const auto ub_val = value_of(ub);
204 const bool is_lb_inf = lb_val == NEGATIVE_INFTY;
205 const bool is_ub_inf = ub_val == INFTY;
206 if (unlikely(is_ub_inf && is_lb_inf)) {
207 return ret_type(identity_constrain(x, ub, lb));
208 } else if (unlikely(is_ub_inf)) {
209 return ret_type(lb_constrain(identity_constrain(x, ub), lb, lp));
210 } else if (unlikely(is_lb_inf)) {
211 return ret_type(ub_constrain(identity_constrain(x, lb), ub, lp));
212 } else {
213 check_less("lub_constrain", "lb", lb_val, ub_val);
214 arena_t<T> arena_x = x;
215 auto neg_abs_x = to_arena(-(value_of(arena_x).array()).abs());
216 auto diff = ub_val - lb_val;
217 lp += (log(diff) + (neg_abs_x - (2.0 * log1p_exp(neg_abs_x)))).sum();
218 auto inv_logit_x = to_arena(inv_logit(value_of(arena_x).array()));
219 arena_t<ret_type> ret = diff * inv_logit_x + lb_val;
220 reverse_pass_callback(
221 [arena_x, ub, lb, ret, lp, diff, inv_logit_x]() mutable {
222 if (!is_constant<T>::value) {
223 forward_as<arena_t<promote_scalar_t<var, T>>>(arena_x).adj().array()
224 += ret.adj().array() * diff * inv_logit_x * (1.0 - inv_logit_x)
225 + lp.adj() * (1.0 - 2.0 * inv_logit_x);
226 }
227 if (!is_constant<L>::value && !is_constant<U>::value) {
228 const auto lp_calc = lp.adj() * ret.size();
229 const auto one_over_diff = 1.0 / diff;
230 forward_as<var>(lb).adj()
231 += (ret.adj().array() * (1.0 - inv_logit_x)).sum()
232 + -one_over_diff * lp_calc;
233 forward_as<var>(ub).adj() += (ret.adj().array() * inv_logit_x).sum()
234 + one_over_diff * lp_calc;
235 } else if (!is_constant<L>::value) {
236 forward_as<var>(lb).adj()
237 += (ret.adj().array() * (1.0 - inv_logit_x)).sum()
238 + -(1.0 / diff) * lp.adj() * ret.size();
239 } else if (!is_constant<U>::value) {
240 forward_as<var>(ub).adj() += (ret.adj().array() * inv_logit_x).sum()
241 + (1.0 / diff) * lp.adj() * ret.size();
242 }
243 });
244 return ret_type(ret);
245 }
246}
247
252template <typename T, typename L, typename U,
253 require_all_matrix_t<T, L>* = nullptr,
254 require_stan_scalar_t<U>* = nullptr,
255 require_var_t<return_type_t<T, L, U>>* = nullptr>
256inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
257 using std::exp;
258 using ret_type = return_var_matrix_t<T, T, L, U>;
259 const auto ub_val = value_of(ub);
260 const bool is_ub_inf = ub_val == INFTY;
261 if (unlikely(is_ub_inf)) {
262 return eval(lb_constrain(identity_constrain(x, ub), lb));
263 } else {
264 arena_t<T> arena_x = x;
265 arena_t<L> arena_lb = lb;
266 const auto lb_val = value_of(arena_lb).array().eval();
267 check_less("lub_constrain", "lb", lb_val, ub_val);
268 auto is_lb_inf = to_arena((lb_val == NEGATIVE_INFTY));
269 auto diff = to_arena(ub_val - lb_val);
270 auto inv_logit_x = to_arena(inv_logit(value_of(arena_x).array()));
271 arena_t<ret_type> ret = (is_lb_inf).select(
272 ub_val - value_of(arena_x).array().exp(), diff * inv_logit_x + lb_val);
273 reverse_pass_callback([arena_x, ub, arena_lb, ret, diff, inv_logit_x,
274 is_lb_inf]() mutable {
275 using T_var = arena_t<promote_scalar_t<var, T>>;
276 using L_var = arena_t<promote_scalar_t<var, L>>;
277 if (!is_constant<T>::value) {
278 forward_as<T_var>(arena_x).adj().array() += (is_lb_inf).select(
279 ret.adj().array() * -value_of(arena_x).array().exp(),
280 ret.adj().array() * diff * inv_logit_x * (1.0 - inv_logit_x));
281 }
282 if (!is_constant<U>::value) {
283 forward_as<var>(ub).adj()
284 += (is_lb_inf)
285 .select(ret.adj().array(), ret.adj().array() * inv_logit_x)
286 .sum();
287 }
288 if (!is_constant<L>::value) {
289 forward_as<L_var>(arena_lb).adj().array()
290 += (is_lb_inf).select(0, ret.adj().array() * (1.0 - inv_logit_x));
291 }
292 });
293 return ret_type(ret);
294 }
295}
296
301template <typename T, typename L, typename U,
302 require_all_matrix_t<T, L>* = nullptr,
303 require_stan_scalar_t<U>* = nullptr,
304 require_var_t<return_type_t<T, L, U>>* = nullptr>
305inline auto lub_constrain(const T& x, const L& lb, const U& ub,
306 std::decay_t<return_type_t<T, L, U>>& lp) {
307 using std::exp;
308 using ret_type = return_var_matrix_t<T, T, L, U>;
309 const auto ub_val = value_of(ub);
310 const bool is_ub_inf = ub_val == INFTY;
311 if (unlikely(is_ub_inf)) {
312 return ret_type(lb_constrain(identity_constrain(x, ub), lb, lp));
313 } else {
314 arena_t<T> arena_x = x;
315 arena_t<L> arena_lb = lb;
316 const auto arena_x_val = to_arena(value_of(arena_x).array());
317 const auto lb_val = value_of(arena_lb).array().eval();
318 check_less("lub_constrain", "lb", lb_val, ub_val);
319 auto is_lb_inf = to_arena((lb_val == NEGATIVE_INFTY));
320 auto diff = to_arena(ub_val - lb_val);
321 auto neg_abs_x = to_arena(-arena_x_val.abs());
322 auto inv_logit_x = to_arena(inv_logit(arena_x_val));
323 arena_t<ret_type> ret = (is_lb_inf).select(ub_val - arena_x_val.exp(),
324 diff * inv_logit_x + lb_val);
325 lp += (is_lb_inf)
326 .select(arena_x_val,
327 log(diff) + (neg_abs_x - (2.0 * log1p_exp(neg_abs_x))))
328 .sum();
329 reverse_pass_callback([arena_x, arena_x_val, ub, arena_lb, ret, lp, diff,
330 inv_logit_x, is_lb_inf]() mutable {
331 using T_var = arena_t<promote_scalar_t<var, T>>;
332 using L_var = arena_t<promote_scalar_t<var, L>>;
333 const auto lp_adj = lp.adj();
334 if (!is_constant<T>::value) {
335 const auto x_sign = arena_x_val.sign().eval();
336 forward_as<T_var>(arena_x).adj().array() += (is_lb_inf).select(
337 ret.adj().array() * -arena_x_val.exp() + lp_adj,
338 ret.adj().array() * diff * inv_logit_x * (1.0 - inv_logit_x)
339 + lp.adj() * (1.0 - 2.0 * inv_logit_x));
340 }
341 if (!is_constant<L>::value) {
342 forward_as<L_var>(arena_lb).adj().array()
343 += (is_lb_inf).select(0, ret.adj().array() * (1.0 - inv_logit_x)
344 + -(1.0 / diff) * lp_adj);
345 }
346 if (!is_constant<U>::value) {
347 forward_as<var>(ub).adj()
348 += (is_lb_inf)
349 .select(ret.adj().array(), ret.adj().array() * inv_logit_x
350 + (1.0 / diff) * lp_adj)
351 .sum();
352 }
353 });
354 return ret_type(ret);
355 }
356}
357
362template <typename T, typename L, typename U,
363 require_all_matrix_t<T, U>* = nullptr,
364 require_stan_scalar_t<L>* = nullptr,
365 require_var_t<return_type_t<T, L, U>>* = nullptr>
366inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
367 using std::exp;
368 using ret_type = return_var_matrix_t<T, T, L, U>;
369 const auto lb_val = value_of(lb);
370 const bool is_lb_inf = lb_val == NEGATIVE_INFTY;
371 if (unlikely(is_lb_inf)) {
372 return eval(ub_constrain(identity_constrain(x, lb), ub));
373 } else {
374 arena_t<T> arena_x = x;
375 auto arena_x_val = to_arena(value_of(arena_x));
376 arena_t<U> arena_ub = ub;
377 const auto ub_val = value_of(arena_ub).array().eval();
378 check_less("lub_constrain", "lb", lb_val, ub_val);
379 auto is_ub_inf = to_arena((ub_val == INFTY));
380 auto diff = to_arena(ub_val - lb_val);
381 auto inv_logit_x = to_arena(inv_logit(arena_x_val.array()));
382 arena_t<ret_type> ret = (is_ub_inf).select(
383 arena_x_val.array().exp() + lb_val, diff * inv_logit_x + lb_val);
384 reverse_pass_callback([arena_x, arena_x_val, arena_ub, lb, ret, is_ub_inf,
385 inv_logit_x, diff]() mutable {
386 using T_var = arena_t<promote_scalar_t<var, T>>;
387 using U_var = arena_t<promote_scalar_t<var, U>>;
388 if (!is_constant<T>::value) {
389 forward_as<T_var>(arena_x).adj().array() += (is_ub_inf).select(
390 ret.adj().array() * arena_x_val.array().exp(),
391 ret.adj().array() * diff * inv_logit_x * (1.0 - inv_logit_x));
392 }
393 if (!is_constant<L>::value) {
394 forward_as<var>(lb).adj()
395 += (is_ub_inf)
396 .select(ret.adj().array(),
397 ret.adj().array() * (1.0 - inv_logit_x))
398 .sum();
399 }
400 if (!is_constant<U>::value) {
401 forward_as<U_var>(arena_ub).adj().array()
402 += (is_ub_inf).select(0, ret.adj().array() * inv_logit_x);
403 }
404 });
405 return ret_type(ret);
406 }
407}
408
413template <typename T, typename L, typename U,
414 require_all_matrix_t<T, U>* = nullptr,
415 require_stan_scalar_t<L>* = nullptr,
416 require_var_t<return_type_t<T, L, U>>* = nullptr>
417inline auto lub_constrain(const T& x, const L& lb, const U& ub,
418 std::decay_t<return_type_t<T, L, U>>& lp) {
419 using std::exp;
420 using ret_type = return_var_matrix_t<T, T, L, U>;
421 const auto lb_val = value_of(lb);
422 const bool is_lb_inf = lb_val == NEGATIVE_INFTY;
423 if (unlikely(is_lb_inf)) {
424 return eval(ub_constrain(identity_constrain(x, lb), ub, lp));
425 } else {
426 arena_t<T> arena_x = x;
427 auto arena_x_val = to_arena(value_of(arena_x));
428 arena_t<U> arena_ub = ub;
429 const auto& ub_val = to_ref(value_of(arena_ub));
430 check_less("lub_constrain", "lb", lb_val, ub_val);
431 auto is_ub_inf = to_arena((ub_val.array() == INFTY));
432 auto diff = to_arena(ub_val.array() - lb_val);
433 auto neg_abs_x = to_arena(-(arena_x_val.array()).abs());
434 lp += (is_ub_inf)
435 .select(arena_x_val.array(),
436 log(diff) + (neg_abs_x - (2.0 * log1p_exp(neg_abs_x))))
437 .sum();
438 auto inv_logit_x = to_arena(inv_logit(arena_x_val.array()));
439 arena_t<ret_type> ret = (is_ub_inf).select(
440 arena_x_val.array().exp() + lb_val, diff * inv_logit_x + lb_val);
441 reverse_pass_callback([arena_x, arena_x_val, diff, inv_logit_x, arena_ub,
442 lb, ret, lp, is_ub_inf]() mutable {
443 using T_var = arena_t<promote_scalar_t<var, T>>;
444 using U_var = arena_t<promote_scalar_t<var, U>>;
445 const auto lp_adj = lp.adj();
446 if (!is_constant<T>::value) {
447 forward_as<T_var>(arena_x).adj().array() += (is_ub_inf).select(
448 ret.adj().array() * arena_x_val.array().exp() + lp_adj,
449 ret.adj().array() * diff * inv_logit_x * (1.0 - inv_logit_x)
450 + lp.adj() * (1.0 - 2.0 * inv_logit_x));
451 }
452 if (!is_constant<L>::value) {
453 forward_as<var>(lb).adj()
454 += (is_ub_inf)
455 .select(ret.adj().array(),
456 ret.adj().array() * (1.0 - inv_logit_x)
457 + -(1.0 / diff) * lp_adj)
458 .sum();
459 }
460 if (!is_constant<U>::value) {
461 forward_as<U_var>(arena_ub).adj().array() += (is_ub_inf).select(
462 0, ret.adj().array() * inv_logit_x + (1.0 / diff) * lp_adj);
463 }
464 });
465 return ret_type(ret);
466 }
467}
468
472template <typename T, typename L, typename U,
473 require_all_matrix_t<T, L, U>* = nullptr,
474 require_var_t<return_type_t<T, L, U>>* = nullptr>
475inline auto lub_constrain(const T& x, const L& lb, const U& ub) {
476 using std::exp;
477 using ret_type = return_var_matrix_t<T, T, L, U>;
478 arena_t<T> arena_x = x;
479 auto arena_x_val = value_of(arena_x);
480 arena_t<L> arena_lb = lb;
481 arena_t<U> arena_ub = ub;
482 auto lb_val = value_of(arena_lb).array();
483 auto ub_val = value_of(arena_ub).array();
484 check_less("lub_constrain", "lb", lb_val, ub_val);
485 auto inv_logit_x = to_arena(inv_logit(arena_x_val.array()));
486 auto is_lb_inf = to_arena((lb_val == NEGATIVE_INFTY));
487 auto is_ub_inf = to_arena((ub_val == INFTY));
488 auto is_lb_ub_inf = to_arena(is_lb_inf && is_ub_inf);
489 auto diff = to_arena(ub_val - lb_val);
490 // if both, identity, then lb_inf -> ub_constrain, then ub_inf -> lb_constrain
491 arena_t<ret_type> ret
492 = (is_lb_ub_inf)
493 .select(arena_x_val.array(),
494 (is_lb_inf).select(
495 ub_val - arena_x.val().array().exp(),
496 (is_ub_inf).select(arena_x_val.array().exp() + lb_val,
497 diff * inv_logit_x + lb_val)));
498 reverse_pass_callback([arena_x, arena_x_val, inv_logit_x, arena_ub, arena_lb,
499 diff, ret, is_ub_inf, is_lb_inf,
500 is_lb_ub_inf]() mutable {
501 using T_var = arena_t<promote_scalar_t<var, T>>;
502 using L_var = arena_t<promote_scalar_t<var, L>>;
503 using U_var = arena_t<promote_scalar_t<var, U>>;
504 // The most likely case is neither of them are infinity
505 const bool is_none_inf = !(is_lb_inf.any() || is_ub_inf.any());
506 if (is_none_inf) {
507 if (!is_constant<T>::value) {
508 forward_as<T_var>(arena_x).adj().array()
509 += ret.adj().array() * diff * inv_logit_x * (1.0 - inv_logit_x);
510 }
511 if (!is_constant<L>::value) {
512 forward_as<L_var>(arena_lb).adj().array()
513 += ret.adj().array() * (1.0 - inv_logit_x);
514 }
515 if (!is_constant<U>::value) {
516 forward_as<U_var>(arena_ub).adj().array()
517 += ret.adj().array() * inv_logit_x;
518 }
519 } else {
520 if (!is_constant<T>::value) {
521 forward_as<T_var>(arena_x).adj().array()
522 += (is_lb_ub_inf)
523 .select(
524 ret.adj().array(),
525 (is_lb_inf).select(
526 ret.adj().array() * -arena_x_val.array().exp(),
527 (is_ub_inf).select(
528 ret.adj().array() * arena_x_val.array().exp(),
529 ret.adj().array() * diff * inv_logit_x
530 * (1.0 - inv_logit_x))));
531 }
532 if (!is_constant<U>::value) {
533 forward_as<U_var>(arena_ub).adj().array() += (is_ub_inf).select(
534 0, (is_lb_inf).select(ret.adj().array(),
535 ret.adj().array() * inv_logit_x));
536 }
537 if (!is_constant<L>::value) {
538 forward_as<L_var>(arena_lb).adj().array() += (is_lb_inf).select(
539 0, (is_ub_inf).select(ret.adj().array(),
540 ret.adj().array() * (1.0 - inv_logit_x)));
541 }
542 }
543 });
544 return ret_type(ret);
545}
546
550template <typename T, typename L, typename U,
551 require_all_matrix_t<T, L, U>* = nullptr,
552 require_var_t<return_type_t<T, L, U>>* = nullptr>
553inline auto lub_constrain(const T& x, const L& lb, const U& ub,
554 return_type_t<T, L, U>& lp) {
555 using std::exp;
556 using ret_type = return_var_matrix_t<T, T, L, U>;
557 arena_t<T> arena_x = x;
558 auto arena_x_val = value_of(arena_x);
559 arena_t<L> arena_lb = lb;
560 arena_t<U> arena_ub = ub;
561 auto lb_val = value_of(arena_lb).array();
562 auto ub_val = value_of(arena_ub).array();
563 check_less("lub_constrain", "lb", lb_val, ub_val);
564 auto inv_logit_x = to_arena(inv_logit(arena_x_val.array()));
565 auto is_lb_inf = to_arena((lb_val == NEGATIVE_INFTY));
566 auto is_ub_inf = to_arena((ub_val == INFTY));
567 auto is_lb_ub_inf = to_arena(is_lb_inf && is_ub_inf);
568 auto diff = to_arena(ub_val - lb_val);
569 // if both, identity, then lb_inf -> ub_constrain, then ub_inf -> lb_constrain
570 arena_t<ret_type> ret
571 = (is_lb_ub_inf)
572 .select(arena_x_val.array(),
573 (is_lb_inf).select(
574 ub_val - arena_x.val().array().exp(),
575 (is_ub_inf).select(arena_x_val.array().exp() + lb_val,
576 diff * inv_logit_x + lb_val)));
577 auto neg_abs_x = to_arena(-(arena_x_val.array()).abs());
578 lp += (is_lb_ub_inf)
579 .select(
580 0,
581 (is_lb_inf || is_ub_inf)
582 .select(
583 arena_x_val.array(),
584 log(diff) + (neg_abs_x - (2.0 * log1p_exp(neg_abs_x)))))
585 .sum();
586 reverse_pass_callback([arena_x, arena_x_val, inv_logit_x, arena_ub, arena_lb,
587 diff, ret, is_ub_inf, is_lb_inf, is_lb_ub_inf,
588 lp]() mutable {
589 using T_var = arena_t<promote_scalar_t<var, T>>;
590 using L_var = arena_t<promote_scalar_t<var, L>>;
591 using U_var = arena_t<promote_scalar_t<var, U>>;
592 const auto lp_adj = lp.adj();
593 // The most likely case is neither of them are infinity
594 const bool is_none_inf = !(is_lb_inf.any() || is_ub_inf.any());
595 if (is_none_inf) {
596 if (!is_constant<T>::value) {
597 forward_as<T_var>(arena_x).adj().array()
598 += ret.adj().array() * diff * inv_logit_x * (1.0 - inv_logit_x)
599 + lp.adj() * (1.0 - 2.0 * inv_logit_x);
600 }
601 if (!is_constant<L>::value) {
602 forward_as<L_var>(arena_lb).adj().array()
603 += ret.adj().array() * (1.0 - inv_logit_x) + -(1.0 / diff) * lp_adj;
604 }
605 if (!is_constant<U>::value) {
606 forward_as<U_var>(arena_ub).adj().array()
607 += ret.adj().array() * inv_logit_x + (1.0 / diff) * lp_adj;
608 }
609 } else {
610 if (!is_constant<T>::value) {
611 forward_as<T_var>(arena_x).adj().array()
612 += (is_lb_ub_inf)
613 .select(
614 ret.adj().array(),
615 (is_lb_inf).select(
616 ret.adj().array() * -arena_x_val.array().exp()
617 + lp_adj,
618 (is_ub_inf).select(
619 ret.adj().array() * arena_x_val.array().exp()
620 + lp_adj,
621 ret.adj().array() * diff * inv_logit_x
622 * (1.0 - inv_logit_x)
623 + lp.adj() * (1.0 - 2.0 * inv_logit_x))));
624 }
625 if (!is_constant<L>::value) {
626 forward_as<L_var>(arena_lb).adj().array() += (is_lb_inf).select(
627 0, (is_ub_inf).select(ret.adj().array(),
628 ret.adj().array() * (1.0 - inv_logit_x)
629 + -(1.0 / diff) * lp_adj));
630 }
631 if (!is_constant<U>::value) {
632 forward_as<U_var>(arena_ub).adj().array() += (is_ub_inf).select(
633 0, (is_lb_inf).select(
634 ret.adj().array(),
635 ret.adj().array() * inv_logit_x + (1.0 / diff) * lp_adj));
636 }
637 }
638 });
639 return ret_type(ret);
640}
641
642} // namespace math
643} // namespace stan
644
645#endif
unlikely
#define unlikely(x)
Definition compiler_attributes.hpp:31
stan::require_all_matrix_t
require_all_t< is_matrix< std::decay_t< Types > >... > require_all_matrix_t
Require all of the types satisfy is_matrix.
Definition is_matrix.hpp:38
stan::math::select
select_< as_operation_cl_t< T_condition >, as_operation_cl_t< T_then >, as_operation_cl_t< T_else > > select(T_condition &&condition, T_then &&then, T_else &&els)
Selection operation on kernel generator expressions.
Definition select.hpp:148
stan::require_stan_scalar_t
require_t< is_stan_scalar< std::decay_t< T > > > require_stan_scalar_t
Require type satisfies is_stan_scalar.
Definition is_stan_scalar.hpp:39
stan::return_type_t
typename return_type< Ts... >::type return_type_t
Convenience type for the return type of the specified template parameters.
Definition return_type.hpp:218
stan::require_var_t
require_t< is_var< std::decay_t< T > > > require_var_t
Require type satisfies is_var.
Definition is_var.hpp:24
stan::math::abs
fvar< T > abs(const fvar< T > &x)
Definition abs.hpp:15
stan::math::make_callback_var
var_value< plain_type_t< T > > make_callback_var(T &&value, F &&functor)
Creates a new var initialized with a callback_vari with a given value and reverse-pass callback funct...
Definition callback_vari.hpp:61
stan::math::reverse_pass_callback
void reverse_pass_callback(F &&functor)
Puts a callback on the autodiff stack to be called in reverse pass.
Definition reverse_pass_callback.hpp:38
stan::math::eval
T eval(T &&arg)
Inputs which have a plain_type equal to the own time are forwarded unmodified (for Eigen expressions ...
Definition eval.hpp:20
stan::math::value_of
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition value_of.hpp:18
stan::math::log
fvar< T > log(const fvar< T > &x)
Definition log.hpp:18
stan::math::NEGATIVE_INFTY
static constexpr double NEGATIVE_INFTY
Negative infinity.
Definition constants.hpp:51
stan::math::log1p_exp
fvar< T > log1p_exp(const fvar< T > &x)
Definition log1p_exp.hpp:14
stan::math::to_arena
arena_t< T > to_arena(const T &a)
Converts given argument into a type that either has any dynamic allocation on AD stack or schedules i...
Definition to_arena.hpp:25
stan::math::to_ref
ref_type_t< T && > to_ref(T &&a)
This evaluates expensive Eigen expressions.
Definition to_ref.hpp:17
stan::math::lb_constrain
auto lb_constrain(T &&x, L &&lb)
Return the lower-bounded value for the specified unconstrained input and specified lower bound.
Definition lb_constrain.hpp:32
stan::math::sum
auto sum(const std::vector< T > &m)
Return the sum of the entries of the specified standard vector.
Definition sum.hpp:23
stan::math::ub_constrain
auto ub_constrain(T &&x, U &&ub)
Return the upper-bounded value for the specified unconstrained matrix and upper bound.
Definition ub_constrain.hpp:32
stan::math::lub_constrain
matrix_cl< double > lub_constrain(const T &x, const L &lb, const U &ub)
Return the lower and upper-bounded matrix derived by transforming the specified free matrix given the...
Definition lub_constrain.hpp:35
stan::math::check_less
void check_less(const char *function, const char *name, const T_y &y, const T_high &high, Idxs... idxs)
Throw an exception if y is not strictly less than high.
Definition check_less.hpp:35
stan::math::inv_logit
fvar< T > inv_logit(const fvar< T > &x)
Returns the inverse logit function applied to the argument.
Definition inv_logit.hpp:20
stan::math::identity_constrain
auto identity_constrain(T &&x, Types &&...)
Returns the result of applying the identity constraint transform to the input.
Definition identity_constrain.hpp:23
stan::math::INFTY
static constexpr double INFTY
Positive infinity.
Definition constants.hpp:46
stan::arena_t
typename internal::arena_type_impl< std::decay_t< T > >::type arena_t
Determines a type that can be used in place of T that does any dynamic allocations on the AD stack.
Definition arena_type.hpp:58
stan::return_var_matrix_t
std::conditional_t< is_any_var_matrix< ReturnType, Types... >::value, stan::math::var_value< stan::math::promote_scalar_t< double, plain_type_t< ReturnType > > >, stan::math::promote_scalar_t< stan::math::var_value< double >, plain_type_t< ReturnType > > > return_var_matrix_t
Given an Eigen type and several inputs, determine if a matrix should be var<Matrix> or Matrix<var>.
Definition return_var_matrix.hpp:27
stan
The lgamma implementation in stan-math is based on either the reentrant safe lgamma_r implementation ...
Definition unit_vector_constrain.hpp:15
lub_constrain.hpp
lb_constrain.hpp
ub_constrain.hpp
stan::is_constant
Metaprogramming struct to detect whether a given type is constant in the mathematical sense (not the ...
Definition is_constant.hpp:30