SLIP  1.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
linear_least_squares.hpp
Go to the documentation of this file.
1 /*
2  * Copyright(c):
3  * Signal Image and Communications (SIC) Department
4  * http://www.sic.sp2mi.univ-poitiers.fr/
5  * - University of Poitiers, France http://www.univ-poitiers.fr
6  * - XLIM Institute UMR CNRS 7252 http://www.xlim.fr/
7  *
8  * and
9  *
10  * D2 Fluid, Thermic and Combustion
11  * - University of Poitiers, France http://www.univ-poitiers.fr
12  * - PPRIME Institute - UPR CNRS 3346 http://www.pprime.fr
13  * - ISAE-ENSMA http://www.ensma.fr
14  *
15  * Contributor(s):
16  * The SLIP team,
17  * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
18  * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>,
19  * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>,
20  * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>,
21  * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>,
22  * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
23  *
24  * Description:
25  * The Simple Library of Image Processing (SLIP) is a new image processing
26  * library. It is written in the C++ language following as much as possible
27  * the ISO/ANSI C++ standard. It is consequently compatible with any system
28  * satisfying the ANSI C++ complience. It works on different Unix , Linux ,
29  * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that
30  * was created by the Signal, Image and Communications (SIC) departement of
31  * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic
32  * and Combustion departement of the P', UPR 3346 CNRS Institute of the
33  * University of Poitiers.
34  *
35  * The SLIP Library source code has been registered to the APP (French Agency
36  * for the Protection of Programs) by the University of Poitiers and CNRS,
37  * under registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.
38 
39  * http://www.sic.sp2mi.univ-poitiers.fr/slip/
40  *
41  * This software is governed by the CeCILL-C license under French law and
42  * abiding by the rules of distribution of free software. You can use,
43  * modify and/ or redistribute the software under the terms of the CeCILL-C
44  * license as circulated by CEA, CNRS and INRIA at the following URL
45  * http://www.cecill.info.
46  * As a counterpart to the access to the source code and rights to copy,
47  * modify and redistribute granted by the license, users are provided only
48  * with a limited warranty and the software's author, the holder of the
49  * economic rights, and the successive licensors have only limited
50  * liability.
51  *
52  * In this respect, the user's attention is drawn to the risks associated
53  * with loading, using, modifying and/or developing or reproducing the
54  * software by the user in light of its specific status of free software,
55  * that may mean that it is complicated to manipulate, and that also
56  * therefore means that it is reserved for developers and experienced
57  * professionals having in-depth computer knowledge. Users are therefore
58  * encouraged to load and test the software's suitability as regards their
59  * requirements in conditions enabling the security of their systems and/or
60  * data to be ensured and, more generally, to use and operate it in the
61  * same conditions as regards security.
62  *
63  * The fact that you are presently reading this means that you have had
64  * knowledge of the CeCILL-C license and that you accept its terms.
65  */
66 
67 
68 
75 #ifndef SLIP_LINEAR_LEAST_SQUARES_HPP
76 #define SLIP_LINEAR_LEAST_SQUARES_HPP
77 
78 #include "Vector.hpp"
79 #include "Matrix.hpp"
80 #include "Array2d.hpp"
81 #include "linear_algebra_svd.hpp"
82 #include "polynomial_algo.hpp"
83 
84 namespace slip
85 {
86 
88  /* @{ */
89 
150  template<typename RandomAccessIterator1,
151  typename RandomAccessIterator2,
152  typename RandomAccessIterator3,
153  typename RandomAccessIterator4>
154  inline
155  double svd_least_square(RandomAccessIterator1 x_first,
156  RandomAccessIterator1 x_last,
157  RandomAccessIterator2 y_first,
158  RandomAccessIterator3 s_first,
159  RandomAccessIterator4 p_first,
160  RandomAccessIterator4 p_last,
161  slip::EvalBasis<typename std::iterator_traits<RandomAccessIterator1>::value_type,RandomAccessIterator2>& basis_fun)
162 {
163  typedef typename std::iterator_traits<RandomAccessIterator2>::difference_type _Distance;
164  _Distance x_size = x_last - x_first;
165  _Distance p_size = p_last - p_first;
166 
167  assert(x_first != x_last);
168  assert(y_first != (y_first + x_size));
169  assert(s_first != (s_first + x_size));
170  assert(p_first != p_last);
171 
172  typedef typename std::iterator_traits<RandomAccessIterator2>::value_type value_type;
173 
174  //Vector B such that bi = yi/si
175  slip::Vector<value_type> B(x_size);
176  //Fitting matrix A
177  slip::Matrix<value_type> A(x_size,p_size);
178 
179  //slip::EvalPowerBasis<value_type,RandomAccessIterator1> basis_fun;
180  //-----------------------------------------------
181  // Computes coefficients of the fitting matrix A
182  //-----------------------------------------------
183  for (_Distance i = 0; i < x_size; ++i)
184  {
185  //computes Fj(x_first[i])
186  basis_fun(x_first[i],(p_size - 1),A.row_begin(i),A.row_end(i));
187  //divides the aij by si
188  double tmp = value_type(1.0) / s_first[i];
189  std::transform(A.row_begin(i),A.row_end(i),
190  A.row_begin(i),std::bind2nd(std::multiplies<value_type>(),tmp));
191  //computes bi = yi/si
192  B[i] = y_first[i] * tmp;
193  }
194  //-----------------------------------------------
195  // svd solve
196  //-----------------------------------------------
197  slip::svd_solve(A,p_first,p_last,B.begin(),B.end());
198 
199  //-----------------------------------------------
200  // Evaluate chi-square
201  //-----------------------------------------------
202  value_type chisq = value_type(0);
203  for (_Distance i = 0; i < x_size; ++i)
204  {
205  value_type err = (B[i] - std::inner_product(p_first,p_last,A.row_begin(i),value_type(0)));
206  chisq += (err * err);
207 
208  }
209  return double(chisq);
210 }
211 
264  template<typename Vector1,
265  typename Vector2,
266  typename Vector3,
267  typename Vector4>
268  inline
269  double svd_least_square(Vector1& X,
270  Vector2& Y,
271  Vector3& S,
272  Vector4& P,
273  slip::EvalBasis<typename Vector1::value_type,typename Vector2::iterator>& basis_fun)
274 {
275  return slip::svd_least_square<typename Vector1::iterator,
276  typename Vector2::iterator,
277  typename Vector3::iterator,
278  typename Vector4::iterator>(X.begin(),X.end(),Y.begin(),S.begin(),P.begin(),P.end(),basis_fun);
279 }
280 
281 
282 
339  template<typename RandomAccessIterator1,
340  typename RandomAccessIterator2,
341  typename RandomAccessIterator3>
342  inline
343  double svd_least_square(RandomAccessIterator1 x_first,
344  RandomAccessIterator1 x_last,
345  RandomAccessIterator2 y_first,
346  RandomAccessIterator3 p_first,
347  RandomAccessIterator3 p_last,
348  slip::EvalBasis<typename std::iterator_traits<RandomAccessIterator1>::value_type,RandomAccessIterator2>& basis_fun)
349 {
350  typedef typename std::iterator_traits<RandomAccessIterator2>::difference_type _Distance;
351  _Distance x_size = x_last - x_first;
352  _Distance p_size = p_last - p_first;
353 
354  assert(x_first != x_last);
355  assert(y_first != (y_first + x_size));
356  assert(p_first != p_last);
357 
358  typedef typename std::iterator_traits<RandomAccessIterator3>::value_type value_type;
359 
360  //Fitting matrix A
361  slip::Matrix<value_type> A(x_size,p_size);
362 
363  //-----------------------------------------------
364  // Computes coefficients of the fitting matrix A
365  //-----------------------------------------------
366  for (_Distance i = 0; i < x_size; ++i)
367  {
368  //computes Fj(x_first[i])
369  basis_fun(x_first[i],(p_size - 1),A.row_begin(i),A.row_end(i));
370  }
371 
372  //-----------------------------------------------
373  // svd solve
374  //-----------------------------------------------
375  slip::svd_solve(A,p_first,p_last,y_first,y_first+x_size);
376 
377  //-----------------------------------------------
378  // Evaluate chi-square
379  //-----------------------------------------------
380  value_type chisq = value_type(0);
381  for (_Distance i = 0; i < x_size; ++i)
382  {
383  value_type err = (y_first[i] - std::inner_product(p_first,p_last,A.row_begin(i),value_type(0)));
384  chisq += (err * err);
385 
386  }
387  return double(chisq);
388 }
389 
390 
439  template<typename Vector1,
440  typename Vector2,
441  typename Vector3>
442  inline
443  double svd_least_square(Vector1& X,
444  Vector2& Y,
445  Vector3& P,
446  slip::EvalBasis<typename Vector1::value_type,typename Vector2::iterator>& basis_fun)
447 {
448  return slip::svd_least_square<typename Vector1::iterator,
449  typename Vector2::iterator,
450  typename Vector3::iterator>(X.begin(),X.end(),Y.begin(),P.begin(),P.end(),basis_fun);
451 }
452 /* @} */
453 
455  /* @{ */
456 
544  template<typename RandomAccessIterator1,
545  typename RandomAccessIterator2,
546  typename RandomAccessIterator3,
547  typename RandomAccessIterator4>
548  inline
549  double svd_least_square_nd(RandomAccessIterator1 x_first,
550  RandomAccessIterator1 x_last,
551  RandomAccessIterator2 y_first,
552  RandomAccessIterator3 s_first,
553  RandomAccessIterator4 p_first,
554  RandomAccessIterator4 p_last,
555  const std::size_t order,
556  slip::EvalBasis<typename std::iterator_traits<RandomAccessIterator1>::value_type,RandomAccessIterator2>& basis_fun)
557 {
558  typedef typename std::iterator_traits<RandomAccessIterator2>::difference_type _Distance;
559  _Distance x_size = x_last - x_first;
560  _Distance p_size = p_last - p_first;
561 
562  assert(x_first != x_last);
563  assert(y_first != (y_first + x_size));
564  assert(s_first != (s_first + x_size));
565  assert(p_first != p_last);
566 
567  typedef typename std::iterator_traits<RandomAccessIterator2>::value_type value_type;
568 
569  //Vector B such that bi = yi/si
570  slip::Vector<value_type> B(x_size);
571  //Fitting matrix A
572  slip::Matrix<value_type> A(x_size,p_size);
573 
574  //-----------------------------------------------
575  // Computes coefficients of the fitting matrix A
576  //-----------------------------------------------
577  for (_Distance i = 0; i < x_size; ++i)
578  {
579  //computes Fj(x_first[i])
580 
581  basis_fun(x_first[i],order,A.row_begin(i),A.row_end(i));
582 
583  //divides the aij by si
584  double tmp = value_type(1.0) / s_first[i];
585  std::transform(A.row_begin(i),A.row_end(i),
586  A.row_begin(i),std::bind2nd(std::multiplies<value_type>(),tmp));
587  //computes bi = yi/si
588  B[i] = y_first[i] * tmp;
589  }
590 
591  //-----------------------------------------------
592  // svd solve
593  //-----------------------------------------------
594  slip::svd_solve(A,p_first,p_last,B.begin(),B.end());
595 
596  //-----------------------------------------------
597  // Evaluate chi-square
598  //-----------------------------------------------
599  value_type chisq = value_type(0);
600  for (_Distance i = 0; i < x_size; ++i)
601  {
602  value_type err = (B[i] - std::inner_product(p_first,p_last,A.row_begin(i),value_type(0)));
603  chisq += (err * err);
604 
605  }
606  return double(chisq);
607 }
608 
609 
693  template<typename RandomAccessIterator1,
694  typename RandomAccessIterator2,
695  typename RandomAccessIterator3>
696  inline
697  double svd_least_square_nd(RandomAccessIterator1 x_first,
698  RandomAccessIterator1 x_last,
699  RandomAccessIterator2 y_first,
700  RandomAccessIterator3 p_first,
701  RandomAccessIterator3 p_last,
702  const std::size_t order,
703  slip::EvalBasis<typename std::iterator_traits<RandomAccessIterator1>::value_type,RandomAccessIterator2>& basis_fun)
704 {
705  typedef typename std::iterator_traits<RandomAccessIterator2>::difference_type _Distance;
706  _Distance x_size = x_last - x_first;
707  _Distance p_size = p_last - p_first;
708 
709  assert(x_first != x_last);
710  assert(y_first != (y_first + x_size));
711  assert(p_first != p_last);
712 
713  typedef typename std::iterator_traits<RandomAccessIterator2>::value_type value_type;
714 
715  //Fitting matrix A
716  slip::Matrix<value_type> A(x_size,p_size);
717 
718  //-----------------------------------------------
719  // Computes coefficients of the fitting matrix A
720  //-----------------------------------------------
721  for (_Distance i = 0; i < x_size; ++i)
722  {
723  //computes Fj(x_first[i])
724 
725  basis_fun(x_first[i],order,A.row_begin(i),A.row_end(i));
726  }
727 
728  //-----------------------------------------------
729  // svd solve
730  //-----------------------------------------------
731  slip::svd_solve(A,p_first,p_last,y_first,y_first+x_size);
732 
733 
734  //-----------------------------------------------
735  // Evaluate chi-square
736  //-----------------------------------------------
737  value_type chisq = value_type(0);
738  for (_Distance i = 0; i < x_size; ++i)
739  {
740  value_type err = (y_first[i] - std::inner_product(p_first,p_last,A.row_begin(i),value_type(0)));
741  chisq += (err * err);
742 
743  }
744  return double(chisq);
745 }
746 
822  template<typename Vector1,
823  typename Vector2,
824  typename Vector3>
825  inline
826  double svd_least_square_nd(Vector1& X,
827  Vector2& Y,
828  Vector3& P,
829  const std::size_t order,
830  slip::EvalBasis<typename Vector1::value_type,typename Vector2::iterator>& basis_fun)
831 {
832  return slip::svd_least_square_nd<typename Vector1::iterator,
833  typename Vector2::iterator,
834  typename Vector3::iterator>(X.begin(),X.end(),Y.begin(),P.begin(),P.end(),order,basis_fun);
835 }
836 /* @} */
837 }//slip::
838 
839 
840 #endif //SLIP_LINEAR_LEAST_SQUARES_HPP
Array2d.hpp
Provides a class to manipulate 2d dynamic and generic arrays.
slip::Vector::end
iterator end()
Returns a read/write iterator that points one past the last element in the Vector. Iteration is done in ordinary element order.
Definition: Vector.hpp:1481
slip::svd_solve
void svd_solve(const Matrix1 &U, RandomAccessIterator1 W_first, RandomAccessIterator1 W_last, const Matrix2 &V, RandomAccessIterator2 B_first, RandomAccessIterator2 B_last, RandomAccessIterator3 X_first, RandomAccessIterator3 X_last)
Solve Linear (USV^H).X = B where U, S and V are provided by a svd decomposition.
Definition: linear_algebra_svd.hpp:950
slip::svd_least_square_nd
double svd_least_square_nd(RandomAccessIterator1 x_first, RandomAccessIterator1 x_last, RandomAccessIterator2 y_first, RandomAccessIterator3 s_first, RandomAccessIterator4 p_first, RandomAccessIterator4 p_last, const std::size_t order, slip::EvalBasis< typename std::iterator_traits< RandomAccessIterator1 >::value_type, RandomAccessIterator2 > &basis_fun)
nd Linear Least Square fitting using SVD.
Definition: linear_least_squares.hpp:549
polynomial_algo.hpp
Provides some polynomial algorithms.
slip::Vector
Numerical Vector class. This container statisfies the RandomAccessContainer concepts of the Standard ...
Definition: Vector.hpp:104
slip::Matrix::row_end
row_iterator row_end(const size_type row)
Returns a read/write iterator that points one past the end element of the row row in the Matrix...
slip::Matrix
Numerical matrix class. This container statisfies the BidirectionnalContainer concepts of the STL...
Definition: GenericMultiComponent2d.hpp:132
slip::inner_product
std::iterator_traits< RandomAccessIterator1 >::value_type inner_product(RandomAccessIterator1 first1, RandomAccessIterator1 last1, RandomAccessIterator2 first2, typename std::iterator_traits< RandomAccessIterator1 >::value_type init=typename std::iterator_traits< RandomAccessIterator1 >::value_type())
Computes the inner_product of two ranges X and Y: .
Definition: linear_algebra.hpp:323
slip::svd_least_square
double svd_least_square(RandomAccessIterator1 x_first, RandomAccessIterator1 x_last, RandomAccessIterator2 y_first, RandomAccessIterator3 s_first, RandomAccessIterator4 p_first, RandomAccessIterator4 p_last, slip::EvalBasis< typename std::iterator_traits< RandomAccessIterator1 >::value_type, RandomAccessIterator2 > &basis_fun)
Linear Least Square fitting using SVD.
Definition: linear_least_squares.hpp:155
Matrix.hpp
Provides a class to manipulate Numerical Matrix.
slip::Matrix::row_begin
row_iterator row_begin(const size_type row)
Returns a read/write iterator that points to the first element of the row row in the Matrix...
linear_algebra_svd.hpp
Provides some Singular Value Decomposition (SVD) algorithms.
Vector.hpp
Provides a class to manipulate numerical vectors.
slip::Vector::begin
iterator begin()
Returns a read/write iterator that points to the first element in the Vector. Iteration is done in or...
Definition: Vector.hpp:1474
Generated on Tue Sep 30 2014 10:05:57 for SLIP by   doxygen 1.8.6