CMSSW Reference Manual

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