Reference Manual

Go to the documentation of this file.
00001 //   COCOA class implementation file
00002 //Id:  MatrixMeschach.C
00003 //CAT: Model
00004 //
00005 //   History: v1.0 
00006 //   Pedro Arce
00007 #include <iomanip>
00008 #include <cmath>  
00009 
00010 
00011 #include "Alignment/CocoaUtilities/interface/ALIUtils.h"
00012 #include "Alignment/CocoaFit/interface/MatrixMeschach.h"
00013 
00014 
00015 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00016 MatrixMeschach::MatrixMeschach()
00017 {
00018   //-  std::cout << "creating matrix0 " << std::endl;
00019 }
00020 
00021 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00022 MatrixMeschach::~MatrixMeschach()
00023 {
00024   //-  std::cout << "deleting matrix " << std::endl;
00025   M_FREE(_Mat);
00026 }
00027 
00028 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00029 MatrixMeschach::MatrixMeschach( ALIint NoLin, ALIint NoCol )
00030 {
00031   //-  std::cout << "creating matrix " << std::endl;
00032   _NoLines = NoLin; 
00033   _NoColumns = NoCol;
00034   _Mat = m_get( NoLin, NoCol );
00035   //_data = new ALIdouble[NoCol][NoLin];
00036 }
00037 
00038 
00039 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00040 MatrixMeschach::MatrixMeschach( const MatrixMeschach& mat )
00041 {
00042   //-  std::cout << "creating matrixc " << std::endl;
00043   _NoLines = mat._Mat->m; 
00044   _NoColumns = mat._Mat->n;
00045   _Mat = m_get( mat.NoLines(), mat.NoColumns() );
00046   // std::cout <<  "copy const" << mat._Mat << _Mat << NoLines() << NoColumns() << Mat()->m << Mat()->n <<std::endl;
00047   copy( mat );
00048 }
00049  
00050 
00051 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00052 void MatrixMeschach::copy( const MatrixMeschach& mat) 
00053 {
00054   //  if( ALIUtils::debug >= 5) std::cout <<  "copy matrix" << mat._Mat << " " << _Mat << " L " << mat.NoLines() << " C " << mat.NoColumns() << " l " <<  mat.Mat()->m << " c " << mat.Mat()->n <<std::endl;
00055 
00056   for( ALIint lin=0; lin < _NoLines; lin++ ) {
00057     for( ALIint col=0;  col < _NoColumns; col++ ) {
00058       _Mat->me[lin][col] = mat.MatNonConst()->me[lin][col];
00059     } 
00060   }
00061   //   m_copy( mat._Mat, _Mat );
00062 } 
00063 
00064 
00065 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00066 MatrixMeschach& MatrixMeschach::operator=( const MatrixMeschach& mat )
00067 {
00068   if ( this == 0 ) {
00069       std::cerr << "EXITING: trying to use 'operator=' with an MatrixMeschach for which memory space is not reserved!!!!" << std::endl;
00070       std::exception();
00071   }
00072          
00073   if ( mat._Mat != _Mat ) {
00074       _NoLines = mat._Mat->m;
00075       _NoColumns = mat._Mat->n; 
00076       M_FREE( _Mat );
00077       _Mat = m_get( mat.NoLines(), mat.NoColumns() );
00078       copy( mat );
00079   }
00080    if(ALIUtils::debug >= 9) std::cout <<  "operator=" << mat._Mat << _Mat << NoLines() << NoColumns() << Mat()->m << Mat()->n <<std::endl;
00081   return *this;
00082 }
00083 
00084 
00085 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00086 void MatrixMeschach::operator*=( const MatrixMeschach& mat )
00087 {
00088   //  std::cout << " multiply matrices " << std::endl;
00089 
00090   if( _NoColumns != mat.NoLines() ){
00091     std::cerr << " !! Trying two multiply two matrices when the number of columns of first one is not equal to number of files of second one " << std::endl;
00092     std::cout << " multiplying matrix " << _NoLines  << " x " << _NoColumns << " and " << mat.NoLines() << " x " << mat.NoColumns() << " results in " << _NoLines << " x " << mat.NoColumns() << std::endl;
00093   }
00094   _NoColumns = mat.NoColumns();
00095 
00096   MAT* tempmat = m_get( _NoColumns, _NoLines );
00097   m_transp( _Mat, tempmat); 
00098   //  M_FREE( _Mat );
00099   _Mat = m_get( _NoLines, mat.NoColumns() );
00100   mtrm_mlt( tempmat, mat._Mat, _Mat);
00101 
00102   //  _NoColumns = mat.NoColumns(); 
00103   //  M_FREE(tempmat);
00104 }
00105 
00106 
00107 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00108 void MatrixMeschach::operator+=( const MatrixMeschach& mat )
00109 {
00110   if(_NoLines != mat._NoLines  || _NoColumns != mat._NoColumns ) {
00111     std::cerr << "!!!! cannot sum two matrices with different size" << std::endl
00112          << "MATRIX 1:" << _NoLines << " X " << _NoColumns << std::endl
00113          << "MATRIX 2:" << mat._NoLines << " X " << mat._NoColumns << std::endl;
00114   }
00115   MAT* tempmat = m_get( _NoColumns, _NoLines );
00116   m_copy( _Mat, tempmat); 
00117   M_FREE( _Mat );
00118   _Mat = m_get( _NoLines, _NoColumns );
00119   m_add( tempmat, mat._Mat, _Mat);
00120   M_FREE(tempmat);
00121 
00122 }
00123 
00124 
00125 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00126 void MatrixMeschach::operator*=( const ALIdouble num )
00127 {
00128   for (ALIuint ii=0; ii<_Mat->m; ii++) {
00129       for (ALIuint jj=0; jj<_Mat->n; jj++) {
00130           _Mat->me[ii][jj] *= num;
00131       }
00132   }
00133 
00134 }
00135 
00136 
00137 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00138 MatrixMeschach operator*( const MatrixMeschach& mat1, const MatrixMeschach& mat2 )
00139 {
00140   MatrixMeschach mat1copy = mat1;
00141   if( mat1copy.NoColumns() != mat2.NoLines() ){
00142     std::cerr << " !! Trying two multiply two matrices when the number of columns of first one is not equal to number of files of second one " << std::endl;
00143     std::cout << " multiplying matrix " << mat1copy.NoLines() << " x " << mat1copy.NoColumns() << " and " << mat2.NoLines() << " x " << mat2.NoColumns() << " results in " << mat1copy.NoLines() << " x " << mat2.NoColumns() << std::endl;
00144   }
00145   mat1copy.setNoColumns( mat2.NoColumns() );
00146 
00147   MAT* tempmat = m_get( mat1copy.NoColumns(), mat1copy.NoLines() );
00148   m_transp( mat1copy.MatNonConst(), tempmat);
00149  
00150   //M_FREE( _Mat );
00151   mat1copy.setMat( m_get( mat1copy.NoLines(), mat2.NoColumns() ) );
00152   mtrm_mlt( tempmat, mat2.MatNonConst(), mat1copy.MatNonConst());
00153 
00154   free(tempmat);
00155 
00156   return mat1copy;
00157 
00158 
00159   MatrixMeschach* matout = new MatrixMeschach( mat1copy );
00160 
00161   /*  if(ALIUtils::debug >= 9) std::cout << "add" << mat1copy.NoLines() << mat1copy.NoColumns()
00162        << mat2.NoLines() << mat2.NoColumns()  
00163        << matout.NoLines() << matout.NoColumns() << std::endl;
00164   if(ALIUtils::debug >= 9) std::cout << "addM" << mat1copy.Mat()->m << mat1copy.Mat()->n
00165                                      << mat2.Mat()->m << mat2.Mat()->n
00166        << matout.Mat()->m << matout.Mat()->n << std::endl;
00167   */
00168   //  return MatrixMeschach( matout );
00169   return ( *matout );
00170 
00171 }
00172  
00173 
00174 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00175 MatrixMeschach operator+( const MatrixMeschach& mat1, const MatrixMeschach& mat2 )
00176 {
00177   MatrixMeschach matout( mat1 );
00178   matout += mat2;
00179    if(ALIUtils::debug >= 9) std::cout << "add" << mat1.NoLines() << mat1.NoColumns()
00180        << mat2.NoLines() << mat2.NoColumns()  
00181        << matout.NoLines() << matout.NoColumns() << std::endl;
00182    if(ALIUtils::debug >= 9) std::cout << "addM" << mat1.Mat()->m << mat1.Mat()->n
00183        << mat2.Mat()->m << mat2.Mat()->n
00184        << matout.Mat()->m << matout.Mat()->n << std::endl;
00185   return MatrixMeschach( matout );
00186 
00187 }
00188  
00189 
00190 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00191 MatrixMeschach operator-( const MatrixMeschach& mat1, const MatrixMeschach& mat2 )
00192 {
00193   MatrixMeschach matout( mat1 );
00194   MatrixMeschach matout2( mat2 );
00195   matout += (-1 * matout2);
00196   return MatrixMeschach( matout );
00197 }
00198 
00199 
00200 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00201 MatrixMeschach operator*( const ALIdouble doub, const MatrixMeschach& mat )
00202 {
00203   MatrixMeschach matout( mat );
00204   matout *= doub;
00205   return matout;
00206 }
00207 
00208  
00209 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00210 MatrixMeschach operator*( const MatrixMeschach& mat, const ALIdouble doub )
00211 {
00212   return doub*mat;
00213 }
00214 
00215 
00216 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00217 void MatrixMeschach::transpose()
00218 {
00219   //   if(ALIUtils::debug >= 9)
00220   /*  std::cout << "transposed"  <<_NoLines<<_NoColumns;
00221   MAT* tempmat = m_get(_NoColumns, _NoLines);
00222   m_transp( _Mat, tempmat );
00223   std::cout << "transposed"  <<_NoLines<<_NoColumns;
00224   M_FREE( _Mat );
00225   _Mat = m_get(_NoColumns, _NoLines);
00226   std::cout << "transposed"  <<_NoLines<<_NoColumns;
00227   m_copy( tempmat, _Mat );
00228   std::cout << "transposed"  <<_NoLines<<_NoColumns;
00229   int ntemp = _NoLines;
00230   _NoLines = _NoColumns;
00231   _NoColumns = ntemp;
00232   M_FREE(tempmat);
00233   */
00234 
00235   //-  std::cout << "transposed"  <<_NoLines<<_NoColumns;
00236   MAT* tempmat = m_get(_NoColumns, _NoLines);
00237   m_transp( _Mat, tempmat );
00238   //-  std::cout << "transposed"  <<_NoLines<<_NoColumns;
00239   M_FREE( _Mat );
00240   _Mat = m_get(_NoColumns, _NoLines);
00241   //- std::cout << "transposed"  <<_NoLines<<_NoColumns;
00242   for( ALIint lin=0; lin < _NoColumns; lin++ ) {
00243     for( ALIint col=0;  col < _NoLines; col++ ) {
00244       //-  std::cout << "setting mat "  << lin << " " << col << std::endl;
00245       _Mat->me[lin][col] = tempmat->me[lin][col];
00246     }
00247   }
00248   //  m_copy( tempmat, _Mat );
00249   //-  std::cout << "transposed"  <<_NoLines<<_NoColumns;
00250   int ntemp = _NoLines;
00251   _NoLines = _NoColumns;
00252   _NoColumns = ntemp;
00253   M_FREE(tempmat);
00254  
00255 }
00256 
00257 
00258 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00259 void MatrixMeschach::inverse()
00260 {
00261   if(ALIUtils::debug >= 9) std::cout << "inverse" << _NoLines << "C" << _NoColumns << std::endl;
00262    MAT* tempmat = m_get(_NoLines, _NoColumns);
00263    ALIdouble factor = 1000.;
00264    (*this) *= 1./factor;
00265    m_inverse( _Mat, tempmat );
00266    m_copy( tempmat, _Mat );
00267    (*this) *= 1./factor;
00268   M_FREE(tempmat);
00269 }
00270 
00271 
00272 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00273 void MatrixMeschach::AddData( ALIuint lin, ALIuint col, ALIdouble data ) 
00274 {
00275   if ( lin >= _Mat->m || col >= _Mat->n ) {
00276       std::cerr << "EXITING: matrix has only " << _NoLines << " lines and "
00277            << _NoColumns << " columns " << std::endl;
00278       std::cerr << "EXITING: matrix has only " << _Mat->m << " lines and "
00279            << _Mat->n << " columns " << std::endl;
00280       std::cerr << " You tried to add data in line " << lin << " and column "
00281            << col << std::endl;
00282       std::exception();
00283   }
00284   _Mat->me[lin][col] = data;
00285 
00286 } 
00287 
00288 
00289 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00290 ALIdouble MatrixMeschach::operator () (int i, int j) const 
00291 {
00292   return _Mat->me[i][j];
00293 }  
00294 
00295 
00296 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00297 //@@ 
00298 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00299 void MatrixMeschach::EliminateLines( ALIint lin_first, ALIint lin_last )
00300 {
00301   //-  return;
00302   if ( lin_last < lin_first ) {
00303        std::cerr << "EXITING: cannot Eliminate Lines in matrix if first line is " << 
00304       lin_first << " and lastt line is " << lin_last << std::endl;
00305       //t          std::exception();
00306        return;
00307   }
00308   ALIint dif = (lin_last - lin_first) + 1; 
00309   ALIint newANolin = NoLines() - dif; 
00310   ALIint newANocol = NoColumns();
00311   MatrixMeschach newmat( newANolin, newANocol );
00312   ALIint iin = 0;
00313   for ( ALIint ii=0; ii<NoLines(); ii++) {
00314       if( ii < lin_first  || ii > lin_last ) { 
00315           for ( ALIint jj=0; jj<NoColumns(); jj++) {
00316               newmat.AddData(iin, jj, (*this)(ii,jj) );
00317                 if(ALIUtils::debug >= 9) std::cout << iin << jj << "nmat" << newmat.Mat()->me[iin][jj] << std::endl;
00318           }
00319           iin++;
00320       }
00321   }
00322   M_FREE( _Mat );
00323   _Mat = m_get( newANolin, newANocol );
00324   copy( newmat );
00325   _NoLines = _Mat->m; 
00326   _NoColumns = _Mat->n;
00327   Dump( "newnewmat" );
00328 
00329 }
00330 
00331 
00332 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00333 void MatrixMeschach::EliminateColumns( ALIint lin_first, ALIint lin_last )
00334 {
00335   //-  return;
00336   if ( lin_last < lin_first ) {
00337       std::cerr << "EXITING: cannot Eliminate Lines in matrix if first line is " << 
00338       lin_first << " and lastt line is " << lin_last << std::endl;
00339       //t      std::exception();
00340       return;
00341   }
00342   ALIint dif = (lin_last - lin_first) + 1; 
00343   ALIint newANolin = NoLines();
00344   ALIint newANocol = NoColumns() - dif;
00345   MatrixMeschach newmat( newANolin, newANocol );
00346   ALIint jjn = 0;
00347   for ( ALIint jj=0; jj<NoColumns(); jj++) {
00348       if( jj < lin_first  || jj > lin_last ) { 
00349           for ( ALIint ii=0; ii<NoLines(); ii++) {
00350               newmat.AddData( ii, jjn, (*this)(ii,jj) );
00351                if(ALIUtils::debug >= 9) std::cout << ii << jjn << "nmat" << newmat.Mat()->me[ii][jjn] << std::endl;
00352           }
00353           jjn++;
00354       }
00355   }
00356   M_FREE( _Mat );
00357   _Mat = m_get( newANolin, newANocol );
00358   copy( newmat );
00359   _NoLines = _Mat->m; 
00360   _NoColumns = _Mat->n;
00361   Dump( "newnewmat" );
00362 }
00363 
00364 
00365 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00366 void MatrixMeschach::Dump( const ALIstring& mtext )
00367 {
00368   ostrDump( std::cout, mtext); 
00369 }
00370 
00371 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00372 void MatrixMeschach::ostrDump( std::ostream& fout, const ALIstring& mtext )
00373 {
00374   fout << "DUMPM@@@@@    " << mtext << "    @@@@@" << std::endl;
00375   fout << "Matrix is (_Mat)" << _Mat->m << "x" << _Mat->n << std::endl;
00376   fout << "Matrix is " << _NoLines << "x" << _NoColumns << std::endl;
00377   for (ALIuint ii=0; ii<_Mat->m; ii++) {
00378     for (ALIuint jj=0; jj<_Mat->n; jj++) {
00379       fout << std::setw(8) << _Mat->me[ii][jj] << " ";
00380     }
00381     fout << std::endl;
00382   }
00383   //  m_output(_Mat);
00384 
00385 }
00386  
00387 
00388 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
00389 void MatrixMeschach::SetCorrelation(ALIint i1, ALIint i2, ALIdouble corr)
00390 {  
00391   AddData(i1,i2,corr * sqrt( (*this)(i1,i1)*(*this)(i2,i2) ) );
00392   AddData(i2,i1,corr * sqrt( (*this)(i1,i1)*(*this)(i2,i2) ) );
00393    if(ALIUtils::debug >= 9) std::cout << i1<< i2<< corr << "CORR" << (*this)(i1,i1) << " " << (*this)(i2,i2) << std::endl;
00394   //-  std::cout << (*this)(i1,i1)*(*this)(i2,i2)  << std::endl;
00395   //- std::cout << sqrt( (*this)(i1,i1)*(*this)(i2,i2) ) << std::endl;
00396    if(ALIUtils::debug >= 9) std::cout << corr * sqrt( (*this)(i1,i1)*(*this)(i2,i2) ) << std::endl;
00397 
00398 }
00399 
00400 
00401 MatrixMeschach* MatrixByMatrix( const MatrixMeschach& mat1, const MatrixMeschach& mat2 )
00402 {
00403  MatrixMeschach* matout = new MatrixMeschach( mat1 );
00404   if( matout->NoColumns() != mat2.NoLines() ){
00405     std::cerr << " !! Trying two multiply two matrices when the number of columns of first one is not equal to number of files of second one " << std::endl;
00406     //    std::cout << " multiplying matrix " << matout->NoLines() << " x " << matout->NoColumns() << " and " << mat2.NoLines() << " x " << mat2.NoColumns() << " results in " << matout->NoLines() << " x " << mat2.NoColumns() << std::endl;
00407   }
00408   matout->setNoColumns( mat2.NoColumns() );
00409 
00410   MAT* tempmat = m_get( matout->NoColumns(), matout->NoLines() );
00411   m_transp( matout->MatNonConst(), tempmat); 
00412   //M_FREE( _Mat );
00413   matout->setMat( m_get( matout->NoLines(), mat2.NoColumns() ) );
00414   mtrm_mlt( tempmat, mat2.MatNonConst(), matout->MatNonConst());
00415 
00416 
00417   /*  if(ALIUtils::debug >= 9) std::cout << "add" << matout->NoLines() << matout->NoColumns()
00418        << mat2.NoLines() << mat2.NoColumns()  
00419        << matout->NoLines() << matout->NoColumns() << std::endl;
00420   if(ALIUtils::debug >= 9) std::cout << "addM" << matout->Mat()->m << matout->Mat()->n
00421                                      << mat2.Mat()->m << mat2.Mat()->n
00422        << matout->Mat()->m << matout->Mat()->n << std::endl;
00423   */
00424   //  return MatrixMeschach( matout );
00425   return ( matout );
00426 
00427 }
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/Alignment/CocoaFit/src/MatrixMeschach.cc
