Stan  1.3
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
dot_product.hpp
Go to the documentation of this file.
1 #ifndef __STAN__AGRAD__REV__MATRIX__DOT_PRODUCT_HPP__
2 #define __STAN__AGRAD__REV__MATRIX__DOT_PRODUCT_HPP__
3 
4 #include <vector>
5 #include <stan/math/matrix/Eigen.hpp>
6 #include <stan/math/matrix/typedefs.hpp>
7 #include <stan/math/matrix/validate_vector.hpp>
8 #include <stan/math/matrix/validate_matching_sizes.hpp>
9 #include <stan/agrad/rev/var.hpp>
10 #include <stan/agrad/rev/vari.hpp>
11 #include <stan/agrad/rev/matrix/typedefs.hpp>
12 
13 namespace stan {
14  namespace agrad {
15 
16  namespace {
17  class dot_product_vv_vari : public vari {
18  protected:
19  vari** v1_;
20  vari** v2_;
21  size_t length_;
22  inline static double var_dot(const var* v1, const var* v2,
23  size_t length) {
24  double result = 0;
25  for (size_t i = 0; i < length; i++)
26  result += v1[i].vi_->val_ * v2[i].vi_->val_;
27  return result;
28  }
29  template<typename Derived1,typename Derived2>
30  inline static double var_dot(const Eigen::DenseBase<Derived1> &v1,
31  const Eigen::DenseBase<Derived2> &v2) {
32  double result = 0;
33  for (int i = 0; i < v1.size(); i++)
34  result += v1[i].vi_->val_ * v2[i].vi_->val_;
35  return result;
36  }
37  inline static double var_dot(vari** v1, vari** v2, size_t length) {
38  double result = 0;
39  for (size_t i = 0; i < length; ++i)
40  result += v1[i]->val_ * v2[i]->val_;
41  return result;
42  }
43  public:
44  dot_product_vv_vari(vari** v1, vari** v2, size_t length)
45  : vari(var_dot(v1,v2,length)),
46  v1_(v1),
47  v2_(v2),
48  length_(length) {
49 
50  }
51  dot_product_vv_vari(const var* v1, const var* v2, size_t length,
52  dot_product_vv_vari* shared_v1 = NULL,
53  dot_product_vv_vari* shared_v2 = NULL) :
54  vari(var_dot(v1, v2, length)), length_(length) {
55  if (shared_v1 == NULL) {
56  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
57  for (size_t i = 0; i < length_; i++)
58  v1_[i] = v1[i].vi_;
59  }
60  else {
61  v1_ = shared_v1->v1_;
62  }
63  if (shared_v2 == NULL) {
64  v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
65  for (size_t i = 0; i < length_; i++)
66  v2_[i] = v2[i].vi_;
67  }
68  else {
69  v2_ = shared_v2->v2_;
70  }
71  }
72  template<typename Derived1,typename Derived2>
73  dot_product_vv_vari(const Eigen::DenseBase<Derived1> &v1,
74  const Eigen::DenseBase<Derived2> &v2,
75  dot_product_vv_vari* shared_v1 = NULL,
76  dot_product_vv_vari* shared_v2 = NULL) :
77  vari(var_dot(v1, v2)), length_(v1.size()) {
78  if (shared_v1 == NULL) {
79  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
80  for (size_t i = 0; i < length_; i++)
81  v1_[i] = v1[i].vi_;
82  }
83  else {
84  v1_ = shared_v1->v1_;
85  }
86  if (shared_v2 == NULL) {
87  v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
88  for (size_t i = 0; i < length_; i++)
89  v2_[i] = v2[i].vi_;
90  }
91  else {
92  v2_ = shared_v2->v2_;
93  }
94  }
95  template<int R1,int C1,int R2,int C2>
96  dot_product_vv_vari(const Eigen::Matrix<var,R1,C1> &v1,
97  const Eigen::Matrix<var,R2,C2> &v2,
98  dot_product_vv_vari* shared_v1 = NULL,
99  dot_product_vv_vari* shared_v2 = NULL) :
100  vari(var_dot(v1, v2)), length_(v1.size()) {
101  if (shared_v1 == NULL) {
102  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
103  for (size_t i = 0; i < length_; i++)
104  v1_[i] = v1[i].vi_;
105  }
106  else {
107  v1_ = shared_v1->v1_;
108  }
109  if (shared_v2 == NULL) {
110  v2_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
111  for (size_t i = 0; i < length_; i++)
112  v2_[i] = v2[i].vi_;
113  }
114  else {
115  v2_ = shared_v2->v2_;
116  }
117  }
118  virtual void chain() {
119  for (size_t i = 0; i < length_; i++) {
120  v1_[i]->adj_ += adj_ * v2_[i]->val_;
121  v2_[i]->adj_ += adj_ * v1_[i]->val_;
122  }
123  }
124  };
125 
126  class dot_product_vd_vari : public vari {
127  protected:
128  vari** v1_;
129  double* v2_;
130  size_t length_;
131  inline static double var_dot(const var* v1, const double* v2,
132  size_t length) {
133  double result = 0;
134  for (size_t i = 0; i < length; i++)
135  result += v1[i].vi_->val_ * v2[i];
136  return result;
137  }
138  template<typename Derived1,typename Derived2>
139  inline static double var_dot(const Eigen::DenseBase<Derived1> &v1,
140  const Eigen::DenseBase<Derived2> &v2) {
141  double result = 0;
142  for (int i = 0; i < v1.size(); i++)
143  result += v1[i].vi_->val_ * v2[i];
144  return result;
145  }
146  public:
147  dot_product_vd_vari(const var* v1, const double* v2, size_t length,
148  dot_product_vd_vari *shared_v1 = NULL,
149  dot_product_vd_vari *shared_v2 = NULL) :
150  vari(var_dot(v1, v2, length)), length_(length) {
151  if (shared_v1 == NULL) {
152  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
153  for (size_t i = 0; i < length_; i++)
154  v1_[i] = v1[i].vi_;
155  } else {
156  v1_ = shared_v1->v1_;
157  }
158  if (shared_v2 == NULL) {
159  v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
160  for (size_t i = 0; i < length_; i++)
161  v2_[i] = v2[i];
162  } else {
163  v2_ = shared_v2->v2_;
164  }
165  }
166  template<typename Derived1,typename Derived2>
167  dot_product_vd_vari(const Eigen::DenseBase<Derived1> &v1,
168  const Eigen::DenseBase<Derived2> &v2,
169  dot_product_vd_vari *shared_v1 = NULL,
170  dot_product_vd_vari *shared_v2 = NULL) :
171  vari(var_dot(v1, v2)), length_(v1.size()) {
172  if (shared_v1 == NULL) {
173  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
174  for (size_t i = 0; i < length_; i++)
175  v1_[i] = v1[i].vi_;
176  } else {
177  v1_ = shared_v1->v1_;
178  }
179  if (shared_v2 == NULL) {
180  v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
181  for (size_t i = 0; i < length_; i++)
182  v2_[i] = v2[i];
183  } else {
184  v2_ = shared_v2->v2_;
185  }
186  }
187  template<int R1,int C1,int R2,int C2>
188  dot_product_vd_vari(const Eigen::Matrix<var,R1,C1> &v1,
189  const Eigen::Matrix<double,R2,C2> &v2,
190  dot_product_vd_vari *shared_v1 = NULL,
191  dot_product_vd_vari *shared_v2 = NULL) :
192  vari(var_dot(v1, v2)), length_(v1.size()) {
193  if (shared_v1 == NULL) {
194  v1_ = (vari**)memalloc_.alloc(length_*sizeof(vari*));
195  for (size_t i = 0; i < length_; i++)
196  v1_[i] = v1[i].vi_;
197  } else {
198  v1_ = shared_v1->v1_;
199  }
200  if (shared_v2 == NULL) {
201  v2_ = (double*)memalloc_.alloc(length_*sizeof(double));
202  for (size_t i = 0; i < length_; i++)
203  v2_[i] = v2[i];
204  } else {
205  v2_ = shared_v2->v2_;
206  }
207  }
208  virtual void chain() {
209  for (size_t i = 0; i < length_; i++) {
210  v1_[i]->adj_ += adj_ * v2_[i];
211  }
212  }
213  };
214  }
215 
224  template<int R1,int C1,int R2, int C2>
225  inline var dot_product(const Eigen::Matrix<var, R1, C1>& v1,
226  const Eigen::Matrix<var, R2, C2>& v2) {
227  stan::math::validate_vector(v1,"dot_product");
228  stan::math::validate_vector(v2,"dot_product");
229  stan::math::validate_matching_sizes(v1,v2,"dot_product");
230  return var(new dot_product_vv_vari(v1,v2));
231  }
241  template<int R1,int C1,int R2, int C2>
242  inline var dot_product(const Eigen::Matrix<var, R1, C1>& v1,
243  const Eigen::Matrix<double, R2, C2>& v2) {
244  stan::math::validate_vector(v1,"dot_product");
245  stan::math::validate_vector(v2,"dot_product");
246  stan::math::validate_matching_sizes(v1,v2,"dot_product");
247  return var(new dot_product_vd_vari(v1,v2));
248  }
258  template<int R1,int C1,int R2, int C2>
259  inline var dot_product(const Eigen::Matrix<double, R1, C1>& v1,
260  const Eigen::Matrix<var, R2, C2>& v2) {
261  stan::math::validate_vector(v1,"dot_product");
262  stan::math::validate_vector(v2,"dot_product");
263  stan::math::validate_matching_sizes(v1,v2,"dot_product");
264  return var(new dot_product_vd_vari(v2,v1));
265  }
274  inline var dot_product(const var* v1, const var* v2, size_t length) {
275  return var(new dot_product_vv_vari(v1, v2, length));
276  }
285  inline var dot_product(const var* v1, const double* v2, size_t length) {
286  return var(new dot_product_vd_vari(v1, v2, length));
287  }
296  inline var dot_product(const double* v1, const var* v2, size_t length) {
297  return var(new dot_product_vd_vari(v2, v1, length));
298  }
307  inline var dot_product(const std::vector<var>& v1,
308  const std::vector<var>& v2) {
309  stan::math::validate_matching_sizes(v1,v2,"dot_product");
310  return var(new dot_product_vv_vari(&v1[0], &v2[0], v1.size()));
311  }
320  inline var dot_product(const std::vector<var>& v1,
321  const std::vector<double>& v2) {
322  stan::math::validate_matching_sizes(v1,v2,"dot_product");
323  return var(new dot_product_vd_vari(&v1[0], &v2[0], v1.size()));
324  }
333  inline var dot_product(const std::vector<double>& v1,
334  const std::vector<var>& v2) {
335  stan::math::validate_matching_sizes(v1,v2,"dot_product");
336  return var(new dot_product_vd_vari(&v2[0], &v1[0], v1.size()));
337  }
338 
339  template<int R1,int C1,int R2, int C2>
340  inline Eigen::Matrix<var, 1, C1>
341  columns_dot_product(const Eigen::Matrix<var, R1, C1>& v1,
342  const Eigen::Matrix<var, R2, C2>& v2) {
343  stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
344  Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
345  for (size_type j = 0; j < v1.cols(); ++j) {
346  ret(j) = var(new dot_product_vv_vari(v1.col(j),v2.col(j)));
347  }
348  return ret;
349  }
350 
351  template<int R1,int C1,int R2, int C2>
352  inline Eigen::Matrix<var, 1, C1>
353  columns_dot_product(const Eigen::Matrix<var, R1, C1>& v1,
354  const Eigen::Matrix<double, R2, C2>& v2) {
355  stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
356  Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
357  for (size_type j = 0; j < v1.cols(); ++j) {
358  ret(j) = var(new dot_product_vd_vari(v1.col(j),v2.col(j)));
359  }
360  return ret;
361  }
362 
363  template<int R1,int C1,int R2, int C2>
364  inline Eigen::Matrix<var, 1, C1>
365  columns_dot_product(const Eigen::Matrix<double, R1, C1>& v1,
366  const Eigen::Matrix<var, R2, C2>& v2) {
367  stan::math::validate_matching_sizes(v1,v2,"columns_dot_product");
368  Eigen::Matrix<var, 1, C1> ret(1,v1.cols());
369  for (size_type j = 0; j < v1.cols(); ++j) {
370  ret(j) = var(new dot_product_vd_vari(v2.col(j),v1.col(j)));
371  }
372  return ret;
373  }
374 
375  template<int R1,int C1,int R2, int C2>
376  inline Eigen::Matrix<var, R1, 1>
377  rows_dot_product(const Eigen::Matrix<var, R1, C1>& v1,
378  const Eigen::Matrix<var, R2, C2>& v2) {
379  stan::math::validate_matching_sizes(v1,v2,"rows_dot_product");
380  Eigen::Matrix<var, R1, 1> ret(v1.rows(),1);
381  for (size_type j = 0; j < v1.rows(); ++j) {
382  ret(j) = var(new dot_product_vv_vari(v1.row(j),v2.row(j)));
383  }
384  return ret;
385  }
386 
387  template<int R1,int C1,int R2, int C2>
388  inline Eigen::Matrix<var, R1, 1>
389  rows_dot_product(const Eigen::Matrix<var, R1, C1>& v1,
390  const Eigen::Matrix<double, R2, C2>& v2) {
391  stan::math::validate_matching_sizes(v1,v2,"rows_dot_product");
392  Eigen::Matrix<var, R1, 1> ret(v1.rows(),1);
393  for (size_type j = 0; j < v1.rows(); ++j) {
394  ret(j) = var(new dot_product_vd_vari(v1.row(j),v2.row(j)));
395  }
396  return ret;
397  }
398 
399  template<int R1,int C1,int R2, int C2>
400  inline Eigen::Matrix<var, R1, 1>
401  rows_dot_product(const Eigen::Matrix<double, R1, C1>& v1,
402  const Eigen::Matrix<var, R2, C2>& v2) {
403  stan::math::validate_matching_sizes(v1,v2,"rows_dot_product");
404  Eigen::Matrix<var, R1, 1> ret(v1.rows(),1);
405  for (size_type j = 0; j < v1.rows(); ++j) {
406  ret(j) = var(new dot_product_vd_vari(v2.row(j),v1.row(j)));
407  }
408  return ret;
409  }
410  }
411 }
412 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.