Changeset 804 in Sophya for trunk/SophyaLib/TArray/sopemtx.cc


Ignore:
Timestamp:
Apr 3, 2000, 7:36:01 PM (25 years ago)
Author:
ansari
Message:

Amelioation / debugging de la classe TArray<T> - TVector et TMatrix

heritent maintenant de TArray<T> - Classe RCMatrix rendu prive au fichier
sopemtx.cc - linfit.cc integre a sopemtx.cc

Reza 03/04/2000

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/TArray/sopemtx.cc

    r772 r804  
    1414// -------------------------------------------------------------
    1515////////////////////////////////////////////////////////////////
     16/////////////////////////////////////////////////////////////////////////
     17// Classe de lignes/colonnes de matrices
     18enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
     19template <class T>
     20class TMatrixRC {
     21public:
     22  TMatrixRC();
     23  TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
     24  virtual ~TMatrixRC() {}
     25
     26  // Acces aux rangees et colonnes de matrices
     27  static TMatrixRC<T> Row(TMatrix<T> & m, uint_4 r);
     28  static TMatrixRC<T> Col(TMatrix<T> & m, uint_4 c);
     29  static TMatrixRC<T> Diag(TMatrix<T> & m);
     30
     31  // ---- A virer   $CHECK$ Reza 03/2000
     32  //  int_4 Next();
     33  //  int_4 Prev();
     34  int_4 SetCol(int_4 c);
     35  int_4 SetRow(int_4 r);
     36  int_4 SetDiag();
     37
     38  static uint_4 Step(const TMatrix<T>& m, TRCKind rckind);
     39  static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);
     40
     41  TRCKind Kind() const { return kind; }
     42  uint_4 NElts() const;
     43  T& operator()(uint_4 i);
     44  T  operator()(uint_4 i) const;
     45
     46
     47   TMatrixRC<T>& operator = (const TMatrixRC<T>& rc);
     48
     49// ---- A virer   $CHECK$ Reza 03/2000
     50//   TVector<T> GetVect() const;
     51//   TMatrixRC<T>& operator += (const TMatrixRC<T>& rc);
     52//   TMatrixRC<T>& operator -= (const TMatrixRC<T>& rc);
     53
     54  TMatrixRC<T>& operator *= (T x);
     55  TMatrixRC<T>& operator /= (T x);
     56
     57//   TMatrixRC<T>& operator -= (T x);
     58//   TMatrixRC<T>& operator += (T x);
     59 
     60
     61  TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, uint_4 first=0);
     62  TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, uint_4 first=0);
     63
     64  uint_4 IMaxAbs(uint_4 first=0);
     65
     66  static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
     67
     68  inline static double Abs_Value(uint_1 v) {return (double) v;}
     69  inline static double Abs_Value(uint_2 v) {return (double) v;}
     70  inline static double Abs_Value(int_2 v)  {return (v>0)? (double) v: (double) -v;}
     71  inline static double Abs_Value(int_4 v)  {return (v>0)? (double) v: (double) -v;}
     72  inline static double Abs_Value(int_8 v)  {return (v>0)? (double) v: (double) -v;}
     73  inline static double Abs_Value(uint_4 v) {return (double) v;}
     74  inline static double Abs_Value(uint_8 v) {return (double) v;}
     75  inline static double Abs_Value(r_4 v)    {return (double) fabsf(v);}
     76  inline static double Abs_Value(r_8 v)    {return fabs(v);}
     77  inline static double Abs_Value(complex<float> v)
     78                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
     79  inline static double Abs_Value(complex<double> v)
     80                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
     81
     82protected:
     83  TMatrix<T>* matrix;
     84  T*          data;
     85  int_4       index;
     86  uint_4      step;
     87  TRCKind     kind;
     88};
     89
     90
     91
     92template <class T>
     93inline T operator * (const TMatrixRC<T>& a, const TMatrixRC<T>& b)
     94  {
     95  if ( a.NElts() != b.NElts() )
     96    throw(SzMismatchError("TMatrixRC::operator * size mismatch\n"));
     97  if ( a.Kind() != b.Kind() )
     98    throw(SzMismatchError("TMatrixRC::operator * type mismatch\n"));
     99  T sum = 0;
     100  for(uint_4 i=0; i<a.NElts(); i++) sum += a(i)*b(i);
     101  return sum;
     102  }
     103
     104
     105template <class T>
     106inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
     107  { switch (rckind) { case TmatrixRow  : return m.Step(m.RowsKA());
     108                      case TmatrixCol  : return m.Step(m.ColsKA());
     109                      case TmatrixDiag : return (m.Step(m.RowsKA())+m.Step(m.ColsKA())); }
     110    return 0; }
     111
     112template <class T>
     113inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
     114  { switch (rckind) { case TmatrixRow  : return const_cast<T *>(&(m(index,0)));
     115                      case TmatrixCol  : return const_cast<T *>(&(m(0,index)));
     116                      case TmatrixDiag : return const_cast<T *>(m.Data()); }
     117    return NULL; }
     118
     119template <class T> inline uint_4 TMatrixRC<T>::NElts() const
     120  { if (!matrix) return 0;
     121    switch (kind) { case TmatrixRow  : return matrix->NCols();
     122                    case TmatrixCol  : return matrix->NRows();
     123                    case TmatrixDiag : return matrix->NCols(); }
     124    return 0; }
     125
     126template <class T>
     127inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
     128template <class T>
     129inline T  TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}
     130
     131////////////////////////////////////////////////////////////////
     132// Typedef pour simplifier et compatibilite Peida
     133typedef TMatrixRC<r_8> MatrixRC;
     134
    16135
    17136template <class T> TMatrixRC<T>::TMatrixRC()
    … …  
    52171
    53172
    54 template <class T> int_4 TMatrixRC<T>::Next()
    55 {
    56 if (!matrix || kind==TmatrixDiag) return -1;
    57 index++;
    58 if(kind == TmatrixRow) {
    59   if(index > (int_4)matrix->NRows()) {index = (int_4)matrix->NRows(); return -1;}
    60   data += matrix->NCols();
    61 } else {
    62   if (index > (int_4)matrix->NCols()) {index = (int_4)matrix->NCols(); return -1;}
    63   data++;
    64 }
    65 return index;
    66 }
    67 
    68 template <class T> int_4 TMatrixRC<T>::Prev()
    69 {
    70 if (!matrix || kind == TmatrixDiag) return -1;
    71 index--;
    72 if(index < 0) {index = 0; return -1;}
    73 if(kind == TmatrixRow) data -= matrix->NCols();
    74 else data--;
    75 return index;
    76 }
    77 
     173// ---- A virer   $CHECK$ Reza 03/2000
     174// template <class T> int_4 TMatrixRC<T>::Next()
     175// {
     176// if (!matrix || kind==TmatrixDiag) return -1;
     177// index++;
     178// if(kind == TmatrixRow) {
     179//   if(index > (int_4)matrix->NRows()) {index = (int_4)matrix->NRows(); return -1;}
     180//   data += matrix->NCols();
     181// } else {
     182//   if (index > (int_4)matrix->NCols()) {index = (int_4)matrix->NCols(); return -1;}
     183//   data++;
     184// }
     185// return index;
     186// }
     187
     188// template <class T> int_4 TMatrixRC<T>::Prev()
     189// {
     190// if (!matrix || kind == TmatrixDiag) return -1;
     191// index--;
     192// if(index < 0) {index = 0; return -1;}
     193// if(kind == TmatrixRow) data -= matrix->NCols();
     194// else data--;
     195// return index;
     196// }
     197
     198 
    78199template <class T> int_4 TMatrixRC<T>::SetCol(int_4 c)
    79200{
    … …  
    121242}
    122243
    123 template <class T> TVector<T> TMatrixRC<T>::GetVect() const
    124 {
    125 TVector<T> v(NElts());
    126 for (uint_4 i=0; i<NElts(); i++) v(i) = (*this)(i);
    127 return v;
    128 }
    129 
    130 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (const TMatrixRC<T>& rc)
    131 {
    132 if ( NElts() != rc.NElts() )
    133   throw(SzMismatchError("TMatrixRC::operator+= size mismatch\n"));
    134 if ( kind != rc.kind )
    135   throw(SzMismatchError("TMatrixRC::operator+= type mismatch\n"));
    136 for (uint_4 i=0; i<NElts(); i++) (*this)(i) += rc(i);
    137 return *this;
    138 }
    139 
    140 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (const TMatrixRC<T>& rc)
    141 {
    142 if( NElts() != rc.NElts() )
    143   throw(SzMismatchError("TMatrixRC::operator-= size mismatch\n"));
    144 if( kind != rc.kind )
    145   throw(SzMismatchError("TMatrixRC::operator-= type mismatch\n"));
    146 for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= rc(i);
    147 return *this;
    148 }
     244// ---- A virer   $CHECK$ Reza 03/2000
     245// template <class T> TVector<T> TMatrixRC<T>::GetVect() const
     246// {
     247// TVector<T> v(NElts());
     248// for (uint_4 i=0; i<NElts(); i++) v(i) = (*this)(i);
     249// return v;
     250// }
     251
     252// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (const TMatrixRC<T>& rc)
     253// {
     254// if ( NElts() != rc.NElts() )
     255//   throw(SzMismatchError("TMatrixRC::operator+= size mismatch\n"));
     256// if ( kind != rc.kind )
     257//   throw(SzMismatchError("TMatrixRC::operator+= type mismatch\n"));
     258// for (uint_4 i=0; i<NElts(); i++) (*this)(i) += rc(i);
     259// return *this;
     260// }
     261
     262// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (const TMatrixRC<T>& rc)
     263// {
     264// if( NElts() != rc.NElts() )
     265//   throw(SzMismatchError("TMatrixRC::operator-= size mismatch\n"));
     266// if( kind != rc.kind )
     267//   throw(SzMismatchError("TMatrixRC::operator-= type mismatch\n"));
     268// for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= rc(i);
     269// return *this;
     270// }
    149271
    150272
    … …  
    162284
    163285
    164 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (T x)
    165 {
    166 for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= x;
    167 return *this;
    168 }
    169 
    170 template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (T x)
    171 {
    172 for(uint_4 i=0; i<NElts(); i++) (*this)(i) += x;
    173 return *this;
    174 }
     286// ---- A virer   $CHECK$ Reza 03/2000
     287// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (T x)
     288// {
     289// for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= x;
     290// return *this;
     291// }
     292
     293// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (T x)
     294// {
     295// for(uint_4 i=0; i<NElts(); i++) (*this)(i) += x;
     296// return *this;
     297// }
    175298
    176299template <class T>
    … …  
    331454
    332455
     456LinFitter::LinFitter()
     457{
     458}
     459
     460LinFitter::~LinFitter()
     461{
     462}
     463
     464double LinFitter::LinFit(const Vector& x, const Vector& y, int nf,
     465                         double (*f)(int, double), Vector& c)
     466{
     467  int n = x.NElts();
     468  if (n != y.NElts()) THROW(sizeMismatchErr);
     469 
     470  Matrix fx(nf, n);
     471  for (int i=0; i<nf; i++)
     472    for (int j=0; j<n; j++)
     473      fx(i,j) = f(i,x(j));
     474
     475  return LinFit(fx,y,c);
     476}
     477
     478
     479
     480double LinFitter::LinFit(const Matrix& fx, const Vector& y, Vector& c)
     481{
     482  int n = y.NElts();
     483  if ( n != fx.NCol()) THROW(sizeMismatchErr);
     484
     485  int nf = fx.NRows();
     486
     487  Matrix a(nf,nf);
     488
     489  for (int j=0; j<nf; j++)
     490    for (int k=j; k<nf; k++)
     491      a(j,k) = a(k,j) = TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), j)
     492
     493        * TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), k); /* $CHECK$ Reza 10/3/2000 */
     494   
     495  c = fx * y;
     496
     497  if (SimpleMatrixOperation<r_8>::GausPiv(a,c) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
     498
     499  double xi2 = 0;
     500
     501  for (int k=0; k<n; k++) {
     502    double x = 0;
     503    for (int i=0; i<nf; i++)
     504      x += c(i) * fx(i,k);
     505    x -= y(k);
     506    xi2 += x*x;
     507  }
     508  return xi2;
     509}
     510
     511
     512
     513double LinFitter::LinFit(const Vector& x, const Vector& y, const Vector& errY2, int nf,
     514                         double (*f)(int, double), Vector& c, Vector& errC)
     515{
     516  int n = x.NElts();
     517  if (n != y.NElts()) THROW(sizeMismatchErr);
     518 
     519  Matrix fx(nf, n);
     520  for (int i=0; i<nf; i++)
     521    for (int j=0; j<n; j++)
     522      fx(i,j) = f(i,x(j));
     523
     524  return LinFit(fx,y,errY2,c,errC);
     525}
     526
     527
     528double LinFitter::LinFit(const Matrix& fx, const Vector& y, const Vector& errY2,
     529           Vector& c, Vector& errC)
     530{
     531  int n = y.NElts();
     532  if ( n != errY2.NElts()) THROW(sizeMismatchErr);
     533  if ( n != fx.NCol()) THROW(sizeMismatchErr);
     534
     535  int nf = fx.NRows();
     536
     537  Matrix a(nf,nf);
     538
     539  c.Realloc(nf);
     540  errC.Realloc(nf);
     541
     542  for (int j=0; j<nf; j++)
     543    for (int k=j; k<nf; k++) {
     544      double x=0;
     545      for (int l=0; l<n; l++)
     546        x += fx(j,l) * fx(k,l) / errY2(l);        // Matrice a inverser
     547      a(j,k) = a(k,j) = x;
     548    }
     549 
     550  Matrix d(nf,nf+1);
     551  for (int k=0; k<nf; k++) {
     552    double x=0;
     553    for (int l=0; l<n; l++)
     554      x += fx(k,l) * y(l) / errY2(l);             // Second membre 1ere colonne
     555    d(k,0) = x;
     556    for (int m=1; m<=nf; m++)
     557      d(k,m) = (k==m) ? 1 : 0;                // Reste second membre = Id.
     558  }
     559
     560  if (SimpleMatrixOperation<r_8>::GausPiv(a,d) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
     561
     562  for (int l=0; l<nf; l++) {
     563    c(l) = d(l,0);                            // Parametres = 1ere colonne
     564    errC(l) = d(l,l+1);                        // Erreurs = diag inverse.
     565  }
     566
     567  double xi2 = 0;
     568
     569  for (int jj=0; jj<n; jj++) {
     570    double x = 0;
     571    for (int ii=0; ii<nf; ii++)
     572      x += c(ii) * fx(ii,jj);
     573    x -= y(jj);
     574    xi2 += x*x/errY2(jj);
     575  }
     576  return xi2;
     577}
     578
     579
     580
    333581
    334582///////////////////////////////////////////////////////////////
    335583#ifdef __CXX_PRAGMA_TEMPLATES__
    336584// Instances gestion lignes/colonnes
    337 #pragma define_template TMatrixRC<uint_1>
    338 #pragma define_template TMatrixRC<uint_2>
    339 #pragma define_template TMatrixRC<int_2>
    340585#pragma define_template TMatrixRC<int_4>
    341 #pragma define_template TMatrixRC<int_8>
    342 #pragma define_template TMatrixRC<uint_4>
    343 #pragma define_template TMatrixRC<uint_8>
    344586#pragma define_template TMatrixRC<r_4>
    345587#pragma define_template TMatrixRC<r_8>
    346 #pragma define_template TMatrixRC< complex<r_4> >
    347 #pragma define_template TMatrixRC< complex<r_8> >
    348 // #pragma define_template SimpleMatrixOperation<r_4>
     588#pragma define_template SimpleMatrixOperation<int_4>
     589#pragma define_template SimpleMatrixOperation<r_4>
    349590#pragma define_template SimpleMatrixOperation<r_8>
    350591#endif
    … …  
    352593#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
    353594// Instances gestion lignes/colonnes
    354 template class TMatrixRC<uint_1>;
    355 template class TMatrixRC<uint_2>;
    356 template class TMatrixRC<int_2>;
    357595template class TMatrixRC<int_4>;
    358 template class TMatrixRC<int_8>;
    359 template class TMatrixRC<uint_4>;
    360 template class TMatrixRC<uint_8>;
    361596template class TMatrixRC<r_4>;
    362597template class TMatrixRC<r_8>;
    363 template class TMatrixRC< complex<r_4> >;
    364 template class TMatrixRC< complex<r_8> >;
     598template class SimpleMatrixOperation<int_4>;
    365599template class SimpleMatrixOperation<r_4>;
    366600template class SimpleMatrixOperation<r_8>;
Note: See TracChangeset for help on using the changeset viewer.