Stan  1.1
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
matrix.hpp
Go to the documentation of this file.
1 #ifndef __STAN__MATH__MATRIX_HPP__
2 #define __STAN__MATH__MATRIX_HPP__
3 
4 #include <stdarg.h>
5 #include <stdexcept>
6 #include <ostream>
7 #include <vector>
8 
9 #include <boost/math/tools/promotion.hpp>
10 
11 #define EIGEN_DENSEBASE_PLUGIN "stan/math/EigenDenseBaseAddons.hpp"
12 #include <Eigen/Dense>
13 #include <Eigen/QR>
14 
15 #include <stan/math/boost_error_handling.hpp>
16 
17 namespace stan {
18 
19  namespace math {
20 
21  // from input type F to output type T
22 
23  // scalar, F != T (base template)
24  template <typename F, typename T>
25  struct promoter {
26  inline static void promote(const F& u, T& t) {
27  t = u;
28  }
29  inline static T promote_to(const F& u) {
30  return u;
31  }
32  };
33  // scalar, F == T
34  template <typename T>
35  struct promoter<T,T> {
36  inline static void promote(const T& u, T& t) {
37  t = u;
38  }
39  inline static T promote_to(const T& u) {
40  return u;
41  }
42  };
43 
44  // std::vector, F != T
45  template <typename F, typename T>
46  struct promoter<std::vector<F>, std::vector<T> > {
47  inline static void promote(const std::vector<F>& u,
48  std::vector<T>& t) {
49  t.resize(u.size());
50  for (size_t i = 0; i < u.size(); ++i)
51  promoter<F,T>::promote(u[i],t[i]);
52  }
53  inline static std::vector<T>
54  promote_to(const std::vector<F>& u) {
55  std::vector<T> t;
56  promoter<std::vector<F>,std::vector<T> >::promote(u,t);
57  return t;
58  }
59  };
60  // std::vector, F == T
61  template <typename T>
62  struct promoter<std::vector<T>, std::vector<T> > {
63  inline static void promote(const std::vector<T>& u,
64  std::vector<T>& t) {
65  t = u;
66  }
67  inline static std::vector<T> promote_to(const std::vector<T>& u) {
68  return u;
69  }
70  };
71 
72  // Eigen::Matrix, F != T
73  template <typename F, typename T, int R, int C>
74  struct promoter<Eigen::Matrix<F,R,C>, Eigen::Matrix<T,R,C> > {
75  inline static void promote(const Eigen::Matrix<F,R,C>& u,
76  Eigen::Matrix<T,R,C>& t) {
77  t.resize(u.rows(), u.cols());
78  for (int i = 0; i < u.size(); ++i)
79  promoter<F,T>::promote(u(i),t(i));
80  }
81  inline static Eigen::Matrix<T,R,C>
82  promote_to(const Eigen::Matrix<F,R,C>& u) {
83  Eigen::Matrix<T,R,C> t;
84  promoter<Eigen::Matrix<F,R,C>,Eigen::Matrix<T,R,C> >::promote(u,t);
85  return t;
86  }
87  };
88  // Eigen::Matrix, F == T
89  template <typename T, int R, int C>
90  struct promoter<Eigen::Matrix<T,R,C>, Eigen::Matrix<T,R,C> > {
91  inline static void promote(const Eigen::Matrix<T,R,C>& u,
92  Eigen::Matrix<T,R,C>& t) {
93  t = u;
94  }
95  inline static Eigen::Matrix<T,R,C> promote_to(const Eigen::Matrix<T,R,C>& u) {
96  return u;
97  }
98  };
99 
100  template <typename T1, typename T2>
101  struct common_type {
102  typedef typename boost::math::tools::promote_args<T1,T2>::type type;
103  };
104 
105  template <typename T1, typename T2>
106  struct common_type<std::vector<T1>, std::vector<T2> > {
107  typedef std::vector<typename common_type<T1,T2>::type> type;
108  };
109 
110  template <typename T1, typename T2, int R, int C>
111  struct common_type<Eigen::Matrix<T1,R,C>, Eigen::Matrix<T2,R,C> > {
112  typedef Eigen::Matrix<typename common_type<T1,T2>::type,R,C> type;
113  };
114 
115  template <typename T1, typename T2, typename F>
116  inline
117  typename common_type<T1,T2>::type
118  promote_common(const F& u) {
119  return promoter<F, typename common_type<T1,T2>::type>
120  ::promote_to(u);
121  }
122 
123 
124 
125 
131  template <typename T>
132  struct array_builder {
133  std::vector<T> x_;
134  array_builder() : x_() { }
135  template <typename F>
136  array_builder& add(const F& u) {
137  T t;
138  promoter<F,T>::promote(u,t);
139  x_.push_back(t);
140  return *this;
141  }
142  std::vector<T> array() {
143  return x_;
144  }
145  };
146 
150  typedef
151  Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic>
152  matrix_d;
153 
157  typedef
158  Eigen::Matrix<double,Eigen::Dynamic,1>
159  vector_d;
160 
164  typedef
165  Eigen::Matrix<double,1,Eigen::Dynamic>
166  row_vector_d;
167 
168  namespace {
169 
170  template <typename T>
171  void resize(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& x,
172  const std::vector<size_t>& dims,
173  size_t pos) {
174  x.resize(dims[pos],dims[pos+1]);
175  }
176 
177  template <typename T>
178  void resize(Eigen::Matrix<T,Eigen::Dynamic,1>& x,
179  const std::vector<size_t>& dims,
180  size_t pos) {
181  x.resize(dims[pos]);
182  }
183 
184  template <typename T>
185  void resize(Eigen::Matrix<T,1,Eigen::Dynamic>& x,
186  const std::vector<size_t>& dims,
187  size_t pos) {
188  x.resize(dims[pos]);
189  }
190 
191 
192  void resize(double /*x*/,
193  const std::vector<size_t>& /*dims*/,
194  size_t /*pos*/) {
195  // no-op
196  }
197 
198  template <typename T>
199  void resize(std::vector<T>& x,
200  const std::vector<size_t>& dims,
201  size_t pos) {
202  x.resize(dims[pos]);
203  ++pos;
204  if (pos >= dims.size()) return; // skips lowest loop to scalar
205  for (size_t i = 0; i < x.size(); ++i)
206  resize(x[i],dims,pos);
207  }
208 
209  }
210 
220  template <typename T>
221  inline void resize(T& x, std::vector<size_t> dims) {
222  resize(x,dims,0U);
223  }
224 
225  // polymorphic gets with bounds checking
226 
227 
228  namespace {
229 
230  void raise_range_error(size_t max,
231  size_t i,
232  const char* msg,
233  size_t idx) {
234  std::stringstream s;
235  s << "INDEX OPERATOR [] OUT OF BOUNDS"
236  << "; index=" << i
237  << "; lower bound=1"
238  << "; upper bound=" << max
239  << "; index position=" << idx
240  << "; " << msg
241  << std::endl;
242  throw std::out_of_range(s.str());
243  }
244 
245  inline
246  void check_range(size_t max,
247  size_t i,
248  const char* msg,
249  size_t idx) {
250 #ifndef NDEBUG
251  if (i < 1 || i > max)
252  raise_range_error(max,i,msg,idx);
253 #endif
254  }
255 
256  }
257 
272  template <typename T>
273  inline
274  T& get_base1(std::vector<T>& x,
275  size_t i,
276  const char* error_msg,
277  size_t idx) {
278  check_range(x.size(),i,error_msg,idx);
279  return x[i - 1];
280  }
281 
297  template <typename T>
298  inline
299  T& get_base1(std::vector<std::vector<T> >& x,
300  size_t i1,
301  size_t i2,
302  const char* error_msg,
303  size_t idx) {
304  check_range(x.size(),i1,error_msg,idx);
305  return get_base1(x[i1 - 1],i2,error_msg,idx+1);
306  }
307 
324  template <typename T>
325  inline
326  T& get_base1(std::vector<std::vector<std::vector<T> > >& x,
327  size_t i1,
328  size_t i2,
329  size_t i3,
330  const char* error_msg,
331  size_t idx) {
332  check_range(x.size(),i1,error_msg,idx);
333  return get_base1(x[i1 - 1],i2,i3,error_msg,idx+1);
334  }
335 
353  template <typename T>
354  inline
355  T& get_base1(std::vector<std::vector<std::vector<std::vector<T> > > >& x,
356  size_t i1,
357  size_t i2,
358  size_t i3,
359  size_t i4,
360  const char* error_msg,
361  size_t idx) {
362  check_range(x.size(),i1,error_msg,idx);
363  return get_base1(x[i1 - 1],i2,i3,i4,error_msg,idx+1);
364  }
365 
384  template <typename T>
385  inline
386  T& get_base1(std::vector<std::vector<std::vector<std::vector<std::vector<T> > > > >& x,
387  size_t i1,
388  size_t i2,
389  size_t i3,
390  size_t i4,
391  size_t i5,
392  const char* error_msg,
393  size_t idx) {
394  check_range(x.size(),i1,error_msg,idx);
395  return get_base1(x[i1 - 1],i2,i3,i4,i5,error_msg,idx+1);
396  }
397 
417  template <typename T>
418  inline
419  T& get_base1(std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<T> > > > > >& x,
420  size_t i1,
421  size_t i2,
422  size_t i3,
423  size_t i4,
424  size_t i5,
425  size_t i6,
426  const char* error_msg,
427  size_t idx) {
428  check_range(x.size(),i1,error_msg,idx);
429  return get_base1(x[i1 - 1],i2,i3,i4,i5,i6,error_msg,idx+1);
430  }
431 
432 
453  template <typename T>
454  inline
455  T& get_base1(std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<T> > > > > > >& x,
456  size_t i1,
457  size_t i2,
458  size_t i3,
459  size_t i4,
460  size_t i5,
461  size_t i6,
462  size_t i7,
463  const char* error_msg,
464  size_t idx) {
465  check_range(x.size(),i1,error_msg,idx);
466  return get_base1(x[i1 - 1],i2,i3,i4,i5,i6,i7,error_msg,idx+1);
467  }
468 
469 
491  template <typename T>
492  inline
493  T& get_base1(std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<std::vector<T> > > > > > > >& x,
494  size_t i1,
495  size_t i2,
496  size_t i3,
497  size_t i4,
498  size_t i5,
499  size_t i6,
500  size_t i7,
501  size_t i8,
502  const char* error_msg,
503  size_t idx) {
504  check_range(x.size(),i1,error_msg,idx);
505  return get_base1(x[i1 - 1],i2,i3,i4,i5,i6,i7,i8,error_msg,idx+1);
506  }
507 
508 
509 
529  template <typename T>
530  inline
531  Eigen::Matrix<T,1,Eigen::Dynamic>
532  get_base1(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& x,
533  size_t m,
534  const char* error_msg,
535  size_t idx) {
536  check_range(x.rows(),m,error_msg,idx);
537  return x.row(m - 1);
538  }
539 
556  template <typename T>
557  inline
558  T& get_base1(Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& x,
559  size_t m,
560  size_t n,
561  const char* error_msg,
562  size_t idx) {
563  check_range(x.rows(),m,error_msg,idx);
564  check_range(x.cols(),n,error_msg,idx + 1);
565  return x(m - 1, n - 1);
566  }
567 
582  template <typename T>
583  inline
584  T& get_base1(Eigen::Matrix<T,Eigen::Dynamic,1>& x,
585  size_t m,
586  const char* error_msg,
587  size_t idx) {
588  check_range(x.size(),m,error_msg,idx);
589  return x(m - 1);
590 
591  }
592 
607  template <typename T>
608  inline
609  T& get_base1(Eigen::Matrix<T,1,Eigen::Dynamic>& x,
610  size_t n,
611  const char* error_msg,
612  size_t idx) {
613  check_range(x.size(),n,error_msg,idx);
614  return x(n - 1);
615  }
616 
617 
618  // int returns
619 
620  template <typename T, int R, int C>
621  inline
622  size_t
623  rows(const Eigen::Matrix<T,R,C>& m) {
624  return m.rows();
625  }
626  template <typename T, int R, int C>
627  inline
628  size_t
629  cols(const Eigen::Matrix<T,R,C>& m) {
630  return m.cols();
631  }
632 
633  template <typename T1, typename T2>
634  inline
635  void validate_less_or_equal(const T1& x, const T2& y,
636  const char* x_name, const char* y_name,
637  const char* fun_name) {
638  if (x <= y) return;
639  std::stringstream ss;
640  ss << "require " << x_name << " <= " << y_name
641  << " in " << fun_name
642  << "; found " << x_name << "=" << x
643  << ", " << y_name << "=" << y;
644  throw std::domain_error(ss.str());
645  }
646  template <typename T1, typename T2>
647  inline
648  void validate_less(const T1& x, const T2& y,
649  const char* x_name, const char* y_name,
650  const char* fun_name) {
651  if (x < y) return;
652  std::stringstream ss;
653  ss << "require " << x_name << " < " << y_name
654  << " in " << fun_name
655  << "; found " << x_name << "=" << x
656  << ", " << y_name << "=" << y;
657  throw std::domain_error(ss.str());
658  }
659  template <typename T1, typename T2>
660  inline
661  void validate_greater_or_equal(const T1& x, const T2& y,
662  const char* x_name, const char* y_name,
663  const char* fun_name) {
664  if (x >= y) return;
665  std::stringstream ss;
666  ss << "require " << x_name << " >= " << y_name
667  << " in " << fun_name
668  << "; found " << x_name << "=" << x
669  << ", " << y_name << "=" << y;
670  throw std::domain_error(ss.str());
671  }
672  template <typename T1, typename T2>
673  inline
674  void validate_greater(const T1& x, const T2& y,
675  const char* x_name, const char* y_name,
676  const char* fun_name) {
677  if (x > y) return;
678  std::stringstream ss;
679  ss << "require " << x_name << " > " << y_name
680  << " in " << fun_name
681  << "; found " << x_name << "=" << x
682  << ", " << y_name << "=" << y;
683  throw std::domain_error(ss.str());
684  }
685 
686 
687 
688  template <typename T, int R, int C>
689  void validate_column_index(const Eigen::Matrix<T,R,C>& m,
690  size_t j,
691  const char* msg) {
692  if (j > 0 && j <= static_cast<size_t>(m.cols())) return;
693  std::stringstream ss;
694  ss << "require 0 < column index <= number of columns in" << msg;
695  ss << " found cols()=" << m.cols()
696  << "; index j=" << j;
697  throw std::domain_error(ss.str());
698  }
699 
700  template <typename T, int R, int C>
701  void validate_row_index(const Eigen::Matrix<T,R,C>& m,
702  size_t i,
703  const char* msg) {
704  if (i > 0 && i <= static_cast<size_t>(m.rows())) return;
705  std::stringstream ss;
706  ss << "require 0 < row index <= number of rows in" << msg;
707  ss << " found rows()=" << m.rows()
708  << "; index i=" << i;
709  throw std::domain_error(ss.str());
710  }
711 
712  template <typename T, int R, int C>
713  void validate_square(const Eigen::Matrix<T,R,C>& x,
714  const char* msg) {
715  if (x.rows() == x.cols()) return;
716  std::stringstream ss;
717  ss << "error in call to " << msg
718  << "; require square matrix, but found"
719  << " rows=" << x.rows()
720  << "; cols=" << x.cols();
721  throw std::domain_error(ss.str());
722  }
723 
724  template <typename T, int R, int C>
725  void validate_symmetric(const Eigen::Matrix<T,R,C>& x,
726  const char* msg) {
727  // tolerance = 1E-8
728  validate_square(x,msg);
729  for (int i = 0; i < x.rows(); ++i) {
730  for (int j = 0; j < x.cols(); ++j) {
731  if (x(i,j) != x(j,i)) {
732  std::stringstream ss;
733  ss << "error in call to " << msg
734  << "; require symmetric matrix, but found"
735  << "; x[" << i << "," << j << "]=" << x(i,j)
736  << "; x[" << j << "," << i << "]=" << x(j,i);
737  throw std::domain_error(ss.str());
738  }
739  }
740  }
741  }
742 
743  template <typename T1, int R1, int C1, typename T2, int R2, int C2>
744  inline void validate_matching_dims(const Eigen::Matrix<T1,R1,C1>& x1,
745  const Eigen::Matrix<T2,R2,C2>& x2,
746  const char* msg) {
747  if (x1.rows() == x2.rows()
748  && x1.cols() == x2.cols()) return;
749  std::stringstream ss;
750  ss << "error in call to " << msg
751  << "; require matching dimensions, but found"
752  << " arg1(rows=" << x1.rows() << ",cols=" << x1.cols() << ");"
753  << " arg2(rows=" << x2.rows() << ",cols=" << x2.cols() << ")";
754  throw std::domain_error(ss.str());
755  }
756 
757  template <typename T1, typename T2>
758  inline void validate_matching_sizes(const std::vector<T1>& x1,
759  const std::vector<T2>& x2,
760  const char* msg) {
761  if (x1.size() == x2.size()) return;
762  std::stringstream ss;
763  ss << "error in call to " << msg
764  << "; require matching sizes, but found"
765  << " arg1(size=" << x1.size() << ");"
766  << " arg2(size=" << x2.size() << ");";
767  throw std::domain_error(ss.str());
768  }
769 
770  template <typename T1, int R1, int C1, typename T2, int R2, int C2>
771  inline void validate_matching_sizes(const Eigen::Matrix<T1,R1,C1>& x1,
772  const Eigen::Matrix<T2,R2,C2>& x2,
773  const char* msg) {
774  if (x1.size() == x2.size()) return;
775  std::stringstream ss;
776  ss << "error in call to " << msg
777  << "; require matching sizes, but found"
778  << " arg1(rows=" << x1.rows() << ",cols=" << x1.cols()
779  << ",size=" << (x1.rows() * x1.cols()) << ");"
780  << " arg2(rows=" << x2.rows() << ",cols=" << x2.cols()
781  << ",size=" << (x2.rows() * x2.cols()) << ")";
782  throw std::domain_error(ss.str());
783  }
784 
785  template <typename T1, int R1, int C1, typename T2, int R2, int C2>
786  inline void validate_multiplicable(const Eigen::Matrix<T1,R1,C1>& x1,
787  const Eigen::Matrix<T2,R2,C2>& x2,
788  const char* msg) {
789  if (x1.cols() == x2.rows()) return;
790  std::stringstream ss;
791  ss << "error in call to " << msg
792  << "; require cols of arg1 to match rows of arg2, but found "
793  << " arg1 rows=" << x1.rows() << " arg1 cols=" << x1.cols()
794  << " arg2 rows=" << x2.rows() << " arg2 cols=" << x2.cols();
795  throw std::domain_error(ss.str());
796  }
797 
798  template <typename T>
799  inline void validate_nonzero_size(const T& x, const char* msg) {
800  if (x.size() > 0) return;
801  std::stringstream ss;
802  ss << "require non-zero size for " << msg
803  << "found size=" << x.size();
804  throw std::domain_error(ss.str());
805  }
806 
807  template <typename T, int R, int C>
808  inline void validate_vector(const Eigen::Matrix<T,R,C>& x,
809  const char* msg) {
810  if (x.rows() == 1 || x.cols() == 1) return;
811  std::stringstream ss;
812  ss << "error in " << msg
813  << "; require vector, found "
814  << " rows=" << x.rows() << "cols=" << x.cols();
815  throw std::domain_error(ss.str());
816  }
817 
818  // scalar returns
819 
827  template <typename T,int R, int C>
828  inline T determinant(const Eigen::Matrix<T,R,C>& m) {
829  stan::math::validate_square(m,"determinant");
830  return m.determinant();
831  }
832 
840  template <typename T,int R, int C>
841  inline T log_determinant(const Eigen::Matrix<T,R,C>& m) {
842  stan::math::validate_square(m,"log_determinant");
843  return m.colPivHouseholderQr().logAbsDeterminant();
844  }
845 
846 
854  template <int R, int C>
855  inline double dot_self(const Eigen::Matrix<double, R, C>& v) {
856  validate_vector(v,"dot_self");
857  double sum = 0.0;
858  for (int i = 0; i < v.size(); ++i)
859  sum += v(i) * v(i);
860  return sum;
861  }
862 
868  template<typename T,int R,int C>
869  inline Eigen::Matrix<T,1,C>
870  columns_dot_self(const Eigen::Matrix<T,R,C>& x) {
871  return x.colwise().squaredNorm();
872  }
873 
883  template<int R1,int C1,int R2, int C2>
884  inline Eigen::Matrix<double, 1, C1>
885  columns_dot_product(const Eigen::Matrix<double, R1, C1>& v1,
886  const Eigen::Matrix<double, R2, C2>& v2) {
887  validate_matching_sizes(v1,v2,"columns_dot_product");
888  Eigen::Matrix<double, 1, C1> ret(1,v1.cols());
889  size_t j;
890  for (j = 0; j < v1.cols(); ++j) {
891  ret(j) = v1.col(j).dot(v2.col(j));
892  }
893  return ret;
894  }
895 
905  template<int R1,int C1,int R2, int C2>
906  inline double dot_product(const Eigen::Matrix<double, R1, C1>& v1,
907  const Eigen::Matrix<double, R2, C2>& v2) {
908  validate_vector(v1,"dot_product");
909  validate_vector(v2,"dot_product");
910  validate_matching_sizes(v1,v2,"dot_product");
911  double sum = 0.0;
912  for (int i = 0; i < v1.size(); ++i)
913  sum += v1[i] * v2[i];
914  return sum;
915  }
922  inline double dot_product(const double* v1, const double* v2,
923  size_t length) {
924  double result = 0;
925  for (size_t i = 0; i < length; i++)
926  result += v1[i] * v2[i];
927  return result;
928  }
935  inline double dot_product(const std::vector<double>& v1,
936  const std::vector<double>& v2) {
937  validate_matching_sizes(v1,v2,"dot_product");
938  return dot_product(&v1[0], &v2[0], v1.size());
939  }
940 
948  inline int min(const std::vector<int>& x) {
949  if (x.size() == 0)
950  throw std::domain_error("error: cannot take min of empty int vector");
951  int min = x[0];
952  for (size_t i = 1; i < x.size(); ++i)
953  if (x[i] < min)
954  min = x[i];
955  return min;
956  }
957 
965  template <typename T>
966  inline T min(const std::vector<T>& x) {
967  if (x.size() == 0)
968  return std::numeric_limits<T>::infinity();
969  T min = x[0];
970  for (size_t i = 1; i < x.size(); ++i)
971  if (x[i] < min)
972  min = x[i];
973  return min;
974  }
975 
982  template <typename T, int R, int C>
983  inline T min(const Eigen::Matrix<T,R,C>& m) {
984  if (m.size() == 0)
985  return std::numeric_limits<double>::infinity();
986  return m.minCoeff();
987  }
988 
989 
990 
999  inline int max(const std::vector<int>& x) {
1000  if (x.size() == 0)
1001  throw std::domain_error("error: cannot take max of empty int vector");
1002  int max = x[0];
1003  for (size_t i = 1; i < x.size(); ++i)
1004  if (x[i] > max)
1005  max = x[i];
1006  return max;
1007  }
1008 
1016  template <typename T>
1017  inline T max(const std::vector<T>& x) {
1018  if (x.size() == 0)
1019  return -std::numeric_limits<T>::infinity();
1020  T max = x[0];
1021  for (size_t i = 1; i < x.size(); ++i)
1022  if (x[i] > max)
1023  max = x[i];
1024  return max;
1025  }
1026 
1033  template <typename T, int R, int C>
1034  inline T max(const Eigen::Matrix<T,R,C>& m) {
1035  if (m.size() == 0)
1036  return -std::numeric_limits<double>::infinity();
1037  return m.maxCoeff();
1038  }
1039 
1048  template <typename T>
1049  inline
1050  typename boost::math::tools::promote_args<T>::type
1051  mean(const std::vector<T>& v) {
1052  validate_nonzero_size(v,"mean");
1053  T sum(v[0]);
1054  for (size_t i = 1; i < v.size(); ++i)
1055  sum += v[i];
1056  return sum / v.size();
1057  }
1058 
1065  template <typename T, int R, int C>
1066  inline
1067  typename boost::math::tools::promote_args<T>::type
1068  mean(const Eigen::Matrix<T,R,C>& m) {
1069  validate_nonzero_size(m,"mean");
1070  return m.mean();
1071  }
1072 
1081  template <typename T>
1082  inline
1083  typename boost::math::tools::promote_args<T>::type
1084  variance(const std::vector<T>& v) {
1085  validate_nonzero_size(v,"variance");
1086  if (v.size() == 1)
1087  return 0.0;
1088  T v_mean(mean(v));
1089  T sum_sq_diff(0);
1090  for (size_t i = 0; i < v.size(); ++i) {
1091  T diff = v[i] - v_mean;
1092  sum_sq_diff += diff * diff;
1093  }
1094  return sum_sq_diff / (v.size() - 1);
1095  }
1096 
1103  template <typename T, int R, int C>
1104  inline
1105  typename boost::math::tools::promote_args<T>::type
1106  variance(const Eigen::Matrix<T,R,C>& m) {
1107  validate_nonzero_size(m,"variance");
1108  if (m.size() == 1)
1109  return 0.0;
1110  typename boost::math::tools::promote_args<T>::type
1111  mn(mean(m));
1112  typename boost::math::tools::promote_args<T>::type
1113  sum_sq_diff(0);
1114  for (int i = 0; i < m.size(); ++i) {
1115  typename boost::math::tools::promote_args<T>::type
1116  diff = m(i) - mn;
1117  sum_sq_diff += diff * diff;
1118  }
1119  return sum_sq_diff / (m.size() - 1);
1120  }
1121 
1128  template <typename T>
1129  inline
1130  typename boost::math::tools::promote_args<T>::type
1131  sd(const std::vector<T>& v) {
1132  validate_nonzero_size(v,"sd");
1133  if (v.size() == 1) return 0.0;
1134  return sqrt(variance(v));
1135  }
1136 
1143  template <typename T, int R, int C>
1144  inline
1145  typename boost::math::tools::promote_args<T>::type
1146  sd(const Eigen::Matrix<T,R,C>& m) {
1147  // FIXME: redundant with test in variance; second line saves sqrt
1148  validate_nonzero_size(m,"sd");
1149  if (m.size() == 1) return 0.0;
1150  return sqrt(variance(m));
1151  }
1152 
1161  template <typename T>
1162  inline T sum(const std::vector<T>& xs) {
1163  if (xs.size() == 0) return 0;
1164  T sum(xs[0]);
1165  for (size_t i = 1; i < xs.size(); ++i)
1166  sum += xs[i];
1167  return sum;
1168  }
1169 
1176  template <typename T, int R, int C>
1177  inline double sum(const Eigen::Matrix<T,R,C>& v) {
1178  return v.sum();
1179  }
1180 
1187  template <typename T>
1188  inline T prod(const std::vector<T>& v) {
1189  if (v.size() == 0) return 1;
1190  T product = v[0];
1191  for (size_t i = 1; i < v.size(); ++i)
1192  product *= v[i];
1193  return product;
1194  }
1195 
1202  template <typename T, int R, int C>
1203  inline T prod(const Eigen::Matrix<T,R,C>& v) {
1204  if (v.size() == 0) return 1.0;
1205  return v.prod();
1206  }
1207 
1217  template <typename T>
1218  inline T trace(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
1219  return m.trace();
1220  }
1221 
1228  template<typename T, int Rows, int Cols>
1229  inline Eigen::Matrix<T,Rows,Cols> log(const Eigen::Matrix<T,Rows,Cols>& m) {
1230  return m.array().log().matrix();
1231  }
1232 
1239  template<typename T, int Rows, int Cols>
1240  inline Eigen::Matrix<T,Rows,Cols> exp(const Eigen::Matrix<T,Rows,Cols>& m) {
1241  return m.array().exp().matrix();
1242  }
1243 
1244 
1245  // vector and matrix returns
1246 
1261  template <typename T1, typename T2, int R, int C>
1262  inline
1263  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R,C>
1264  add(const Eigen::Matrix<T1,R,C>& m1,
1265  const Eigen::Matrix<T2,R,C>& m2) {
1266  stan::math::validate_matching_dims(m1,m2,"add");
1267  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R,C>
1268  result(m1.rows(),m1.cols());
1269  for (int i = 0; i < result.size(); ++i)
1270  result(i) = m1(i) + m2(i);
1271  return result;
1272  }
1282  template <typename T1, typename T2, int R, int C>
1283  inline
1284  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R,C>
1285  add(const Eigen::Matrix<T1,R,C>& m,
1286  const T2& c) {
1287  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R,C>
1288  result(m.rows(),m.cols());
1289  for (int i = 0; i < result.size(); ++i)
1290  result(i) = m(i) + c;
1291  return result;
1292  }
1302  template <typename T1, typename T2, int R, int C>
1303  inline
1304  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R,C>
1305  add(const T1& c,
1306  const Eigen::Matrix<T2,R,C>& m) {
1307  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R,C>
1308  result(m.rows(),m.cols());
1309  for (int i = 0; i < result.size(); ++i)
1310  result(i) = c + m(i);
1311  return result;
1312  }
1313 
1314 
1328  template <typename T1, typename T2, int R, int C>
1329  inline
1330  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1331  subtract(const Eigen::Matrix<T1,R,C>& m1,
1332  const Eigen::Matrix<T2,R,C>& m2) {
1333  stan::math::validate_matching_dims(m1,m2,"subtract");
1334  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1335  result(m1.rows(), m1.cols());
1336  for (int i = 0; i < result.size(); ++i)
1337  result(i) = m1(i) - m2(i);
1338  return result;
1339  }
1340 
1341  template <typename T1, typename T2, int R, int C>
1342  inline
1343  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1344  subtract(const T1& c,
1345  const Eigen::Matrix<T2,R,C>& m) {
1346 
1347  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1348  result(m.rows(),m.cols());
1349  for (int i = 0; i < m.size(); ++i)
1350  result(i) = c - m(i);
1351  return result;
1352  }
1353  template <typename T1, typename T2, int R, int C>
1354  inline
1355  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1356  subtract(const Eigen::Matrix<T1,R,C>& m,
1357  const T2& c) {
1358 
1359  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1360  result(m.rows(),m.cols());
1361  for (int i = 0; i < m.size(); ++i)
1362  result(i) = m(i) - c;
1363  return result;
1364  }
1365 
1366 
1367 
1375  template <typename T>
1376  inline
1377  T minus(const T& x) {
1378  return -x;
1379  }
1380 
1381 
1390  template <int R, int C>
1391  inline
1392  Eigen::Matrix<double,R,C>
1393  divide(const Eigen::Matrix<double,R,C>& m,
1394  double c) {
1395  return m / c;
1396  }
1397 
1410  template <typename T1, typename T2, int R, int C>
1411  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1412  elt_multiply(const Eigen::Matrix<T1,R,C>& m1,
1413  const Eigen::Matrix<T2,R,C>& m2) {
1414  stan::math::validate_matching_dims(m1,m2,"elt_multiply");
1415  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1416  result(m1.rows(),m2.cols());
1417  for (int i = 0; i < m1.size(); ++i)
1418  result(i) = m1(i) * m2(i);
1419  return result;
1420  }
1421 
1422  template <typename T1, typename T2, int R1,int C1,int R2,int C2>
1423  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R2, C2>
1424  diag_pre_multiply(const Eigen::Matrix<T1,R1,C1>& m1,
1425  const Eigen::Matrix<T2,R2,C2>& m2) {
1426  if (m1.cols() != 1 && m1.rows() != 1)
1427  throw std::domain_error("m1 must be a vector");
1428  if (m1.size() != m2.rows())
1429  throw std::domain_error("m1 must have same length as m2 has rows");
1430  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R2, C2>
1431  result(m2.rows(),m2.cols());
1432  for (int i = 0; i < m2.rows(); ++i)
1433  for (int j = 0; i < m2.cols(); ++j)
1434  result(i,j) = m1(i) * m2(i,j);
1435  return result;
1436  }
1437 
1438  template <typename T1, typename T2, int R1,int C1,int R2,int C2>
1439  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R1, C1>
1440  diag_post_multiply(const Eigen::Matrix<T1,R1,C1>& m1,
1441  const Eigen::Matrix<T2,R2,C2>& m2) {
1442  if (m2.cols() != 1 && m2.rows() != 1)
1443  throw std::domain_error("m2 must be a vector");
1444  if (m2.size() != m1.cols())
1445  throw std::domain_error("m2 must have same length as m1 has columns");
1446  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R1, C1>
1447  result(m1.rows(),m1.cols());
1448  for (int i = 0; i < m1.rows(); ++i)
1449  for (int j = 0; j < m1.cols(); ++j)
1450  result(i,j) = m2(i) * m1(i,j);
1451  return result;
1452  }
1453 
1454 
1467  template <typename T1, typename T2, int R, int C>
1468  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1469  elt_divide(const Eigen::Matrix<T1,R,C>& m1,
1470  const Eigen::Matrix<T2,R,C>& m2) {
1471  stan::math::validate_matching_dims(m1,m2,"elt_multiply");
1472  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type, R, C>
1473  result(m1.rows(),m2.cols());
1474  for (int i = 0; i < m1.size(); ++i)
1475  result(i) = m1(i) / m2(i);
1476  return result;
1477  }
1478 
1479 
1488  template <int R, int C>
1489  inline
1490  Eigen::Matrix<double,R,C>
1491  multiply(const Eigen::Matrix<double,R,C>& m,
1492  double c) {
1493  return c * m;
1494  }
1503  template <int R, int C>
1504  inline
1505  Eigen::Matrix<double,R,C>
1506  multiply(double c,
1507  const Eigen::Matrix<double,R,C>& m) {
1508  return c * m;
1509  }
1510 
1511 
1522  template<int R1,int C1,int R2,int C2>
1523  inline Eigen::Matrix<double,R1,C2> multiply(const Eigen::Matrix<double,R1,C1>& m1,
1524  const Eigen::Matrix<double,R2,C2>& m2) {
1525 
1526  validate_multiplicable(m1,m2,"multiply");
1527  return m1*m2;
1528  }
1538  template<int C1,int R2>
1539  inline double multiply(const Eigen::Matrix<double,1,C1>& rv,
1540  const Eigen::Matrix<double,R2,1>& v) {
1541  stan::math::validate_matching_sizes(rv,v,"multiply");
1542  if (rv.size() != v.size())
1543  throw std::domain_error("rv.size() != v.size()");
1544  return rv.dot(v);
1545  }
1546 
1555  inline matrix_d
1556  multiply_lower_tri_self_transpose(const matrix_d& L) {
1557  if (L.rows() == 0)
1558  return matrix_d(0,0);
1559  if (L.rows() == 1) {
1560  matrix_d result(1,1);
1561  result(0,0) = L(0,0) * L(0,0);
1562  return result;
1563  }
1564  // FIXME: write custom following agrad/matrix because can't get L_tri into
1565  // multiplication as no template support for tri * tri
1566  matrix_d L_tri = L.transpose().triangularView<Eigen::Upper>();
1567  return L.triangularView<Eigen::Lower>() * L_tri;
1568  }
1569 
1576  inline matrix_d
1577  tcrossprod(const matrix_d& M) {
1578  if (M.rows() == 0)
1579  return matrix_d(0,0);
1580  if (M.rows() == 1)
1581  return M * M.transpose();
1582  matrix_d result(M.rows(),M.rows());
1583  return result
1584  .setZero()
1585  .selfadjointView<Eigen::Upper>()
1586  .rankUpdate(M);
1587  }
1588 
1595  inline matrix_d
1596  crossprod(const matrix_d& M) {
1597  return tcrossprod(M.transpose());
1598  }
1599 
1612  template <typename T>
1613  inline
1614  Eigen::Matrix<T,1,Eigen::Dynamic>
1615  row(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m,
1616  size_t i) {
1617  validate_row_index(m,i,"row");
1618  return m.row(i - 1);
1619  }
1620 
1632  template <typename T>
1633  inline
1634  Eigen::Matrix<T,Eigen::Dynamic,1>
1635  col(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m,
1636  size_t j) {
1637  validate_column_index(m,j,"col");
1638  return m.col(j - 1);
1639  }
1640 
1650  template <typename T>
1651  inline
1652  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1653  block(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m,
1654  size_t i, size_t j, size_t nrows, size_t ncols) {
1655  validate_row_index(m,i,"block");
1656  validate_row_index(m,i+nrows-1,"block");
1657  validate_column_index(m,j,"block");
1658  validate_column_index(m,j+ncols-1,"block");
1659  return m.block(i - 1,j - 1,nrows,ncols);
1660  }
1661 
1662 
1669  template <typename T>
1670  inline
1671  Eigen::Matrix<T,Eigen::Dynamic,1>
1672  diagonal(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
1673  return m.diagonal();
1674  }
1675 
1682  template <typename T>
1683  inline
1684  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1685  diag_matrix(const Eigen::Matrix<T,Eigen::Dynamic,1>& v) {
1686  return v.asDiagonal();
1687  }
1688 
1689  template <typename T, int R, int C>
1690  Eigen::Matrix<T,C,R>
1691  inline
1692  transpose(const Eigen::Matrix<T,R,C>& m) {
1693  return m.transpose();
1694  }
1695 
1701  template <typename T>
1702  inline
1703  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1704  inverse(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
1705  validate_square(m,"matrix inverse");
1706  return m.inverse();
1707  }
1708 
1709 
1716  template <typename T>
1717  inline
1718  Eigen::Matrix<T,Eigen::Dynamic,1>
1719  softmax(const Eigen::Matrix<T,Eigen::Dynamic,1>& v) {
1720  using std::exp;
1721  stan::math::validate_nonzero_size(v,"vector softmax");
1722  Eigen::Matrix<T,Eigen::Dynamic,1> theta(v.size());
1723  T sum(0.0);
1724  T max_v = v.maxCoeff();
1725  for (int i = 0; i < v.size(); ++i) {
1726  theta[i] = exp(v[i] - max_v);
1727  sum += theta[i];
1728  }
1729  for (int i = 0; i < v.size(); ++i)
1730  theta[i] /= sum;
1731  return theta;
1732  }
1733 
1743  template <int TriView, typename T1, typename T2,
1744  int R1, int C1, int R2, int C2>
1745  inline
1746  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,
1747  R1,C2>
1748  mdivide_left_tri(const Eigen::Matrix<T1,R1,C1> &A,
1749  const Eigen::Matrix<T2,R2,C2> &b) {
1750  stan::math::validate_square(A,"mdivide_left_tri");
1751  stan::math::validate_multiplicable(A,b,"mdivide_left_tri");
1752  return promote_common<Eigen::Matrix<T1,R1,C1>,Eigen::Matrix<T2,R1,C1> >(A)
1753  .template triangularView<TriView>()
1754  .solve( promote_common<Eigen::Matrix<T1,R2,C2>,
1755  Eigen::Matrix<T2,R2,C2> >(b) );
1756  }
1764  template<int TriView, typename T,int R1, int C1>
1765  inline
1766  Eigen::Matrix<T,R1,C1>
1767  mdivide_left_tri(const Eigen::Matrix<T,R1,C1> &A) {
1768  stan::math::validate_square(A,"mdivide_left_tri");
1769  int n = A.rows();
1770  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> b;
1771  b.setIdentity(n,n);
1772  A.template triangularView<TriView>().solveInPlace(b);
1773  return b;
1774  }
1775  template <typename T1, typename T2, int R1,int C1,int R2,int C2>
1776  inline
1777  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,
1778  R1,C2>
1779  mdivide_left_tri_low(const Eigen::Matrix<T1,R1,C1> &A,
1780  const Eigen::Matrix<T2,R2,C2> &b) {
1781 // stan::math::validate_square(A,"mdivide_left_tri_low/2");
1782 // stan::math::validate_multiplicable(A,b,"mdivide_left_tri_low");
1783 // return promote_common<Eigen::Matrix<T1,R1,C1>,
1784 // Eigen::Matrix<T2,R1,C1> >(A)
1785 // .template triangularView<Eigen::Lower>()
1786 // .solve( promote_common<Eigen::Matrix<T1,R2,C2>,
1787 // Eigen::Matrix<T2,R2,C2> >(b) );
1788  return mdivide_left_tri<Eigen::Lower,T1,T2,R1,C1,R2,C2>(A,b);
1789  }
1790  template <typename T,int R1, int C1>
1791  inline
1792  Eigen::Matrix<T,R1,C1>
1793  mdivide_left_tri_low(const
1794  Eigen::Matrix<T,R1,C1> &A) {
1795 // stan::math::validate_square(A,"mdivide_left_tri_low/1");
1796 // int n = A.rows();
1797 // Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> b;
1798 // b.setIdentity(n,n);
1799 // A.template triangularView<Eigen::Lower>().solveInPlace(b);
1800 // return b;
1801  return mdivide_left_tri<Eigen::Lower,T,R1,C1>(A);
1802  }
1803 
1812  template <typename T1, typename T2, int R1, int C1, int R2, int C2>
1813  inline
1814  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R1,C2>
1815  mdivide_left(const Eigen::Matrix<T1,R1,C1> &A,
1816  const Eigen::Matrix<T2,R2,C2> &b) {
1817  stan::math::validate_square(A,"mdivide_left");
1818  stan::math::validate_multiplicable(A,b,"mdivide_left");
1819  return promote_common<Eigen::Matrix<T1,R1,C1>,
1820  Eigen::Matrix<T2,R1,C1> >(A)
1821  .lu()
1822  .solve( promote_common<Eigen::Matrix<T1,R2,C2>,
1823  Eigen::Matrix<T2,R2,C2> >(b) );
1824  }
1825 
1835  template <int TriView, typename T1, typename T2,
1836  int R1, int C1, int R2, int C2>
1837  inline
1838  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R1,C2>
1839  mdivide_right_tri(const Eigen::Matrix<T1,R1,C1> &b,
1840  const Eigen::Matrix<T2,R2,C2> &A) {
1841  stan::math::validate_square(A,"mdivide_left_tri_low");
1842  stan::math::validate_multiplicable(b,A,"mdivide_right_tri");
1843  // FIXME: This is nice and general but requires some extra memory and copying.
1844  if (TriView == Eigen::Lower) {
1845  return transpose(mdivide_left_tri<Eigen::Upper>(transpose(A),transpose(b)));
1846  }
1847  else if (TriView == Eigen::Upper) {
1848  return transpose(mdivide_left_tri<Eigen::Lower>(transpose(A),transpose(b)));
1849  }
1850  else {
1851  assert(TriView == Eigen::Lower || TriView == Eigen::Upper);
1852  }
1853 
1854 // return promote_common<Eigen::Matrix<T1,R1,C1>,
1855 // Eigen::Matrix<T2,R1,C1> >(A)
1856 // .template triangularView<TriView>()
1857 // .transpose()
1858 // .solve(promote_common<Eigen::Matrix<T1,R2,C2>,
1859 // Eigen::Matrix<T2,R2,C2> >(b)
1860 // .transpose())
1861 // .transpose();
1862  }
1872  template <typename T1, typename T2, int R1,int C1,int R2,int C2>
1873  inline
1874  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R1,C2>
1875  mdivide_right_tri_low(const Eigen::Matrix<T1,R1,C1> &b,
1876  const Eigen::Matrix<T2,R2,C2> &A) {
1877  return mdivide_right_tri<Eigen::Lower>
1878  (promote_common<Eigen::Matrix<T1,R1,C1>,
1879  Eigen::Matrix<T2,R1,C1> >(b),
1880  promote_common<Eigen::Matrix<T1,R2,C2>,
1881  Eigen::Matrix<T2,R2,C2> >(A));
1882  }
1883 
1884 
1885 
1894  template <typename T1, typename T2, int R1, int C1, int R2, int C2>
1895  inline
1896  Eigen::Matrix<typename boost::math::tools::promote_args<T1,T2>::type,R1,C2>
1897  mdivide_right(const Eigen::Matrix<T1,R1,C1> &b,
1898  const Eigen::Matrix<T2,R2,C2> &A) {
1899  stan::math::validate_square(A,"mdivide_right");
1900  stan::math::validate_multiplicable(b,A,"mdivide_right");
1901  // FIXME: This is nice and general but likely slow.
1902  return transpose(mdivide_left(transpose(A),transpose(b)));
1903 // return promote_common<Eigen::Matrix<T1,R2,C2>,
1904 // Eigen::Matrix<T2,R2,C2> >(A)
1905 // .transpose()
1906 // .lu()
1907 // .solve(promote_common<Eigen::Matrix<T1,R1,C1>,
1908 // Eigen::Matrix<T2,R1,C1> >(b)
1909 // .transpose())
1910 // .transpose();
1911  }
1912 
1922  template <typename T>
1923  Eigen::Matrix<T,Eigen::Dynamic,1>
1924  eigenvalues_sym(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
1925  validate_nonzero_size(m,"eigenvalues_sym");
1926  validate_symmetric(m,"eigenvalues_sym");
1927  Eigen::SelfAdjointEigenSolver<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >
1928  solver(m,Eigen::EigenvaluesOnly);
1929  return solver.eigenvalues();
1930  }
1931 
1932  template <typename T>
1933  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1934  eigenvectors_sym(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
1935  validate_nonzero_size(m,"eigenvectors_sym");
1936  validate_symmetric(m,"eigenvectors_sym");
1937  Eigen::SelfAdjointEigenSolver<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >
1938  solver(m);
1939  return solver.eigenvectors();
1940  }
1941 
1942 
1943 
1944 
1955  template <typename T>
1956  Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>
1957  cholesky_decompose(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
1958  validate_symmetric(m,"cholesky decomposition");
1959  Eigen::LLT<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >llt(m.rows());
1960  llt.compute(m);
1961  return llt.matrixL();
1962  }
1963 
1972  template <typename T>
1973  Eigen::Matrix<T,Eigen::Dynamic,1>
1974  singular_values(const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& m) {
1975  return Eigen::JacobiSVD<Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic> >(m)
1976  .singularValues();
1977  }
1978 
1979  // void eigen_decompose_sym(const matrix_d& m,
1980  // vector_d& eigenvalues,
1981  // matrix_d& eigenvectors) {
1982  // Eigen::SelfAdjointEigenSolver<matrix_d> solver(m);
1983  // eigenvalues = solver.eigenvalues().real();
1984  // eigenvectors = solver.eigenvectors().real();
1985  // }
1986 
1987 
1988  // void svd(const matrix_d& m, matrix_d& u, matrix_d& v, vector_d& s) {
1989  // static const unsigned int THIN_SVD_OPTIONS
1990  // = Eigen::ComputeThinU | Eigen::ComputeThinV;
1991  // Eigen::JacobiSVD<matrix_d> svd(m, THIN_SVD_OPTIONS);
1992  // u = svd.matrixU();
1993  // v = svd.matrixV();
1994  // s = svd.singularValues();
1995  // }
1996 
2008  template <typename T>
2009  inline
2010  std::vector<T>
2011  cumulative_sum(const std::vector<T>& x) {
2012  std::vector<T> result(x.size());
2013  if (x.size() == 0)
2014  return result;
2015  result[0] = x[0];
2016  for (size_t i = 1; i < result.size(); ++i)
2017  result[i] = x[i] + result[i-1];
2018  return result;
2019  }
2034  template <typename T, int R, int C>
2035  inline
2036  Eigen::Matrix<T,R,C>
2037  cumulative_sum(const Eigen::Matrix<T,R,C>& m) {
2038  Eigen::Matrix<T,R,C> result(m.rows(),m.cols());
2039  if (m.size() == 0)
2040  return result;
2041  result(0) = m(0);
2042  for (int i = 1; i < result.size(); ++i)
2043  result(i) = m(i) + result(i-1);
2044  return result;
2045  }
2046 
2047 
2048  // prints used in generator for print() statements in modeling language
2049 
2050  template <typename T>
2051  void stan_print(std::ostream* o, const T& x) {
2052  *o << x;
2053  }
2054 
2055  template <typename T>
2056  void stan_print(std::ostream* o, const std::vector<T>& x) {
2057  *o << '[';
2058  for (int i = 0; i < x.size(); ++i) {
2059  if (i > 0) *o << ',';
2060  stan_print(o,x[i]);
2061  }
2062  *o << ']';
2063  }
2064 
2065  template <typename T>
2066  void stan_print(std::ostream* o, const Eigen::Matrix<T,Eigen::Dynamic,1>& x) {
2067  *o << '[';
2068  for (int i = 0; i < x.size(); ++i) {
2069  if (i > 0) *o << ',';
2070  stan_print(o,x(i));
2071  }
2072  *o << ']';
2073  }
2074 
2075  template <typename T>
2076  void stan_print(std::ostream* o, const Eigen::Matrix<T,1,Eigen::Dynamic>& x) {
2077  *o << '[';
2078  for (int i = 0; i < x.size(); ++i) {
2079  if (i > 0) *o << ',';
2080  stan_print(o,x(i));
2081  }
2082  *o << ']';
2083  }
2084 
2085  template <typename T>
2086  void stan_print(std::ostream* o,
2087  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& x) {
2088  *o << '[';
2089  for (int i = 0; i < x.rows(); ++i) {
2090  if (i > 0) *o << ',';
2091  stan_print(o,x.row(i));
2092  }
2093  *o << ']';
2094  }
2095 
2096  template <typename T>
2097  inline
2098  void
2099  dims(const T& x, std::vector<int>& result) {
2100  /* no op */
2101  }
2102  template <typename T, int R, int C>
2103  inline
2104  void
2105  dims(const Eigen::Matrix<T,R,C>& x,
2106  std::vector<int>& result) {
2107  result.push_back(x.rows());
2108  result.push_back(x.cols());
2109  }
2110  template <typename T>
2111  inline
2112  void
2113  dims(const std::vector<T>& x,
2114  std::vector<int>& result) {
2115  result.push_back(x.size());
2116  if (x.size() > 0)
2117  dims(x[0],result);
2118  }
2119 
2120 
2121  template <typename T>
2122  inline
2123  std::vector<int>
2124  dims(const T& x) {
2125  std::vector<int> result;
2126  dims(x,result);
2127  return result;
2128  }
2129 
2130  template <typename T>
2131  inline
2132  int
2133  size(const std::vector<T>& x) {
2134  return x.size();
2135  }
2136 
2137 
2138 
2139  }
2140 
2141 }
2142 
2143 #endif
2144 
     [ Stan Home Page ] © 2011–2013, Stan Development Team.