Changeset 939 in Sophya for trunk/SophyaLib/TArray


Ignore:
Timestamp:
Apr 14, 2000, 6:15:36 PM (25 years ago)
Author:
ansari
Message:

LinFitter -> passe template cmv 14/4/00

Location:
trunk/SophyaLib/TArray
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/TArray/sopemtx.cc

    r935 r939  
    1010
    1111////////////////////////////////////////////////////////////////
    12 // -------------------------------------------------------------
    13 //   La classe de gestion des lignes et colonnes d'une matrice
    14 // -------------------------------------------------------------
     12// ---------------------------------------------------------- //
     13//  La classe de gestion des lignes et colonnes d'une matrice //
     14// ---------------------------------------------------------- //
    1515////////////////////////////////////////////////////////////////
     16
     17/*!
     18  \class SOPHYA::TMatrixRC
     19  \ingroup TArray
     20  Class for representing rows, columns and diagonal of a matrix.
     21  \sa TMatrix TArray
     22  */
    1623
    1724//! Class of line, column or diagonal of a TMatrix
    … …  
    109116};
    110117
    111 
    112118//! Scalar product of two TMatrixRC
    113119/*!
    … …  
    423429
    424430////////////////////////////////////////////////////////////////
    425 // -------------------------------------------------------------
    426 //   La classe de calcul simple sur les TMatrix
    427 // -------------------------------------------------------------
     431// ---------------------------------------------------------- //
     432//      La classe de calcul simple sur les TMatrix            //
     433// ---------------------------------------------------------- //
    428434////////////////////////////////////////////////////////////////
    429435
    … …  
    542548
    543549////////////////////////////////////////////////////////////////
    544 // -------------------------------------------------------------
    545 //   La classe fit lineaire
    546 // -------------------------------------------------------------
     550// ---------------------------------------------------------- //
     551//               La classe fit lineaire                       //
     552// ---------------------------------------------------------- //
    547553////////////////////////////////////////////////////////////////
    548554
    549 LinFitter::LinFitter()
    550 {
    551 }
    552 
    553 LinFitter::~LinFitter()
    554 {
    555 }
    556 
    557 double LinFitter::LinFit(const Vector& x, const Vector& y, int nf,
    558                          double (*f)(int, double), Vector& c)
    559 {
    560   int n = x.NElts();
    561   if (n != y.NElts()) THROW(sizeMismatchErr);
    562  
    563   Matrix fx(nf, n);
    564   for (int i=0; i<nf; i++)
    565     for (int j=0; j<n; j++)
    566       fx(i,j) = f(i,x(j));
    567 
    568   return LinFit(fx,y,c);
    569 }
    570 
    571 
    572 
    573 double LinFitter::LinFit(const Matrix& fx, const Vector& y, Vector& c)
    574 {
    575   int n = y.NElts();
    576   if ( n != fx.NCol()) THROW(sizeMismatchErr);
    577 
    578   int nf = fx.NRows();
    579 
    580   Matrix a(nf,nf);
    581 
    582   for (int j=0; j<nf; j++)
    583     for (int k=j; k<nf; k++)
    584       a(j,k) = a(k,j) = TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), j)
    585 
    586         * TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), k); /* $CHECK$ Reza 10/3/2000 */
    587    
    588   c = fx * y;
    589 
    590   if (SimpleMatrixOperation<r_8>::GausPiv(a,c) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
    591 
    592   double xi2 = 0;
    593 
    594   for (int k=0; k<n; k++) {
    595     double x = 0;
    596     for (int i=0; i<nf; i++)
    597       x += c(i) * fx(i,k);
    598     x -= y(k);
    599     xi2 += x*x;
     555//! Creator
     556template <class T>
     557LinFitter<T>::LinFitter()
     558{
     559}
     560
     561//! Destructor
     562template <class T>
     563LinFitter<T>::~LinFitter()
     564{
     565}
     566
     567// fit lineaire des y(k) en tant que somme de c(i)f(i,x(k)), i=0..nf-1;
     568//! Linear fitting
     569/*!
     570  Linear fit of y(k) as the sum of \f$ c(i)f(i,x(k)), i=0..nf-1 \f$
     571  \param x : vector of X values
     572  \param y : vector of datas
     573  \param nf: number of functions
     574  \param f : f(i,x(k)), i=0..nf-1
     575  \return c : vector of solutions
     576  \return return chisquare
     577 */
     578template <class T>
     579r_8 LinFitter<T>::LinFit(const TVector<T>& x, const TVector<T>& y,
     580                       uint_4 nf, T (*f)(uint_4,T), TVector<T>& c)
     581{
     582uint_4 n = x.NElts();
     583if (n != y.NElts()) THROW(sizeMismatchErr);
     584 
     585TMatrix<T> fx(nf,n);
     586
     587for(uint_4 i=0; i<nf; i++)
     588  for(uint_4 j=0; j<n; j++) fx(i,j) = f(i,x(j));
     589
     590return LinFit(fx,y,c);
     591}
     592
     593// fit lineaire des y(k) en tant que somme de c(i)f(i,x(k)), i=0..nf-1,
     594// la matrice fx contient les valeurs des f: fx(i,j) = f(i, x(j)).
     595//! Linear fitting
     596/*!
     597  Linear fit of y(k) as the sum of \f$ c(i)f(i,x(k)), i=0..nf-1 \f$.
     598  \param fx : matrix which contains \f$ fx(i,j) = f(i, x(j)) \f$.
     599  \param y : vector of datas
     600  \return c : vector of solutions
     601  \return return chisquare
     602 */
     603template <class T>
     604r_8 LinFitter<T>::LinFit(const TMatrix<T>& fx, const TVector<T>& y, TVector<T>& c)
     605{
     606uint_4 n = y.NElts();
     607if (n != fx.NCol()) THROW(sizeMismatchErr);
     608
     609uint_4 nf = fx.NRows();
     610
     611TMatrix<T> a(nf,nf);
     612
     613for(uint_4 j=0; j<nf; j++)
     614  for(uint_4 k=j; k<nf; k++)
     615     a(j,k) = a(k,j) = TMatrixRC<T>::Row(const_cast<TMatrix<T> &>(fx), j)
     616                     * TMatrixRC<T>::Row(const_cast<TMatrix<T> &>(fx), k);
     617                                /* $CHECK$ Reza 10/3/2000 */
     618
     619c = fx * y;
     620
     621if(SimpleMatrixOperation<T>::GausPiv(a,c)==(T) 0) THROW(singMatxErr);
     622                   /* $CHECK$ Reza 10/3/2000 */
     623
     624r_8 xi2 = 0., ax;
     625T x;
     626for(uint_4 k=0; k<n; k++) {
     627  x = (T) 0;
     628  for(uint_4 i=0; i<nf; i++) x += c(i) * fx(i,k);
     629  x -= y(k);
     630  ax = TMatrixRC<T>::Abs_Value(x);
     631  xi2 += ax*ax;
     632}
     633return xi2;
     634}
     635
     636// fit lineaire des y(k) en tant que somme de c(i)f(i,x(k)), i=0..nf-1,
     637// errY2 contient les carres des erreurs sur les Y.
     638// au retour, errC contient les erreurs sur les coefs.
     639//! Linear fitting with errors
     640/*!
     641  Linear fit with errors of y(k) as the sum of \f$ c(i)f(i,x(k)), i=0..nf-1 \f$.
     642  \param x : vector of X values
     643  \param y : vector of datas
     644  \param errY2 : vector of errors square on Y
     645  \param nf: number of functions
     646  \param f : f(i,x(k)), i=0..nf-1
     647  \return c : vector of solutions
     648  \return errC : vector of errors on solutions C
     649  \return return chisquare
     650 */
     651template <class T>
     652r_8 LinFitter<T>::LinFit(const TVector<T>& x, const TVector<T>& y,
     653             const TVector<T>& errY2, uint_4 nf, T (*f)(uint_4,T),
     654                                  TVector<T>& c, TVector<T>& errC)
     655{
     656uint_4 n = x.NElts();
     657if (n != y.NElts()) THROW(sizeMismatchErr);
     658
     659TMatrix<T> fx(nf, n);
     660for(uint_4 i=0; i<nf; i++)
     661  for(uint_4 j=0; j<n; j++)
     662    fx(i,j) = f(i,x(j));
     663
     664return LinFit(fx,y,errY2,c,errC);
     665}
     666
     667// fit lineaire des y(k) en tant que somme de c(i)f(i,x(k)), i=0..nf-1,
     668// la matrice fx contient les valeurs des f:
     669// fx(i,j) = f(i, x(j)).
     670// errY2 contient les carres des erreurs sur les Y.
     671// au retour, errC contient les erreurs sur les coefs.
     672//! Linear fitting with errors
     673/*!
     674  \param fx : matrix which contains \f$ fx(i,j) = f(i, x(j)) \f$.
     675  \param y : vector of datas
     676  \param errY2 : vector of errors square on Y
     677  \return c : vector of solutions
     678  \return errC : vector of errors on solutions on C
     679  \return return chisquare
     680 */
     681template <class T>
     682r_8 LinFitter<T>::LinFit(const TMatrix<T>& fx, const TVector<T>& y,
     683           const TVector<T>& errY2,TVector<T>& c, TVector<T>& errC)
     684{
     685uint_4 n = y.NElts();
     686if( n != errY2.NElts()) THROW(sizeMismatchErr);
     687if( n != fx.NCol()) THROW(sizeMismatchErr);
     688
     689uint_4 nf = fx.NRows();
     690
     691TMatrix<T> a(nf,nf);
     692
     693c.Realloc(nf);
     694errC.Realloc(nf);
     695
     696for(uint_4 j=0; j<nf; j++)
     697  for(uint_4 k=j; k<nf; k++) {
     698    T x=0;
     699    // Matrice a inverser
     700    for(uint_4 l=0; l<n; l++) x += fx(j,l)*fx(k,l)/errY2(l);
     701    a(j,k) = a(k,j) = x;
     702  }
     703 
     704TMatrix<T> d(nf,nf+1);
     705for(uint_4 k=0; k<nf; k++) {
     706  T x = (T) 0;
     707  // Second membre 1ere colonne
     708  for(uint_4 l=0; l<n; l++) x += fx(k,l)*y(l)/errY2(l);             
     709  d(k,0) = x;
     710  // Reste second membre = Id.
     711  for(uint_4 m=1; m<=nf; m++) d(k,m) = (k==m)?1:0;
     712}
     713
     714if(SimpleMatrixOperation<T>::GausPiv(a,d)==(T) 0) THROW(singMatxErr);
     715                                            /* $CHECK$ Reza 10/3/2000 */
     716
     717for(uint_4 l=0; l<nf; l++) {
     718  c(l) = d(l,0);        // Parametres = 1ere colonne
     719  errC(l) = d(l,l+1);   // Erreurs = diag inverse.
     720}
     721
     722r_8 xi2 = 0., ax;
     723T x;
     724for(uint_4 jj=0; jj<n; jj++) {
     725  x = (T) 0;
     726  for(uint_4 ii=0; ii<nf; ii++) x += c(ii) * fx(ii,jj);
     727    x -= y(jj);
     728    ax = TMatrixRC<T>::Abs_Value(x);
     729    xi2 += ax*ax/TMatrixRC<T>::Abs_Value(errY2(jj));
    600730  }
    601731  return xi2;
    602732}
    603733
    604 
    605 
    606 double LinFitter::LinFit(const Vector& x, const Vector& y, const Vector& errY2, int nf,
    607                          double (*f)(int, double), Vector& c, Vector& errC)
    608 {
    609   int n = x.NElts();
    610   if (n != y.NElts()) THROW(sizeMismatchErr);
    611  
    612   Matrix fx(nf, n);
    613   for (int i=0; i<nf; i++)
    614     for (int j=0; j<n; j++)
    615       fx(i,j) = f(i,x(j));
    616 
    617   return LinFit(fx,y,errY2,c,errC);
    618 }
    619 
    620 
    621 double LinFitter::LinFit(const Matrix& fx, const Vector& y, const Vector& errY2,
    622            Vector& c, Vector& errC)
    623 {
    624   int n = y.NElts();
    625   if ( n != errY2.NElts()) THROW(sizeMismatchErr);
    626   if ( n != fx.NCol()) THROW(sizeMismatchErr);
    627 
    628   int nf = fx.NRows();
    629 
    630   Matrix a(nf,nf);
    631 
    632   c.Realloc(nf);
    633   errC.Realloc(nf);
    634 
    635   for (int j=0; j<nf; j++)
    636     for (int k=j; k<nf; k++) {
    637       double x=0;
    638       for (int l=0; l<n; l++)
    639         x += fx(j,l) * fx(k,l) / errY2(l);        // Matrice a inverser
    640       a(j,k) = a(k,j) = x;
    641     }
    642  
    643   Matrix d(nf,nf+1);
    644   for (int k=0; k<nf; k++) {
    645     double x=0;
    646     for (int l=0; l<n; l++)
    647       x += fx(k,l) * y(l) / errY2(l);             // Second membre 1ere colonne
    648     d(k,0) = x;
    649     for (int m=1; m<=nf; m++)
    650       d(k,m) = (k==m) ? 1 : 0;                // Reste second membre = Id.
    651   }
    652 
    653   if (SimpleMatrixOperation<r_8>::GausPiv(a,d) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
    654 
    655   for (int l=0; l<nf; l++) {
    656     c(l) = d(l,0);                            // Parametres = 1ere colonne
    657     errC(l) = d(l,l+1);                        // Erreurs = diag inverse.
    658   }
    659 
    660   double xi2 = 0;
    661 
    662   for (int jj=0; jj<n; jj++) {
    663     double x = 0;
    664     for (int ii=0; ii<nf; ii++)
    665       x += c(ii) * fx(ii,jj);
    666     x -= y(jj);
    667     xi2 += x*x/errY2(jj);
    668   }
    669   return xi2;
    670 }
    671 
    672 
     734///////////////////////////////////////////////////////////////////////////
     735///////////////////////////////////////////////////////////////////////////
     736///////////////////////////////////////////////////////////////////////////
     737///////////////////////////////////////////////////////////////////////////
    673738void _ZZ_TestTMatrixRC_YY_(TMatrix<r_8> & m)
    674739{
    … …  
    712777#pragma define_template SimpleMatrixOperation< complex<r_4> >
    713778#pragma define_template SimpleMatrixOperation< complex<r_8> >
     779#pragma define_template LinFitter<r_4>
     780#pragma define_template LinFitter<r_8>
     781#pragma define_template LinFitter< complex<r_4> >
     782#pragma define_template LinFitter< complex<r_8> >
    714783#endif
    715784
    … …  
    726795template class SimpleMatrixOperation< complex<r_4> >;
    727796template class SimpleMatrixOperation< complex<r_8> >;
     797template class LinFitter<r_4>;
     798template class LinFitter<r_8>;
     799template class LinFitter< complex<r_4> >;
     800template class LinFitter< complex<r_8> >;
    728801#endif
    729802

  • trunk/SophyaLib/TArray/sopemtx.h

    r935 r939  
    1010////////////////////////////////////////////////////////////////
    1111//------------------------------------------------------------//
    12 //                     Classe TMatrixRC                       //
     12//      La classe de calcul simple sur les TMatrix            //
    1313//------------------------------------------------------------//
    1414////////////////////////////////////////////////////////////////
    … …  
    174174
    175175//!  Class for linear fitting
     176template <class T>
    176177class LinFitter {
     178
    177179public :
    178                 LinFitter();
    179   virtual       ~LinFitter();
    180 
    181   double        LinFit(const Vector& x, const Vector& y, int nf,
    182                        double (*f)(int, double), Vector& c);
    183 // fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1;
    184 
    185   double        LinFit(const Matrix& fx, const Vector& y, Vector& c);
    186 // fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
    187 // la matrice fx contient les valeurs des f:
    188 // fx(i,j) = f(i, x(j)).
     180
     181         LinFitter();
     182virtual ~LinFitter();
     183
     184//! Linear fitting
     185r_8 LinFit(const TVector<T>& x, const TVector<T>& y,
     186           uint_4 nf, T (*f)(uint_4,T), TVector<T>& c);
     187
     188//! Linear fitting
     189r_8 LinFit(const TMatrix<T>& fx, const TVector<T>& y, TVector<T>& c);
    189190                     
    190   double        LinFit(const Vector& x, const Vector& y, const Vector& errY2, int nf,
    191                        double (*f)(int, double), Vector& c, Vector& errC);
    192 // fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
    193 // errY2 contient les carres des erreurs sur les Y.
    194 // au retour, errC contient les erreurs sur les coefs.
    195 
    196   double        LinFit(const Matrix& fx, const Vector& y, const Vector& errY2,
    197                        Vector& c, Vector& errC);
    198 // fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
    199 // la matrice fx contient les valeurs des f:
    200 // fx(i,j) = f(i, x(j)).
    201 // errY2 contient les carres des erreurs sur les Y.
    202 // au retour, errC contient les erreurs sur les coefs.
     191//! Linear fitting with errors
     192r_8 LinFit(const TVector<T>& x, const TVector<T>& y, const TVector<T>& errY2,
     193           uint_4 nf,T (*f)(uint_4,T), TVector<T>& c, TVector<T>& errC);
     194
     195//! Linear fitting with errors
     196r_8 LinFit(const TMatrix<T>& fx, const TVector<T>& y,
     197           const TVector<T>& errY2, TVector<T>& c, TVector<T>& errC);
    203198};
    204199
    205 
    206200} // Fin du namespace
    207201
Note: See TracChangeset for help on using the changeset viewer.