Stan Math Library  2.20.0
reverse mode automatic differentiation
simplex_constrain.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_SIMPLEX_CONSTRAIN_HPP
2 #define STAN_MATH_REV_MAT_FUN_SIMPLEX_CONSTRAIN_HPP
3 
4 #include <stan/math/rev/meta.hpp>
8 #include <tuple>
9 #include <vector>
10 
11 namespace stan {
12 namespace math {
13 
14 namespace internal {
16  int N_;
17  double* diag_; // diagonal of the Jacobian of the operator
18  double* z_;
19 
20  public:
35  template <std::size_t size>
36  Eigen::VectorXd operator()(const std::array<bool, size>& needs_adj,
37  const Eigen::VectorXd& y) {
38  N_ = y.size();
39  diag_ = ChainableStack::instance_->memalloc_.alloc_array<double>(N_);
41 
42  Eigen::Matrix<double, Eigen::Dynamic, 1> x(N_ + 1);
43  double stick_len(1.0);
44  for (int k = 0; k < N_; ++k) {
45  double log_N_minus_k = std::log(N_ - k);
46  z_[k] = inv_logit(y(k) - log_N_minus_k);
47  diag_[k] = stick_len * z_[k] * inv_logit(log_N_minus_k - y(k));
48  x(k) = stick_len * z_[k];
49  stick_len -= x(k);
50  }
51  x(N_) = stick_len;
52  return x;
53  }
54 
55  /*
56  * Compute the result of multiply the transpose of the adjoint vector times
57  * the Jacobian of the simplex_constrain operator.
58  *
59  * @tparam size Number of adjoints to return
60  * @param needs_adj Boolean indicators of if adjoints of arguments will be
61  * needed
62  * @param adj Eigen::VectorXd of adjoints at the output of the softmax
63  * @return Eigen::VectorXd of adjoints propagated through softmax operation
64  */
65  template <std::size_t size>
66  auto multiply_adjoint_jacobian(const std::array<bool, size>& needs_adj,
67  const Eigen::VectorXd& adj) const {
68  Eigen::VectorXd adj_times_jac(N_);
69  double acc = adj(N_);
70 
71  if (N_ > 0) {
72  adj_times_jac(N_ - 1) = diag_[N_ - 1] * (adj(N_ - 1) - acc);
73  for (int n = N_ - 1; --n >= 0;) {
74  acc = adj(n + 1) * z_[n + 1] + (1 - z_[n + 1]) * acc;
75  adj_times_jac(n) = diag_[n] * (adj(n) - acc);
76  }
77  }
78 
79  return std::make_tuple(adj_times_jac);
80  }
81 };
82 } // namespace internal
83 
95 inline Eigen::Matrix<var, Eigen::Dynamic, 1> simplex_constrain(
96  const Eigen::Matrix<var, Eigen::Dynamic, 1>& y) {
97  return adj_jac_apply<internal::simplex_constrain_op>(y);
98 }
99 
100 } // namespace math
101 } // namespace stan
102 #endif
Eigen::VectorXd operator()(const std::array< bool, size > &needs_adj, const Eigen::VectorXd &y)
Return the simplex corresponding to the specified free vector.
fvar< T > log(const fvar< T > &x)
Definition: log.hpp:12
fvar< T > inv_logit(const fvar< T > &x)
Returns the inverse logit function applied to the argument.
Definition: inv_logit.hpp:20
static STAN_THREADS_DEF AutodiffStackStorage * instance_
auto multiply_adjoint_jacobian(const std::array< bool, size > &needs_adj, const Eigen::VectorXd &adj) const
T * alloc_array(size_t n)
Allocate an array on the arena of the specified size to hold values of the specified template paramet...
Eigen::Matrix< T, Eigen::Dynamic, 1 > simplex_constrain(const Eigen::Matrix< T, Eigen::Dynamic, 1 > &y)
Return the simplex corresponding to the specified free vector.

     [ Stan Home Page ] © 2011–2018, Stan Development Team.