File indexing completed on 2024-04-06 11:55:59
0001
0002
0003
0004
0005
0006
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
0016 }
0017
0018
0019 MatrixMeschach::~MatrixMeschach() {
0020
0021 M_FREE(_Mat);
0022 }
0023
0024
0025 MatrixMeschach::MatrixMeschach(ALIint NoLin, ALIint NoCol) {
0026
0027 _NoLines = NoLin;
0028 _NoColumns = NoCol;
0029 _Mat = m_get(NoLin, NoCol);
0030
0031 }
0032
0033
0034 MatrixMeschach::MatrixMeschach(const MatrixMeschach& mat) {
0035
0036 _NoLines = mat._Mat->m;
0037 _NoColumns = mat._Mat->n;
0038 _Mat = m_get(mat.NoLines(), mat.NoColumns());
0039
0040 copy(mat);
0041 }
0042
0043
0044 void MatrixMeschach::copy(const MatrixMeschach& mat) {
0045
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
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
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
0085 _Mat = m_get(_NoLines, mat.NoColumns());
0086 mtrm_mlt(tempmat, mat._Mat, _Mat);
0087
0088
0089
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
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
0143
0144
0145
0146
0147
0148
0149
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
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203 MAT* tempmat = m_get(_NoColumns, _NoLines);
0204 m_transp(_Mat, tempmat);
0205
0206 M_FREE(_Mat);
0207 _Mat = m_get(_NoColumns, _NoLines);
0208
0209 for (ALIint lin = 0; lin < _NoColumns; lin++) {
0210 for (ALIint col = 0; col < _NoLines; col++) {
0211
0212 _Mat->me[lin][col] = tempmat->me[lin][col];
0213 }
0214 }
0215
0216
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
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
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
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
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
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
0340
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
0352 }
0353 matout->setNoColumns(mat2.NoColumns());
0354
0355 MAT* tempmat = m_get(matout->NoColumns(), matout->NoLines());
0356 m_transp(matout->MatNonConst(), tempmat);
0357
0358 matout->setMat(m_get(matout->NoLines(), mat2.NoColumns()));
0359 mtrm_mlt(tempmat, mat2.MatNonConst(), matout->MatNonConst());
0360
0361
0362
0363
0364
0365
0366
0367
0368
0369 return (matout);
0370 }