Stan Math Library  2.10.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 
9 #include <stan/math/rev/core.hpp>
12 #include <boost/utility/enable_if.hpp>
13 #include <boost/type_traits.hpp>
14 #include <vector>
15 
16 namespace stan {
17  namespace math {
18 
19  namespace {
20  template<typename T>
21  struct dot_product_store_type;
22 
23  template<>
24  struct dot_product_store_type<var> {
25  typedef vari** type;
26  };
27 
28  template<>
29  struct dot_product_store_type<double> {
30  typedef double* type;
31  };
32 
33  template<typename T1, typename T2>
34  class dot_product_vari : public vari {
35  protected:
36  typename dot_product_store_type<T1>::type v1_;
37  typename dot_product_store_type<T2>::type v2_;
38  size_t length_;
39 
40  inline static double var_dot(vari** v1, vari** v2,
41  size_t length) {
42  Eigen::VectorXd vd1(length), vd2(length);
43  for (size_t i = 0; i < length; i++) {
44  vd1[i] = v1[i]->val_;
45  vd2[i] = v2[i]->val_;
46  }
47  return vd1.dot(vd2);
48  }
49 
50  inline static double var_dot(const T1* v1, const T2* v2,
51  size_t length) {
53  Eigen::VectorXd vd1(length), vd2(length);
54  for (size_t i = 0; i < length; i++) {
55  vd1[i] = value_of(v1[i]);
56  vd2[i] = value_of(v2[i]);
57  }
58  return vd1.dot(vd2);
59  }
60 
61  template<typename Derived1, typename Derived2>
62  inline static double var_dot(const Eigen::DenseBase<Derived1> &v1,
63  const Eigen::DenseBase<Derived2> &v2) {
66  Eigen::VectorXd vd1(v1.size()), vd2(v1.size());
67  for (int i = 0; i < v1.size(); i++) {
68  vd1[i] = value_of(v1[i]);
69  vd2[i] = value_of(v2[i]);
70  }
71  return vd1.dot(vd2);
72  }
73  inline void chain(vari** v1, vari** v2) {
74  for (size_t i = 0; i < length_; i++) {
75  v1[i]->adj_ += adj_ * v2_[i]->val_;
76  v2[i]->adj_ += adj_ * v1_[i]->val_;
77  }
78  }
79  inline void chain(double* v1, vari** v2) {
80  for (size_t i = 0; i < length_; i++) {
81  v2[i]->adj_ += adj_ * v1_[i];
82  }
83  }
84  inline void chain(vari** v1, double* v2) {
85  for (size_t i = 0; i < length_; i++) {
86  v1[i]->adj_ += adj_ * v2_[i];
87  }
88  }
89  inline void initialize(vari** &mem_v, const var *inv,
90  vari **shared = NULL) {
91  if (shared == NULL) {
92  mem_v = reinterpret_cast<vari**>(ChainableStack::memalloc_
93  .alloc(length_*sizeof(vari*)));
94  for (size_t i = 0; i < length_; i++)
95  mem_v[i] = inv[i].vi_;
96  } else {
97  mem_v = shared;
98  }
99  }
100  template<typename Derived>
101  inline void initialize(vari** &mem_v,
102  const Eigen::DenseBase<Derived> &inv,
103  vari **shared = NULL) {
104  if (shared == NULL) {
105  mem_v = reinterpret_cast<vari**>(ChainableStack::memalloc_
106  .alloc(length_*sizeof(vari*)));
107  for (size_t i = 0; i < length_; i++)
108  mem_v[i] = inv(i).vi_;
109  } else {
110  mem_v = shared;
111  }
112  }
113 
114  inline void initialize(double* &mem_d, const double *ind,
115  double *shared = NULL) {
116  if (shared == NULL) {
117  mem_d = reinterpret_cast<double*>(ChainableStack::memalloc_
118  .alloc(length_*sizeof(double)));
119  for (size_t i = 0; i < length_; i++)
120  mem_d[i] = ind[i];
121  } else {
122  mem_d = shared;
123  }
124  }
125  template<typename Derived>
126  inline void initialize(double* &mem_d,
127  const Eigen::DenseBase<Derived> &ind,
128  double *shared = NULL) {
129  if (shared == NULL) {
130  mem_d = reinterpret_cast<double*>
131  (ChainableStack::memalloc_.alloc(length_*sizeof(double)));
132  for (size_t i = 0; i < length_; i++)
133  mem_d[i] = ind(i);
134  } else {
135  mem_d = shared;
136  }
137  }
138 
139  public:
140  dot_product_vari(typename dot_product_store_type<T1>::type v1,
141  typename dot_product_store_type<T2>::type v2,
142  size_t length)
143  : vari(var_dot(v1, v2, length)), v1_(v1), v2_(v2), length_(length) {}
144 
145  dot_product_vari(const T1* v1, const T2* v2, size_t length,
146  dot_product_vari<T1, T2>* shared_v1 = NULL,
147  dot_product_vari<T1, T2>* shared_v2 = NULL) :
148  vari(var_dot(v1, v2, length)), length_(length) {
149  if (shared_v1 == NULL) {
150  initialize(v1_, v1);
151  } else {
152  initialize(v1_, v1, shared_v1->v1_);
153  }
154  if (shared_v2 == NULL) {
155  initialize(v2_, v2);
156  } else {
157  initialize(v2_, v2, shared_v2->v2_);
158  }
159  }
160  template<typename Derived1, typename Derived2>
161  dot_product_vari(const Eigen::DenseBase<Derived1> &v1,
162  const Eigen::DenseBase<Derived2> &v2,
163  dot_product_vari<T1, T2>* shared_v1 = NULL,
164  dot_product_vari<T1, T2>* shared_v2 = NULL) :
165  vari(var_dot(v1, v2)), length_(v1.size()) {
166  if (shared_v1 == NULL) {
167  initialize(v1_, v1);
168  } else {
169  initialize(v1_, v1, shared_v1->v1_);
170  }
171  if (shared_v2 == NULL) {
172  initialize(v2_, v2);
173  } else {
174  initialize(v2_, v2, shared_v2->v2_);
175  }
176  }
177  template<int R1, int C1, int R2, int C2>
178  dot_product_vari(const Eigen::Matrix<T1, R1, C1> &v1,
179  const Eigen::Matrix<T2, R2, C2> &v2,
180  dot_product_vari<T1, T2>* shared_v1 = NULL,
181  dot_product_vari<T1, T2>* shared_v2 = NULL) :
182  vari(var_dot(v1, v2)), length_(v1.size()) {
183  if (shared_v1 == NULL) {
184  initialize(v1_, v1);
185  } else {
186  initialize(v1_, v1, shared_v1->v1_);
187  }
188  if (shared_v2 == NULL) {
189  initialize(v2_, v2);
190  } else {
191  initialize(v2_, v2, shared_v2->v2_);
192  }
193  }
194  virtual void chain() {
195  chain(v1_, v2_);
196  }
197  };
198  }
199 
208  template<typename T1, int R1, int C1, typename T2, int R2, int C2>
209  inline
210  typename boost::enable_if_c<boost::is_same<T1, var>::value ||
211  boost::is_same<T2, var>::value, var>::type
212  dot_product(const Eigen::Matrix<T1, R1, C1>& v1,
213  const Eigen::Matrix<T2, R2, C2>& v2) {
214  stan::math::check_vector("dot_product", "v1", v1);
215  stan::math::check_vector("dot_product", "v2", v2);
216  stan::math::check_matching_sizes("dot_product",
217  "v1", v1,
218  "v2", v2);
219  return var(new dot_product_vari<T1, T2>(v1, v2));
220  }
229  template<typename T1, typename T2>
230  inline
231  typename boost::enable_if_c<boost::is_same<T1, var>::value ||
232  boost::is_same<T2, var>::value, var>::type
233  dot_product(const T1* v1, const T2* v2, size_t length) {
234  return var(new dot_product_vari<T1, T2>(v1, v2, length));
235  }
236 
245  template<typename T1, typename T2>
246  inline
247  typename boost::enable_if_c<boost::is_same<T1, var>::value ||
248  boost::is_same<T2, var>::value, var>::type
249  dot_product(const std::vector<T1>& v1,
250  const std::vector<T2>& v2) {
251  stan::math::check_matching_sizes("dot_product",
252  "v1", v1,
253  "v2", v2);
254  return var(new dot_product_vari<T1, T2>(&v1[0], &v2[0], v1.size()));
255  }
256 
257  }
258 }
259 #endif
bool check_vector(const char *function, const char *name, const Eigen::Matrix< T, R, C > &x)
Return true if the matrix is either a row vector or column vector.
T value_of(const fvar< T > &v)
Return the value of the specified variable.
Definition: value_of.hpp:16
dot_product_store_type< T1 >::type v1_
Definition: dot_product.hpp:36
size_t length(const std::vector< T > &x)
Definition: length.hpp:10
Independent (input) and dependent (output) variables for gradients.
Definition: var.hpp:31
dot_product_store_type< T2 >::type v2_
Definition: dot_product.hpp:37
void initialize(T &x, const T &v)
Definition: initialize.hpp:17
bool check_matching_sizes(const char *function, const char *name1, const T_y1 &y1, const char *name2, const T_y2 &y2)
Return true if two structures at the same size.
fvar< T > dot_product(const Eigen::Matrix< fvar< T >, R1, C1 > &v1, const Eigen::Matrix< fvar< T >, R2, C2 > &v2)
Definition: dot_product.hpp:20
int size(const std::vector< T > &x)
Return the size of the specified standard vector.
Definition: size.hpp:17
size_t length_
Definition: dot_product.hpp:38
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:15

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