CocoaFit/src/MatrixMeschach.cc

0001 //   COCOA class implementation file
0002 //Id:  MatrixMeschach.C
0003 //CAT: Model
0004 //
0005 //   History: v1.0
0006 //   Pedro Arce
0007 #include <iomanip>
0008 #include <cmath>
0009
0010 #include "Alignment/CocoaUtilities/interface/ALIUtils.h"
0011 #include "Alignment/CocoaFit/interface/MatrixMeschach.h"
0012
0013 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0014 MatrixMeschach::MatrixMeschach() {
0015   //-  std::cout << "creating matrix0 " << std::endl;
0016 }
0017
0018 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0019 MatrixMeschach::~MatrixMeschach() {
0020   //-  std::cout << "deleting matrix " << std::endl;
0021   M_FREE(_Mat);
0022 }
0023
0024 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0025 MatrixMeschach::MatrixMeschach(ALIint NoLin, ALIint NoCol) {
0026   //-  std::cout << "creating matrix " << std::endl;
0027   _NoLines = NoLin;
0028   _NoColumns = NoCol;
0029   _Mat = m_get(NoLin, NoCol);
0030   //_data = new ALIdouble[NoCol][NoLin];
0031 }
0032
0033 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0034 MatrixMeschach::MatrixMeschach(const MatrixMeschach& mat) {
0035   //-  std::cout << "creating matrixc " << std::endl;
0036   _NoLines = mat._Mat->m;
0037   _NoColumns = mat._Mat->n;
0038   _Mat = m_get(mat.NoLines(), mat.NoColumns());
0039   // std::cout <<  "copy const" << mat._Mat << _Mat << NoLines() << NoColumns() << Mat()->m << Mat()->n <<std::endl;
0040   copy(mat);
0041 }
0042
0043 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0044 void MatrixMeschach::copy(const MatrixMeschach& mat) {
0045   //  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;
0046
0047   for (ALIint lin = 0; lin < _NoLines; lin++) {
0048     for (ALIint col = 0; col < _NoColumns; col++) {
0049       _Mat->me[lin][col] = mat.MatNonConst()->me[lin][col];
0050     }
0051   }
0052   //   m_copy( mat._Mat, _Mat );
0053 }
0054
0055 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0056 MatrixMeschach& MatrixMeschach::operator=(const MatrixMeschach& mat) {
0057   if (mat._Mat != _Mat) {
0058     _NoLines = mat._Mat->m;
0059     _NoColumns = mat._Mat->n;
0060     M_FREE(_Mat);
0061     _Mat = m_get(mat.NoLines(), mat.NoColumns());
0062     copy(mat);
0063   }
0064   if (ALIUtils::debug >= 9)
0065     std::cout << "operator=" << mat._Mat << _Mat << NoLines() << NoColumns() << Mat()->m << Mat()->n << std::endl;
0066   return *this;
0067 }
0068
0069 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0070 void MatrixMeschach::operator*=(const MatrixMeschach& mat) {
0071   //  std::cout << " multiply matrices " << std::endl;
0072
0073   if (_NoColumns != mat.NoLines()) {
0074     std::cerr << " !! Trying two multiply two matrices when the number of columns of first one is not equal to number "
0075                  "of files of second one "
0076               << std::endl;
0077     std::cout << " multiplying matrix " << _NoLines << " x " << _NoColumns << " and " << mat.NoLines() << " x "
0078               << mat.NoColumns() << " results in " << _NoLines << " x " << mat.NoColumns() << std::endl;
0079   }
0080   _NoColumns = mat.NoColumns();
0081
0082   MAT* tempmat = m_get(_NoColumns, _NoLines);
0083   m_transp(_Mat, tempmat);
0084   //  M_FREE( _Mat );
0085   _Mat = m_get(_NoLines, mat.NoColumns());
0086   mtrm_mlt(tempmat, mat._Mat, _Mat);
0087
0088   //  _NoColumns = mat.NoColumns();
0089   //  M_FREE(tempmat);
0090 }
0091
0092 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0093 void MatrixMeschach::operator+=(const MatrixMeschach& mat) {
0094   if (_NoLines != mat._NoLines || _NoColumns != mat._NoColumns) {
0095     std::cerr << "!!!! cannot sum two matrices with different size" << std::endl
0096               << "MATRIX 1:" << _NoLines << " X " << _NoColumns << std::endl
0097               << "MATRIX 2:" << mat._NoLines << " X " << mat._NoColumns << std::endl;
0098   }
0099   MAT* tempmat = m_get(_NoColumns, _NoLines);
0100   m_copy(_Mat, tempmat);
0101   M_FREE(_Mat);
0102   _Mat = m_get(_NoLines, _NoColumns);
0103   m_add(tempmat, mat._Mat, _Mat);
0104   M_FREE(tempmat);
0105 }
0106
0107 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0108 void MatrixMeschach::operator*=(const ALIdouble num) {
0109   for (ALIuint ii = 0; ii < _Mat->m; ii++) {
0110     for (ALIuint jj = 0; jj < _Mat->n; jj++) {
0111       _Mat->me[ii][jj] *= num;
0112     }
0113   }
0114 }
0115
0116 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0117 MatrixMeschach operator*(const MatrixMeschach& mat1, const MatrixMeschach& mat2) {
0118   MatrixMeschach mat1copy = mat1;
0119   if (mat1copy.NoColumns() != mat2.NoLines()) {
0120     std::cerr << " !! Trying two multiply two matrices when the number of columns of first one is not equal to number "
0121                  "of files of second one "
0122               << std::endl;
0123     std::cout << " multiplying matrix " << mat1copy.NoLines() << " x " << mat1copy.NoColumns() << " and "
0124               << mat2.NoLines() << " x " << mat2.NoColumns() << " results in " << mat1copy.NoLines() << " x "
0125               << mat2.NoColumns() << std::endl;
0126   }
0127   mat1copy.setNoColumns(mat2.NoColumns());
0128
0129   MAT* tempmat = m_get(mat1copy.NoColumns(), mat1copy.NoLines());
0130   m_transp(mat1copy.MatNonConst(), tempmat);
0131
0132   //M_FREE( _Mat );
0133   mat1copy.setMat(m_get(mat1copy.NoLines(), mat2.NoColumns()));
0134   mtrm_mlt(tempmat, mat2.MatNonConst(), mat1copy.MatNonConst());
0135
0136   free(tempmat);
0137
0138   return mat1copy;
0139
0140   MatrixMeschach* matout = new MatrixMeschach(mat1copy);
0141
0142   /*  if(ALIUtils::debug >= 9) std::cout << "add" << mat1copy.NoLines() << mat1copy.NoColumns()
0143        << mat2.NoLines() << mat2.NoColumns()
0144        << matout.NoLines() << matout.NoColumns() << std::endl;
0145   if(ALIUtils::debug >= 9) std::cout << "addM" << mat1copy.Mat()->m << mat1copy.Mat()->n
0146                      << mat2.Mat()->m << mat2.Mat()->n
0147        << matout.Mat()->m << matout.Mat()->n << std::endl;
0148   */
0149   //  return MatrixMeschach( matout );
0150   return (*matout);
0151 }
0152
0153 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0154 MatrixMeschach operator+(const MatrixMeschach& mat1, const MatrixMeschach& mat2) {
0155   MatrixMeschach matout(mat1);
0156   matout += mat2;
0157   if (ALIUtils::debug >= 9)
0158     std::cout << "add" << mat1.NoLines() << mat1.NoColumns() << mat2.NoLines() << mat2.NoColumns() << matout.NoLines()
0159               << matout.NoColumns() << std::endl;
0160   if (ALIUtils::debug >= 9)
0161     std::cout << "addM" << mat1.Mat()->m << mat1.Mat()->n << mat2.Mat()->m << mat2.Mat()->n << matout.Mat()->m
0162               << matout.Mat()->n << std::endl;
0163   return MatrixMeschach(matout);
0164 }
0165
0166 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0167 MatrixMeschach operator-(const MatrixMeschach& mat1, const MatrixMeschach& mat2) {
0168   MatrixMeschach matout(mat1);
0169   const MatrixMeschach& matout2(mat2);
0170   matout += (-1 * matout2);
0171   return MatrixMeschach(matout);
0172 }
0173
0174 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0175 MatrixMeschach operator*(const ALIdouble doub, const MatrixMeschach& mat) {
0176   MatrixMeschach matout(mat);
0177   matout *= doub;
0178   return matout;
0179 }
0180
0181 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0182 MatrixMeschach operator*(const MatrixMeschach& mat, const ALIdouble doub) { return doub * mat; }
0183
0184 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0185 void MatrixMeschach::transpose() {
0186   //   if(ALIUtils::debug >= 9)
0187   /*  std::cout << "transposed"  <<_NoLines<<_NoColumns;
0188   MAT* tempmat = m_get(_NoColumns, _NoLines);
0189   m_transp( _Mat, tempmat );
0190   std::cout << "transposed"  <<_NoLines<<_NoColumns;
0191   M_FREE( _Mat );
0192   _Mat = m_get(_NoColumns, _NoLines);
0193   std::cout << "transposed"  <<_NoLines<<_NoColumns;
0194   m_copy( tempmat, _Mat );
0195   std::cout << "transposed"  <<_NoLines<<_NoColumns;
0196   int ntemp = _NoLines;
0197   _NoLines = _NoColumns;
0198   _NoColumns = ntemp;
0199   M_FREE(tempmat);
0200   */
0201
0202   //-  std::cout << "transposed"  <<_NoLines<<_NoColumns;
0203   MAT* tempmat = m_get(_NoColumns, _NoLines);
0204   m_transp(_Mat, tempmat);
0205   //-  std::cout << "transposed"  <<_NoLines<<_NoColumns;
0206   M_FREE(_Mat);
0207   _Mat = m_get(_NoColumns, _NoLines);
0208   //- std::cout << "transposed"  <<_NoLines<<_NoColumns;
0209   for (ALIint lin = 0; lin < _NoColumns; lin++) {
0210     for (ALIint col = 0; col < _NoLines; col++) {
0211       //-  std::cout << "setting mat "  << lin << " " << col << std::endl;
0212       _Mat->me[lin][col] = tempmat->me[lin][col];
0213     }
0214   }
0215   //  m_copy( tempmat, _Mat );
0216   //-  std::cout << "transposed"  <<_NoLines<<_NoColumns;
0217   int ntemp = _NoLines;
0218   _NoLines = _NoColumns;
0219   _NoColumns = ntemp;
0220   M_FREE(tempmat);
0221 }
0222
0223 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0224 void MatrixMeschach::inverse() {
0225   if (ALIUtils::debug >= 9)
0226     std::cout << "inverse" << _NoLines << "C" << _NoColumns << std::endl;
0227   MAT* tempmat = m_get(_NoLines, _NoColumns);
0228   ALIdouble factor = 1000.;
0229   (*this) *= 1. / factor;
0230   m_inverse(_Mat, tempmat);
0231   m_copy(tempmat, _Mat);
0232   (*this) *= 1. / factor;
0233   M_FREE(tempmat);
0234 }
0235
0236 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0237 void MatrixMeschach::AddData(ALIuint lin, ALIuint col, ALIdouble data) {
0238   if (lin >= _Mat->m || col >= _Mat->n) {
0239     std::cerr << "EXITING: matrix has only " << _NoLines << " lines and " << _NoColumns << " columns " << std::endl;
0240     std::cerr << "EXITING: matrix has only " << _Mat->m << " lines and " << _Mat->n << " columns " << std::endl;
0241     std::cerr << " You tried to add data in line " << lin << " and column " << col << std::endl;
0242     std::exception();
0243   }
0244   _Mat->me[lin][col] = data;
0245 }
0246
0247 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0248 ALIdouble MatrixMeschach::operator()(int i, int j) const { return _Mat->me[i][j]; }
0249
0250 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0251 //@@
0252 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0253 void MatrixMeschach::EliminateLines(ALIint lin_first, ALIint lin_last) {
0254   //-  return;
0255   if (lin_last < lin_first) {
0256     std::cerr << "EXITING: cannot Eliminate Lines in matrix if first line is " << lin_first << " and lastt line is "
0257               << lin_last << std::endl;
0258     //t          std::exception();
0259     return;
0260   }
0261   ALIint dif = (lin_last - lin_first) + 1;
0262   ALIint newANolin = NoLines() - dif;
0263   ALIint newANocol = NoColumns();
0264   MatrixMeschach newmat(newANolin, newANocol);
0265   ALIint iin = 0;
0266   for (ALIint ii = 0; ii < NoLines(); ii++) {
0267     if (ii < lin_first || ii > lin_last) {
0268       for (ALIint jj = 0; jj < NoColumns(); jj++) {
0269         newmat.AddData(iin, jj, (*this)(ii, jj));
0270         if (ALIUtils::debug >= 9)
0271           std::cout << iin << jj << "nmat" << newmat.Mat()->me[iin][jj] << std::endl;
0272       }
0273       iin++;
0274     }
0275   }
0276   M_FREE(_Mat);
0277   _Mat = m_get(newANolin, newANocol);
0278   copy(newmat);
0279   _NoLines = _Mat->m;
0280   _NoColumns = _Mat->n;
0281   Dump("newnewmat");
0282 }
0283
0284 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0285 void MatrixMeschach::EliminateColumns(ALIint lin_first, ALIint lin_last) {
0286   //-  return;
0287   if (lin_last < lin_first) {
0288     std::cerr << "EXITING: cannot Eliminate Lines in matrix if first line is " << lin_first << " and lastt line is "
0289               << lin_last << std::endl;
0290     //t      std::exception();
0291     return;
0292   }
0293   ALIint dif = (lin_last - lin_first) + 1;
0294   ALIint newANolin = NoLines();
0295   ALIint newANocol = NoColumns() - dif;
0296   MatrixMeschach newmat(newANolin, newANocol);
0297   ALIint jjn = 0;
0298   for (ALIint jj = 0; jj < NoColumns(); jj++) {
0299     if (jj < lin_first || jj > lin_last) {
0300       for (ALIint ii = 0; ii < NoLines(); ii++) {
0301         newmat.AddData(ii, jjn, (*this)(ii, jj));
0302         if (ALIUtils::debug >= 9)
0303           std::cout << ii << jjn << "nmat" << newmat.Mat()->me[ii][jjn] << std::endl;
0304       }
0305       jjn++;
0306     }
0307   }
0308   M_FREE(_Mat);
0309   _Mat = m_get(newANolin, newANocol);
0310   copy(newmat);
0311   _NoLines = _Mat->m;
0312   _NoColumns = _Mat->n;
0313   Dump("newnewmat");
0314 }
0315
0316 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0317 void MatrixMeschach::Dump(const ALIstring& mtext) { ostrDump(std::cout, mtext); }
0318
0319 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0320 void MatrixMeschach::ostrDump(std::ostream& fout, const ALIstring& mtext) {
0321   fout << "DUMPM@@@@@    " << mtext << "    @@@@@" << std::endl;
0322   fout << "Matrix is (_Mat)" << _Mat->m << "x" << _Mat->n << std::endl;
0323   fout << "Matrix is " << _NoLines << "x" << _NoColumns << std::endl;
0324   for (ALIuint ii = 0; ii < _Mat->m; ii++) {
0325     for (ALIuint jj = 0; jj < _Mat->n; jj++) {
0326       fout << std::setw(8) << _Mat->me[ii][jj] << " ";
0327     }
0328     fout << std::endl;
0329   }
0330   //  m_output(_Mat);
0331 }
0332
0333 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
0334 void MatrixMeschach::SetCorrelation(ALIint i1, ALIint i2, ALIdouble corr) {
0335   AddData(i1, i2, corr * sqrt((*this)(i1, i1) * (*this)(i2, i2)));
0336   AddData(i2, i1, corr * sqrt((*this)(i1, i1) * (*this)(i2, i2)));
0337   if (ALIUtils::debug >= 9)
0338     std::cout << i1 << i2 << corr << "CORR" << (*this)(i1, i1) << " " << (*this)(i2, i2) << std::endl;
0339   //-  std::cout << (*this)(i1,i1)*(*this)(i2,i2)  << std::endl;
0340   //- std::cout << sqrt( (*this)(i1,i1)*(*this)(i2,i2) ) << std::endl;
0341   if (ALIUtils::debug >= 9)
0342     std::cout << corr * sqrt((*this)(i1, i1) * (*this)(i2, i2)) << std::endl;
0343 }
0344
0345 MatrixMeschach* MatrixByMatrix(const MatrixMeschach& mat1, const MatrixMeschach& mat2) {
0346   MatrixMeschach* matout = new MatrixMeschach(mat1);
0347   if (matout->NoColumns() != mat2.NoLines()) {
0348     std::cerr << " !! Trying two multiply two matrices when the number of columns of first one is not equal to number "
0349                  "of files of second one "
0350               << std::endl;
0351     //    std::cout << " multiplying matrix " << matout->NoLines() << " x " << matout->NoColumns() << " and " << mat2.NoLines() << " x " << mat2.NoColumns() << " results in " << matout->NoLines() << " x " << mat2.NoColumns() << std::endl;
0352   }
0353   matout->setNoColumns(mat2.NoColumns());
0354
0355   MAT* tempmat = m_get(matout->NoColumns(), matout->NoLines());
0356   m_transp(matout->MatNonConst(), tempmat);
0357   //M_FREE( _Mat );
0358   matout->setMat(m_get(matout->NoLines(), mat2.NoColumns()));
0359   mtrm_mlt(tempmat, mat2.MatNonConst(), matout->MatNonConst());
0360
0361   /*  if(ALIUtils::debug >= 9) std::cout << "add" << matout->NoLines() << matout->NoColumns()
0362        << mat2.NoLines() << mat2.NoColumns()
0363        << matout->NoLines() << matout->NoColumns() << std::endl;
0364   if(ALIUtils::debug >= 9) std::cout << "addM" << matout->Mat()->m << matout->Mat()->n
0365                      << mat2.Mat()->m << mat2.Mat()->n
0366        << matout->Mat()->m << matout->Mat()->n << std::endl;
0367   */
0368   //  return MatrixMeschach( matout );
0369   return (matout);
0370 }