Stan  1.1
probability, sampling & optimization
 All Classes Namespaces Files Functions Variables Typedefs Enumerator Friends Macros Pages
matrix_error_handling.hpp
Go to the documentation of this file.
1 #ifndef __STAN__MATH__MATRIX_ERROR_HANDLING_HPP__
2 #define __STAN__MATH__MATRIX_ERROR_HANDLING_HPP__
3 
4 #include <sstream>
5 #include <stan/math/matrix.hpp>
6 #include <stan/meta/traits.hpp>
7 #include <stan/math/error_handling.hpp>
8 #include <boost/type_traits/common_type.hpp>
9 
10 namespace stan {
11 
12  namespace math {
13 
14  inline
15  void
16  validate_non_negative_index(const std::string& var_name,
17  const std::string& expr,
18  int val) {
19  if (val < 0) {
20  std::stringstream msg;
21  msg << "Found negative dimension size in variable declaration"
22  << "; variable=" << var_name
23  << "; dimension size expression=" << expr
24  << "; expression value=" << val;
25  throw std::invalid_argument(msg.str());
26  }
27  }
28 
29  // FIXME: update warnings
30  template <typename T_size1, typename T_size2, typename T_result,
31  class Policy>
32  inline bool check_size_match(const char* function,
33  T_size1 i,
34  const char* name_i,
35  T_size2 j,
36  const char* name_j,
37  T_result* result,
38  const Policy&) {
39  using stan::math::policies::raise_domain_error;
40  typedef typename boost::common_type<T_size1,T_size2>::type common_type;
41  if (static_cast<common_type>(i) != static_cast<common_type>(j)) {
42  std::ostringstream msg;
43  msg << name_i << " (%1%) and " << name_j << " (" << j << ") must match in size";
44  T_result tmp = policies::raise_domain_error<T_result,T_size1>(function,
45  msg.str().c_str(),
46  i, Policy());
47  if (result != 0)
48  *result = tmp;
49  return false;
50  }
51  return true;
52  }
53 
54  template <typename T_size1, typename T_size2, typename T_result>
55  inline
56  bool check_size_match(const char* function,
57  T_size1 i,
58  const char* name_i,
59  T_size2 j,
60  const char* name_j,
61  T_result* result) {
62  return check_size_match(function,i,name_i,j,name_j,result,default_policy());
63  }
64 
65  template <typename T_size1, typename T_size2>
66  inline
67  bool check_size_match(const char* function,
68  T_size1 i,
69  const char* name_i,
70  T_size2 j,
71  const char* name_j,
72  T_size1* result = 0) {
73  return check_size_match(function,i,name_i,j,name_j,result,default_policy());
74  }
75 
88  template <typename T_y, typename T_result, class Policy>
89  inline bool check_symmetric(const char* function,
90  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
91  const char* name,
92  T_result* result,
93  const Policy&) {
94  typedef
95  typename Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>::size_type
96  size_type;
97  size_type k = y.rows();
98  if (k == 1)
99  return true;
100  for (size_type m = 0; m < k; ++m) {
101  for (size_type n = m + 1; n < k; ++n) {
102  if (fabs(y(m,n) - y(n,m)) > CONSTRAINT_TOLERANCE) {
103  std::ostringstream message;
104  message << name << " is not symmetric. "
105  << name << "[" << m << "," << n << "] is %1%, but "
106  << name << "[" << n << "," << m
107  << "] element is " << y(n,m);
108  T_result tmp
109  = policies::raise_domain_error<T_y>(function,
110  message.str().c_str(),
111  y(m,n), Policy());
112  if (result != 0)
113  *result = tmp;
114  return false;
115  }
116  }
117  }
118  return true;
119  }
120 
121 
122  template <typename T_y, typename T_result>
123  inline bool check_symmetric(const char* function,
124  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
125  const char* name,
126  T_result* result) {
127  return check_symmetric(function,y,name,result,default_policy());
128  }
129 
130  template <typename T>
131  inline bool check_symmetric(const char* function,
132  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y,
133  const char* name,
134  T* result = 0) {
135  return check_symmetric(function,y,name,result,default_policy());
136  }
137 
138 
139 
140 
153  // FIXME: update warnings (message has (0,0) item)
154  template <typename T_y, typename T_result, class Policy>
155  inline bool check_pos_definite(const char* function,
156  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
157  const char* name,
158  T_result* result,
159  const Policy&) {
160  typedef
161  typename Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>::size_type
162  size_type;
163  if (y.rows() == 1 && y(0,0) <= CONSTRAINT_TOLERANCE) {
164  std::ostringstream message;
165  message << name << " is not positive definite. "
166  << name << "(0,0) is %1%.";
167  T_result tmp = policies::raise_domain_error<T_y>(function,
168  message.str().c_str(),
169  y(0,0), Policy());
170  if (result != 0)
171  *result = tmp;
172  return false;
173  }
174  Eigen::LDLT< Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic> > cholesky
175  = y.ldlt();
176  if((cholesky.vectorD().array() <= CONSTRAINT_TOLERANCE).any()) {
177  std::ostringstream message;
178  message << name << " is not positive definite. "
179  << name << "(0,0) is %1%.";
180  T_result tmp = policies::raise_domain_error<T_y>(function,
181  message.str().c_str(),
182  y(0,0), Policy());
183  if (result != 0)
184  *result = tmp;
185  return false;
186  }
187  return true;
188  }
189 
190 
191  template <typename T_y, typename T_result>
192  inline bool check_pos_definite(const char* function,
193  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
194  const char* name,
195  T_result* result) {
196  return check_pos_definite(function,y,name,result,default_policy());
197  }
198 
199 
200  template <typename T>
201  inline bool check_pos_definite(const char* function,
202  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y,
203  const char* name,
204  T* result = 0) {
205  return check_pos_definite(function,y,name,result,default_policy());
206  }
207 
208 
209 
210 
223  // FIXME: update warnings
224  template <typename T_y, typename T_result, class Policy>
225  inline bool check_cov_matrix(const char* function,
226  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
227  const char* name,
228  T_result* result,
229  const Policy&) {
230  if (!check_size_match(function,
231  y.rows(), "Rows of covariance matrix",
232  y.cols(), "columns of covariance matrix",
233  result, Policy()))
234  return false;
235  if (!check_positive(function, y.rows(), "rows", result, Policy()))
236  return false;
237  if (!check_symmetric(function, y, name, result, Policy()))
238  return false;
239  if (!check_pos_definite(function, y, name, result, Policy()))
240  return false;
241  return true;
242  }
243 
244  template <typename T_y, typename T_result>
245  inline bool check_cov_matrix(const char* function,
246  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
247  const char* name,
248  T_result* result) {
249  return check_cov_matrix(function,y,name,result,default_policy());
250  }
251 
252 
253  template <typename T>
254  inline bool check_cov_matrix(const char* function,
255  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y,
256  const char* name,
257  T* result = 0) {
258  return check_cov_matrix(function,y,name,result,default_policy());
259  }
260 
261 
262 
278  // FIXME: update warnings
279  template <typename T_y, typename T_result, class Policy>
280  inline bool check_corr_matrix(const char* function,
281  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
282  const char* name,
283  T_result* result,
284  const Policy&) {
285  if (!check_size_match(function,
286  y.rows(), "Rows of correlation matrix",
287  y.cols(), "columns of correlation matrix",
288  result, Policy()))
289  return false;
290  if (!check_positive(function, y.rows(), "rows", result, Policy()))
291  return false;
292  if (!check_symmetric(function, y, "y", result, Policy()))
293  return false;
294  for (typename Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>::size_type
295  k = 0; k < y.rows(); ++k) {
296  if (fabs(y(k,k) - 1.0) > CONSTRAINT_TOLERANCE) {
297  std::ostringstream message;
298  message << name << " is not a valid correlation matrix. "
299  << name << "(" << k << "," << k
300  << ") is %1%, but should be near 1.0";
301  T_result tmp
302  = policies::raise_domain_error<T_y>(function,
303  message.str().c_str(),
304  y(k,k), Policy());
305  if (result != 0)
306  *result = tmp;
307  return false;
308  }
309  }
310  if (!check_pos_definite(function, y, "y", result, Policy()))
311  return false;
312  return true;
313  }
314 
315 
316  template <typename T_y, typename T_result>
317  inline bool check_corr_matrix(const char* function,
318  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
319  const char* name,
320  T_result* result) {
321  return check_corr_matrix(function,y,name,result,default_policy());
322  }
323 
324  template <typename T>
325  inline bool check_corr_matrix(const char* function,
326  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y,
327  const char* name,
328  T* result = 0) {
329  return check_corr_matrix(function,y,name,result,default_policy());
330  }
331 
343  // FIXME: update warnings
344  template <typename T_covar, typename T_result, class Policy>
345  inline bool check_cov_matrix(const char* function,
346  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
347  T_result* result,
348  const Policy&) {
349  if (!check_size_match(function,
350  Sigma.rows(), "Rows of covariance matrix",
351  Sigma.cols(), "columns of covariance matrix",
352  result, Policy()))
353  return false;
354  if (!check_positive(function, Sigma.rows(), "rows", result, Policy()))
355  return false;
356  if (!check_symmetric(function, Sigma, "Sigma", result, Policy()))
357  return false;
358  return true;
359  }
360  template <typename T_covar, typename T_result>
361  inline bool check_cov_matrix(const char* function,
362  const Eigen::Matrix<T_covar,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
363  T_result* result) {
364  return check_cov_matrix(function,Sigma,result,default_policy());
365 
366  }
367  template <typename T>
368  inline bool check_cov_matrix(const char* function,
369  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& Sigma,
370  T* result = 0) {
371  return check_cov_matrix(function,Sigma,result,default_policy());
372 
373  }
374 
375 
376 
377 
392  template <typename T_prob,
393  typename T_result,
394  class Policy>
395  bool check_simplex(const char* function,
396  const Eigen::Matrix<T_prob,Eigen::Dynamic,1>& theta,
397  const char* name,
398  T_result* result,
399  const Policy&) {
400  typedef typename Eigen::Matrix<T_prob,Eigen::Dynamic,1>::size_type size_t;
401  using stan::math::policies::raise_domain_error;
402  if (theta.size() == 0) {
403  std::string message(name);
404  message += " is not a valid simplex. %1% elements in the vector.";
405  T_result tmp = raise_domain_error<size_t, size_t>(function,
406  message.c_str(),
407  0,
408  Policy());
409  if (result != 0)
410  *result = tmp;
411  return false;
412  }
413  if (fabs(1.0 - theta.sum()) > CONSTRAINT_TOLERANCE) {
414  std::stringstream msg;
415  T_prob sum = theta.sum();
416  msg << "in function check_simplex(%1%), ";
417  msg << name << " is not a valid simplex.";
418  msg << " The sum of the elements should be 1, but is " << sum;
419  T_result tmp = raise_domain_error<T_result,T_prob>(function,
420  msg.str().c_str(),
421  sum,
422  Policy());
423  if (result != 0)
424  *result = tmp;
425  return false;
426  }
427  for (size_t n = 0; n < theta.size(); n++) {
428  if (!(theta[n] >= 0)) {
429  std::ostringstream stream;
430  stream << name << " is not a valid simplex."
431  << " The element at " << n
432  << " is %1%, but should be greater than or equal to 0";
433  T_result tmp
434  = raise_domain_error<T_result,T_prob>(function,
435  stream.str().c_str(),
436  theta[n],
437  Policy());
438  if (result != 0)
439  *result = tmp;
440  return false;
441  }
442  }
443  return true;
444  }
445  template <typename T_y,
446  typename T_result> // = typename T_prob_vector::value_type,
447  inline bool check_simplex(const char* function,
448  const Eigen::Matrix<T_y,Eigen::Dynamic,1>& theta,
449  const char* name,
450  T_result* result) {
451  return check_simplex(function,theta,name,result,default_policy());
452  }
453  template <typename T>
454  inline bool check_simplex(const char* function,
455  const Eigen::Matrix<T,Eigen::Dynamic,1>& theta,
456  const char* name,
457  T* result = 0) {
458  return check_simplex(function,theta,name,result,default_policy());
459  }
460 
461 
462 
463 
478  template <typename T_y, typename T_result, class Policy>
479  bool check_ordered(const char* function,
480  const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
481  const char* name,
482  T_result* result,
483  const Policy&) {
484  using stan::math::policies::raise_domain_error;
485  typedef typename Eigen::Matrix<T_y,Eigen::Dynamic,1>::size_type size_t;
486  if (y.size() == 0) {
487  return true;
488  }
489  for (size_t n = 1; n < y.size(); n++) {
490  if (!(y[n] > y[n-1])) {
491  std::ostringstream stream;
492  stream << name << " is not a valid ordered vector."
493  << " The element at " << n
494  << " is %1%, but should be greater than the previous element, "
495  << y[n-1];
496  T_result tmp = raise_domain_error<T_result,T_y>(function,
497  stream.str().c_str(),
498  y[n],
499  Policy());
500  if (result != 0)
501  *result = tmp;
502  return false;
503  }
504  }
505  return true;
506  }
507  template <typename T_y, typename T_result>
508  bool check_ordered(const char* function,
509  const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
510  const char* name,
511  T_result* result) {
512  return check_ordered(function,y,name,result,default_policy());
513  }
514  template <typename T>
515  bool check_ordered(const char* function,
516  const Eigen::Matrix<T,Eigen::Dynamic,1>& y,
517  const char* name,
518  T* result = 0) {
519  return check_ordered(function,y,name,result,default_policy());
520  }
521 
536  template <typename T_y, typename T_result, class Policy>
537  bool check_positive_ordered(const char* function,
538  const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
539  const char* name,
540  T_result* result,
541  const Policy&) {
542  using stan::math::policies::raise_domain_error;
543  typedef typename Eigen::Matrix<T_y,Eigen::Dynamic,1>::size_type size_t;
544  if (y.size() == 0) {
545  return true;
546  }
547  if (y[0] < 0) {
548  std::ostringstream stream;
549  stream << name << " is not a valid positive_ordered vector."
550  << " The element at 0 is %1%, but should be postive.";
551  T_result tmp = raise_domain_error<T_result,T_y>(function,
552  stream.str().c_str(),
553  y[0],
554  Policy());
555  if (result != 0)
556  *result = tmp;
557  return false;
558  }
559  for (size_t n = 1; n < y.size(); n++) {
560  if (!(y[n] > y[n-1])) {
561  std::ostringstream stream;
562  stream << name << " is not a valid ordered vector."
563  << " The element at " << n
564  << " is %1%, but should be greater than the previous element, "
565  << y[n-1];
566  T_result tmp = raise_domain_error<T_result,T_y>(function,
567  stream.str().c_str(),
568  y[n],
569  Policy());
570  if (result != 0)
571  *result = tmp;
572  return false;
573  }
574  }
575  return true;
576  }
577  template <typename T_y, typename T_result>
578  bool check_positive_ordered(const char* function,
579  const Eigen::Matrix<T_y,Eigen::Dynamic,1>& y,
580  const char* name,
581  T_result* result) {
582  return check_positive_ordered(function,y,name,result,default_policy());
583  }
584  template <typename T>
585  bool check_positive_ordered(const char* function,
586  const Eigen::Matrix<T,Eigen::Dynamic,1>& y,
587  const char* name,
588  T* result = 0) {
589  return check_positive_ordered(function,y,name,result,default_policy());
590  }
591 
592 
593 
594  // error_handling functions for Eigen matrix
595 
596  template <typename T_y, typename T_result, class Policy>
597  inline bool check_not_nan(const char* function,
598  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
599  const char* name,
600  T_result* result,
601  const Policy&) {
602  for (int i = 0; i < y.rows(); i++) {
603  for (int j = 0; j < y.cols(); j++) {
604  if (boost::math::isnan(y(i,j))) {
605  std::ostringstream message;
606  message << name << "[" << i << "," << j
607  << "] is %1%, but must not be nan!";
608  T_result tmp
609  = policies::raise_domain_error<T_y>(function,
610  message.str().c_str(),
611  y(i,j), Policy());
612  if (result != 0)
613  *result = tmp;
614  return false;
615  }
616  }
617  }
618  return true;
619  }
620  template <typename T_y, typename T_result>
621  inline bool check_not_nan(const char* function,
622  const Eigen::Matrix<T_y,Eigen::Dynamic,Eigen::Dynamic>& y,
623  const char* name,
624  T_result* result) {
625  return check_not_nan(function,y,name,result,default_policy());
626  }
627  template <typename T>
628  inline bool check_not_nan(const char* function,
629  const Eigen::Matrix<T,Eigen::Dynamic,Eigen::Dynamic>& y,
630  const char* name,
631  T* result = 0) {
632  return check_not_nan(function,y,name,result,default_policy());
633  }
634 
635  }
636 }
637 
638 #endif
     [ Stan Home Page ] © 2011–2013, Stan Development Team.