Stan Math Library  2.20.0
reverse mode automatic differentiation
dot_product.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_REV_MAT_FUN_DOT_PRODUCT_HPP
2 #define STAN_MATH_REV_MAT_FUN_DOT_PRODUCT_HPP
3 
4 #include <stan/math/rev/meta.hpp>
9 #include <stan/math/rev/core.hpp>
11 #include <type_traits>
12 #include <vector>
13 
14 namespace stan {
15 namespace math {
16 
17 namespace internal {
18 template <typename T>
20 
21 template <>
23  typedef vari** type;
24 };
25 
26 template <>
27 struct dot_product_store_type<double> {
28  typedef double* type;
29 };
30 
31 template <typename T1, typename T2>
32 class dot_product_vari : public vari {
33  protected:
36  size_t length_;
37 
38  inline static double var_dot(vari** v1, vari** v2, size_t length) {
39  Eigen::VectorXd vd1(length), vd2(length);
40  for (size_t i = 0; i < length; i++) {
41  vd1[i] = v1[i]->val_;
42  vd2[i] = v2[i]->val_;
43  }
44  return vd1.dot(vd2);
45  }
46 
47  inline static double var_dot(const T1* v1, const T2* v2, size_t length) {
48  Eigen::VectorXd vd1(length), vd2(length);
49  for (size_t i = 0; i < length; i++) {
50  vd1[i] = value_of(v1[i]);
51  vd2[i] = value_of(v2[i]);
52  }
53  return vd1.dot(vd2);
54  }
55 
56  template <typename Derived1, typename Derived2>
57  inline static double var_dot(const Eigen::DenseBase<Derived1>& v1,
58  const Eigen::DenseBase<Derived2>& v2) {
59  Eigen::VectorXd vd1(v1.size()), vd2(v1.size());
60  for (int i = 0; i < v1.size(); i++) {
61  vd1[i] = value_of(v1[i]);
62  vd2[i] = value_of(v2[i]);
63  }
64  return vd1.dot(vd2);
65  }
66  inline void chain(vari** v1, vari** v2) {
67  for (size_t i = 0; i < length_; i++) {
68  v1[i]->adj_ += adj_ * v2_[i]->val_;
69  v2[i]->adj_ += adj_ * v1_[i]->val_;
70  }
71  }
72  inline void chain(double* v1, vari** v2) {
73  for (size_t i = 0; i < length_; i++) {
74  v2[i]->adj_ += adj_ * v1_[i];
75  }
76  }
77  inline void chain(vari** v1, double* v2) {
78  for (size_t i = 0; i < length_; i++) {
79  v1[i]->adj_ += adj_ * v2_[i];
80  }
81  }
82  inline void initialize(vari**& mem_v, const var* inv,
83  vari** shared = nullptr) {
84  if (shared == nullptr) {
85  mem_v = reinterpret_cast<vari**>(
86  ChainableStack::instance_->memalloc_.alloc(length_ * sizeof(vari*)));
87  for (size_t i = 0; i < length_; i++)
88  mem_v[i] = inv[i].vi_;
89  } else {
90  mem_v = shared;
91  }
92  }
93  template <typename Derived>
94  inline void initialize(vari**& mem_v, const Eigen::DenseBase<Derived>& inv,
95  vari** shared = nullptr) {
96  if (shared == nullptr) {
97  mem_v = reinterpret_cast<vari**>(
98  ChainableStack::instance_->memalloc_.alloc(length_ * sizeof(vari*)));
99  for (size_t i = 0; i < length_; i++)
100  mem_v[i] = inv(i).vi_;
101  } else {
102  mem_v = shared;
103  }
104  }
105 
106  inline void initialize(double*& mem_d, const double* ind,
107  double* shared = nullptr) {
108  if (shared == nullptr) {
109  mem_d = reinterpret_cast<double*>(
110  ChainableStack::instance_->memalloc_.alloc(length_ * sizeof(double)));
111  for (size_t i = 0; i < length_; i++)
112  mem_d[i] = ind[i];
113  } else {
114  mem_d = shared;
115  }
116  }
117  template <typename Derived>
118  inline void initialize(double*& mem_d, const Eigen::DenseBase<Derived>& ind,
119  double* shared = nullptr) {
120  if (shared == nullptr) {
121  mem_d = reinterpret_cast<double*>(
122  ChainableStack::instance_->memalloc_.alloc(length_ * sizeof(double)));
123  for (size_t i = 0; i < length_; i++)
124  mem_d[i] = ind(i);
125  } else {
126  mem_d = shared;
127  }
128  }
129 
130  public:
132  typename dot_product_store_type<T2>::type v2, size_t length)
133  : vari(var_dot(v1, v2, length)), v1_(v1), v2_(v2), length_(length) {}
134 
135  dot_product_vari(const T1* v1, const T2* v2, size_t length,
136  dot_product_vari<T1, T2>* shared_v1 = NULL,
137  dot_product_vari<T1, T2>* shared_v2 = NULL)
138  : vari(var_dot(v1, v2, length)), length_(length) {
139  if (shared_v1 == NULL) {
140  initialize(v1_, v1);
141  } else {
142  initialize(v1_, v1, shared_v1->v1_);
143  }
144  if (shared_v2 == NULL) {
145  initialize(v2_, v2);
146  } else {
147  initialize(v2_, v2, shared_v2->v2_);
148  }
149  }
150  template <typename Derived1, typename Derived2>
151  dot_product_vari(const Eigen::DenseBase<Derived1>& v1,
152  const Eigen::DenseBase<Derived2>& v2,
153  dot_product_vari<T1, T2>* shared_v1 = NULL,
154  dot_product_vari<T1, T2>* shared_v2 = NULL)
155  : vari(var_dot(v1, v2)), length_(v1.size()) {
156  if (shared_v1 == NULL) {
157  initialize(v1_, v1);
158  } else {
159  initialize(v1_, v1, shared_v1->v1_);
160  }
161  if (shared_v2 == NULL) {
162  initialize(v2_, v2);
163  } else {
164  initialize(v2_, v2, shared_v2->v2_);
165  }
166  }
167  template <int R1, int C1, int R2, int C2>
168  dot_product_vari(const Eigen::Matrix<T1, R1, C1>& v1,
169  const Eigen::Matrix<T2, R2, C2>& v2,
170  dot_product_vari<T1, T2>* shared_v1 = NULL,
171  dot_product_vari<T1, T2>* shared_v2 = NULL)
172  : vari(var_dot(v1, v2)), length_(v1.size()) {
173  if (shared_v1 == NULL) {
174  initialize(v1_, v1);
175  } else {
176  initialize(v1_, v1, shared_v1->v1_);
177  }
178  if (shared_v2 == NULL) {
179  initialize(v2_, v2);
180  } else {
181  initialize(v2_, v2, shared_v2->v2_);
182  }
183  }
184  virtual void chain() { chain(v1_, v2_); }
185 };
186 } // namespace internal
187 
196 template <typename T1, int R1, int C1, typename T2, int R2, int C2>
197 inline typename std::enable_if<
198  std::is_same<T1, var>::value || std::is_same<T2, var>::value, var>::type
199 dot_product(const Eigen::Matrix<T1, R1, C1>& v1,
200  const Eigen::Matrix<T2, R2, C2>& v2) {
201  check_vector("dot_product", "v1", v1);
202  check_vector("dot_product", "v2", v2);
203  check_matching_sizes("dot_product", "v1", v1, "v2", v2);
204  return var(new internal::dot_product_vari<T1, T2>(v1, v2));
205 }
214 template <typename T1, typename T2>
215 inline typename std::enable_if<
216  std::is_same<T1, var>::value || std::is_same<T2, var>::value, var>::type
217 dot_product(const T1* v1, const T2* v2, size_t length) {
218  return var(new internal::dot_product_vari<T1, T2>(v1, v2, length));
219 }
220 
229 template <typename T1, typename T2>
230 inline typename std::enable_if<
231  std::is_same<T1, var>::value || std::is_same<T2, var>::value, var>::type
232 dot_product(const std::vector<T1>& v1, const std::vector<T2>& v2) {
233  check_matching_sizes("dot_product", "v1", v1, "v2", v2);
234  return var(new internal::dot_product_vari<T1, T2>(&v1[0], &v2[0], v1.size()));
235 }
236 
237 } // namespace math
238 } // namespace stan
239 #endif
void chain(double *v1, vari **v2)
Definition: dot_product.hpp:72
dot_product_vari(const Eigen::Matrix< T1, R1, C1 > &v1, const Eigen::Matrix< T2, R2, C2 > &v2, dot_product_vari< T1, T2 > *shared_v1=NULL, dot_product_vari< T1, T2 > *shared_v2=NULL)
void initialize(vari **&mem_v, const Eigen::DenseBase< Derived > &inv, vari **shared=nullptr)
Definition: dot_product.hpp:94
void initialize(double *&mem_d, const Eigen::DenseBase< Derived > &ind, double *shared=nullptr)
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:17
void initialize(vari **&mem_v, const var *inv, vari **shared=nullptr)
Definition: dot_product.hpp:82
void chain(vari **v1, vari **v2)
Definition: dot_product.hpp:66
The variable implementation base class.
Definition: vari.hpp:30
void check_vector(const char *function, const char *name, const Eigen::Matrix< T, R, C > &x)
Check if the matrix is either a row vector or column vector.
size_t length(const std::vector< T > &x)
Returns the length of the provided std::vector.
Definition: length.hpp:16
static double var_dot(const T1 *v1, const T2 *v2, size_t length)
Definition: dot_product.hpp:47
static double var_dot(vari **v1, vari **v2, size_t length)
Definition: dot_product.hpp:38
static STAN_THREADS_DEF AutodiffStackStorage * instance_
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:33
dot_product_vari(typename dot_product_store_type< T1 >::type v1, typename dot_product_store_type< T2 >::type v2, size_t length)
dot_product_store_type< T1 >::type v1_
Definition: dot_product.hpp:34
const double val_
The value of this variable.
Definition: vari.hpp:38
void initialize(T &x, const T &v)
Definition: initialize.hpp:15
virtual void chain()
Apply the chain rule to this variable based on the variables on which it depends. ...
dot_product_vari(const Eigen::DenseBase< Derived1 > &v1, const Eigen::DenseBase< Derived2 > &v2, dot_product_vari< T1, T2 > *shared_v1=NULL, dot_product_vari< T1, T2 > *shared_v2=NULL)
fvar< T > dot_product(const Eigen::Matrix< fvar< T >, R1, C1 > &v1, const Eigen::Matrix< fvar< T >, R2, C2 > &v2)
Definition: dot_product.hpp:14
void chain(vari **v1, double *v2)
Definition: dot_product.hpp:77
dot_product_vari(const T1 *v1, const T2 *v2, size_t length, dot_product_vari< T1, T2 > *shared_v1=NULL, dot_product_vari< T1, T2 > *shared_v2=NULL)
dot_product_store_type< T2 >::type v2_
Definition: dot_product.hpp:35
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
void check_matching_sizes(const char *function, const char *name1, const T_y1 &y1, const char *name2, const T_y2 &y2)
Check if two structures at the same size.
void initialize(double *&mem_d, const double *ind, double *shared=nullptr)
double adj_
The adjoint of this variable, which is the partial derivative of this variable with respect to the ro...
Definition: vari.hpp:44
static double var_dot(const Eigen::DenseBase< Derived1 > &v1, const Eigen::DenseBase< Derived2 > &v2)
Definition: dot_product.hpp:57
void * alloc(size_t len)
Return a newly allocated block of memory of the appropriate size managed by the stack allocator...
fvar< T > inv(const fvar< T > &x)
Definition: inv.hpp:12

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