MatrixMeschach.cc

Go to the documentation of this file.

//   COCOA class implementation file
//Id:  MatrixMeschach.C
//CAT: Model
//
//   History: v1.0 
//   Pedro Arce
#include <iomanip>
#include <cmath>  


#include "Alignment/CocoaUtilities/interface/ALIUtils.h"
#include "Alignment/CocoaFit/interface/MatrixMeschach.h"


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
MatrixMeschach::MatrixMeschach()
{
  //-  std::cout << "creating matrix0 " << std::endl;
}

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
MatrixMeschach::~MatrixMeschach()
{
  //-  std::cout << "deleting matrix " << std::endl;
  M_FREE(_Mat);
}

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
MatrixMeschach::MatrixMeschach( ALIint NoLin, ALIint NoCol )
{
  //-  std::cout << "creating matrix " << std::endl;
  _NoLines = NoLin; 
  _NoColumns = NoCol;
  _Mat = m_get( NoLin, NoCol );
  //_data = new ALIdouble[NoCol][NoLin];
}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
MatrixMeschach::MatrixMeschach( const MatrixMeschach& mat )
{
  //-  std::cout << "creating matrixc " << std::endl;
  _NoLines = mat._Mat->m; 
  _NoColumns = mat._Mat->n;
  _Mat = m_get( mat.NoLines(), mat.NoColumns() );
  // std::cout <<  "copy const" << mat._Mat << _Mat << NoLines() << NoColumns() << Mat()->m << Mat()->n <<std::endl;
  copy( mat );
}
 

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::copy( const MatrixMeschach& mat) 
{
  //  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;

  for( ALIint lin=0; lin < _NoLines; lin++ ) {
    for( ALIint col=0;  col < _NoColumns; col++ ) {
      _Mat->me[lin][col] = mat.MatNonConst()->me[lin][col];
    } 
  }
  //   m_copy( mat._Mat, _Mat );
} 


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
MatrixMeschach& MatrixMeschach::operator=( const MatrixMeschach& mat )
{
  if ( mat._Mat != _Mat ) {
      _NoLines = mat._Mat->m;
      _NoColumns = mat._Mat->n; 
      M_FREE( _Mat );
      _Mat = m_get( mat.NoLines(), mat.NoColumns() );
      copy( mat );
  }
   if(ALIUtils::debug >= 9) std::cout <<  "operator=" << mat._Mat << _Mat << NoLines() << NoColumns() << Mat()->m << Mat()->n <<std::endl;
  return *this;
}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::operator*=( const MatrixMeschach& mat )
{
  //  std::cout << " multiply matrices " << std::endl;

  if( _NoColumns != mat.NoLines() ){
    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;
    std::cout << " multiplying matrix " << _NoLines  << " x " << _NoColumns << " and " << mat.NoLines() << " x " << mat.NoColumns() << " results in " << _NoLines << " x " << mat.NoColumns() << std::endl;
  }
  _NoColumns = mat.NoColumns();

  MAT* tempmat = m_get( _NoColumns, _NoLines );
  m_transp( _Mat, tempmat); 
  //  M_FREE( _Mat );
  _Mat = m_get( _NoLines, mat.NoColumns() );
  mtrm_mlt( tempmat, mat._Mat, _Mat);

  //  _NoColumns = mat.NoColumns(); 
  //  M_FREE(tempmat);
}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::operator+=( const MatrixMeschach& mat )
{
  if(_NoLines != mat._NoLines  || _NoColumns != mat._NoColumns ) {
    std::cerr << "!!!! cannot sum two matrices with different size" << std::endl
         << "MATRIX 1:" << _NoLines << " X " << _NoColumns << std::endl
         << "MATRIX 2:" << mat._NoLines << " X " << mat._NoColumns << std::endl;
  }
  MAT* tempmat = m_get( _NoColumns, _NoLines );
  m_copy( _Mat, tempmat); 
  M_FREE( _Mat );
  _Mat = m_get( _NoLines, _NoColumns );
  m_add( tempmat, mat._Mat, _Mat);
  M_FREE(tempmat);

}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::operator*=( const ALIdouble num )
{
  for (ALIuint ii=0; ii<_Mat->m; ii++) {
      for (ALIuint jj=0; jj<_Mat->n; jj++) {
      _Mat->me[ii][jj] *= num;
      }
  }

}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
MatrixMeschach operator*( const MatrixMeschach& mat1, const MatrixMeschach& mat2 )
{
  MatrixMeschach mat1copy = mat1;
  if( mat1copy.NoColumns() != mat2.NoLines() ){
    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;
    std::cout << " multiplying matrix " << mat1copy.NoLines() << " x " << mat1copy.NoColumns() << " and " << mat2.NoLines() << " x " << mat2.NoColumns() << " results in " << mat1copy.NoLines() << " x " << mat2.NoColumns() << std::endl;
  }
  mat1copy.setNoColumns( mat2.NoColumns() );

  MAT* tempmat = m_get( mat1copy.NoColumns(), mat1copy.NoLines() );
  m_transp( mat1copy.MatNonConst(), tempmat);
 
  //M_FREE( _Mat );
  mat1copy.setMat( m_get( mat1copy.NoLines(), mat2.NoColumns() ) );
  mtrm_mlt( tempmat, mat2.MatNonConst(), mat1copy.MatNonConst());

  free(tempmat);

  return mat1copy;


  MatrixMeschach* matout = new MatrixMeschach( mat1copy );

  /*  if(ALIUtils::debug >= 9) std::cout << "add" << mat1copy.NoLines() << mat1copy.NoColumns()
       << mat2.NoLines() << mat2.NoColumns()  
       << matout.NoLines() << matout.NoColumns() << std::endl;
  if(ALIUtils::debug >= 9) std::cout << "addM" << mat1copy.Mat()->m << mat1copy.Mat()->n
                     << mat2.Mat()->m << mat2.Mat()->n
       << matout.Mat()->m << matout.Mat()->n << std::endl;
  */
  //  return MatrixMeschach( matout );
  return ( *matout );

}
 

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
MatrixMeschach operator+( const MatrixMeschach& mat1, const MatrixMeschach& mat2 )
{
  MatrixMeschach matout( mat1 );
  matout += mat2;
   if(ALIUtils::debug >= 9) std::cout << "add" << mat1.NoLines() << mat1.NoColumns()
       << mat2.NoLines() << mat2.NoColumns()  
       << matout.NoLines() << matout.NoColumns() << std::endl;
   if(ALIUtils::debug >= 9) std::cout << "addM" << mat1.Mat()->m << mat1.Mat()->n
       << mat2.Mat()->m << mat2.Mat()->n
       << matout.Mat()->m << matout.Mat()->n << std::endl;
  return MatrixMeschach( matout );

}
 

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
MatrixMeschach operator-( const MatrixMeschach& mat1, const MatrixMeschach& mat2 )
{
  MatrixMeschach matout( mat1 );
  const MatrixMeschach& matout2( mat2 );
  matout += (-1 * matout2);
  return MatrixMeschach( matout );
}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
MatrixMeschach operator*( const ALIdouble doub, const MatrixMeschach& mat )
{
  MatrixMeschach matout( mat );
  matout *= doub;
  return matout;
}

 
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
MatrixMeschach operator*( const MatrixMeschach& mat, const ALIdouble doub )
{
  return doub*mat;
}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::transpose()
{
  //   if(ALIUtils::debug >= 9)
  /*  std::cout << "transposed"  <<_NoLines<<_NoColumns;
  MAT* tempmat = m_get(_NoColumns, _NoLines);
  m_transp( _Mat, tempmat );
  std::cout << "transposed"  <<_NoLines<<_NoColumns;
  M_FREE( _Mat );
  _Mat = m_get(_NoColumns, _NoLines);
  std::cout << "transposed"  <<_NoLines<<_NoColumns;
  m_copy( tempmat, _Mat );
  std::cout << "transposed"  <<_NoLines<<_NoColumns;
  int ntemp = _NoLines;
  _NoLines = _NoColumns;
  _NoColumns = ntemp;
  M_FREE(tempmat);
  */

  //-  std::cout << "transposed"  <<_NoLines<<_NoColumns;
  MAT* tempmat = m_get(_NoColumns, _NoLines);
  m_transp( _Mat, tempmat );
  //-  std::cout << "transposed"  <<_NoLines<<_NoColumns;
  M_FREE( _Mat );
  _Mat = m_get(_NoColumns, _NoLines);
  //- std::cout << "transposed"  <<_NoLines<<_NoColumns;
  for( ALIint lin=0; lin < _NoColumns; lin++ ) {
    for( ALIint col=0;  col < _NoLines; col++ ) {
      //-  std::cout << "setting mat "  << lin << " " << col << std::endl;
      _Mat->me[lin][col] = tempmat->me[lin][col];
    }
  }
  //  m_copy( tempmat, _Mat );
  //-  std::cout << "transposed"  <<_NoLines<<_NoColumns;
  int ntemp = _NoLines;
  _NoLines = _NoColumns;
  _NoColumns = ntemp;
  M_FREE(tempmat);
 
}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::inverse()
{
  if(ALIUtils::debug >= 9) std::cout << "inverse" << _NoLines << "C" << _NoColumns << std::endl;
   MAT* tempmat = m_get(_NoLines, _NoColumns);
   ALIdouble factor = 1000.;
   (*this) *= 1./factor;
   m_inverse( _Mat, tempmat );
   m_copy( tempmat, _Mat );
   (*this) *= 1./factor;
  M_FREE(tempmat);
}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::AddData( ALIuint lin, ALIuint col, ALIdouble data ) 
{
  if ( lin >= _Mat->m || col >= _Mat->n ) {
      std::cerr << "EXITING: matrix has only " << _NoLines << " lines and "
           << _NoColumns << " columns " << std::endl;
      std::cerr << "EXITING: matrix has only " << _Mat->m << " lines and "
           << _Mat->n << " columns " << std::endl;
      std::cerr << " You tried to add data in line " << lin << " and column "
           << col << std::endl;
      std::exception();
  }
  _Mat->me[lin][col] = data;

} 


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
ALIdouble MatrixMeschach::operator () (int i, int j) const 
{
  return _Mat->me[i][j];
}  


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//@@ 
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::EliminateLines( ALIint lin_first, ALIint lin_last )
{
  //-  return;
  if ( lin_last < lin_first ) {
       std::cerr << "EXITING: cannot Eliminate Lines in matrix if first line is " << 
      lin_first << " and lastt line is " << lin_last << std::endl;
      //t          std::exception();
       return;
  }
  ALIint dif = (lin_last - lin_first) + 1; 
  ALIint newANolin = NoLines() - dif; 
  ALIint newANocol = NoColumns();
  MatrixMeschach newmat( newANolin, newANocol );
  ALIint iin = 0;
  for ( ALIint ii=0; ii<NoLines(); ii++) {
      if( ii < lin_first  || ii > lin_last ) { 
          for ( ALIint jj=0; jj<NoColumns(); jj++) {
              newmat.AddData(iin, jj, (*this)(ii,jj) );
            if(ALIUtils::debug >= 9) std::cout << iin << jj << "nmat" << newmat.Mat()->me[iin][jj] << std::endl;
      }
          iin++;
      }
  }
  M_FREE( _Mat );
  _Mat = m_get( newANolin, newANocol );
  copy( newmat );
  _NoLines = _Mat->m; 
  _NoColumns = _Mat->n;
  Dump( "newnewmat" );

}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::EliminateColumns( ALIint lin_first, ALIint lin_last )
{
  //-  return;
  if ( lin_last < lin_first ) {
      std::cerr << "EXITING: cannot Eliminate Lines in matrix if first line is " << 
      lin_first << " and lastt line is " << lin_last << std::endl;
      //t      std::exception();
      return;
  }
  ALIint dif = (lin_last - lin_first) + 1; 
  ALIint newANolin = NoLines();
  ALIint newANocol = NoColumns() - dif;
  MatrixMeschach newmat( newANolin, newANocol );
  ALIint jjn = 0;
  for ( ALIint jj=0; jj<NoColumns(); jj++) {
      if( jj < lin_first  || jj > lin_last ) { 
          for ( ALIint ii=0; ii<NoLines(); ii++) {
              newmat.AddData( ii, jjn, (*this)(ii,jj) );
               if(ALIUtils::debug >= 9) std::cout << ii << jjn << "nmat" << newmat.Mat()->me[ii][jjn] << std::endl;
      }
          jjn++;
      }
  }
  M_FREE( _Mat );
  _Mat = m_get( newANolin, newANocol );
  copy( newmat );
  _NoLines = _Mat->m; 
  _NoColumns = _Mat->n;
  Dump( "newnewmat" );
}


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::Dump( const ALIstring& mtext )
{
  ostrDump( std::cout, mtext); 
}

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::ostrDump( std::ostream& fout, const ALIstring& mtext )
{
  fout << "DUMPM@@@@@    " << mtext << "    @@@@@" << std::endl;
  fout << "Matrix is (_Mat)" << _Mat->m << "x" << _Mat->n << std::endl;
  fout << "Matrix is " << _NoLines << "x" << _NoColumns << std::endl;
  for (ALIuint ii=0; ii<_Mat->m; ii++) {
    for (ALIuint jj=0; jj<_Mat->n; jj++) {
      fout << std::setw(8) << _Mat->me[ii][jj] << " ";
    }
    fout << std::endl;
  }
  //  m_output(_Mat);

}
 

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
void MatrixMeschach::SetCorrelation(ALIint i1, ALIint i2, ALIdouble corr)
{  
  AddData(i1,i2,corr * sqrt( (*this)(i1,i1)*(*this)(i2,i2) ) );
  AddData(i2,i1,corr * sqrt( (*this)(i1,i1)*(*this)(i2,i2) ) );
   if(ALIUtils::debug >= 9) std::cout << i1<< i2<< corr << "CORR" << (*this)(i1,i1) << " " << (*this)(i2,i2) << std::endl;
  //-  std::cout << (*this)(i1,i1)*(*this)(i2,i2)  << std::endl;
  //- std::cout << sqrt( (*this)(i1,i1)*(*this)(i2,i2) ) << std::endl;
   if(ALIUtils::debug >= 9) std::cout << corr * sqrt( (*this)(i1,i1)*(*this)(i2,i2) ) << std::endl;

}


MatrixMeschach* MatrixByMatrix( const MatrixMeschach& mat1, const MatrixMeschach& mat2 )
{
 MatrixMeschach* matout = new MatrixMeschach( mat1 );
  if( matout->NoColumns() != mat2.NoLines() ){
    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;
    //    std::cout << " multiplying matrix " << matout->NoLines() << " x " << matout->NoColumns() << " and " << mat2.NoLines() << " x " << mat2.NoColumns() << " results in " << matout->NoLines() << " x " << mat2.NoColumns() << std::endl;
  }
  matout->setNoColumns( mat2.NoColumns() );

  MAT* tempmat = m_get( matout->NoColumns(), matout->NoLines() );
  m_transp( matout->MatNonConst(), tempmat); 
  //M_FREE( _Mat );
  matout->setMat( m_get( matout->NoLines(), mat2.NoColumns() ) );
  mtrm_mlt( tempmat, mat2.MatNonConst(), matout->MatNonConst());


  /*  if(ALIUtils::debug >= 9) std::cout << "add" << matout->NoLines() << matout->NoColumns()
       << mat2.NoLines() << mat2.NoColumns()  
       << matout->NoLines() << matout->NoColumns() << std::endl;
  if(ALIUtils::debug >= 9) std::cout << "addM" << matout->Mat()->m << matout->Mat()->n
                     << mat2.Mat()->m << mat2.Mat()->n
       << matout->Mat()->m << matout->Mat()->n << std::endl;
  */
  //  return MatrixMeschach( matout );
  return ( matout );

}