CommonAlignmentAlgorithm/src/AlignmentExtendedCorrelationsStore.cc

0001 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentExtendedCorrelationsStore.h"
0002
0003 #include "Alignment/CommonAlignment/interface/Alignable.h"
0004 #include "Alignment/CommonAlignment/interface/AlignmentParameters.h"
0005 #include "Alignment/CommonAlignmentAlgorithm/interface/AlignmentExtendedCorrelationsEntry.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/Utilities/interface/Exception.h"
0009 #include "FWCore/Utilities/interface/isFinite.h"
0010
0011 AlignmentExtendedCorrelationsStore::AlignmentExtendedCorrelationsStore(const edm::ParameterSet& config) {
0012   theMaxUpdates = config.getParameter<int>("MaxUpdates");
0013   theCut = config.getParameter<double>("CutValue");
0014   theWeight = config.getParameter<double>("Weight");
0015
0016   edm::LogInfo("Alignment") << "@SUB=AlignmentExtendedCorrelationsStore::AlignmentExtendedCorrelationsStore"
0017                             << "Created.";
0018 }
0019
0020 void AlignmentExtendedCorrelationsStore::correlations(
0021     Alignable* ap1, Alignable* ap2, AlgebraicSymMatrix& cov, int row, int col) const {
0022   static Alignable* previousAlignable = nullptr;
0023   static ExtendedCorrelationsTable* previousCorrelations;
0024
0025   // Needed by 'resetCorrelations()' to reset the static pointer:
0026   if (ap1 == nullptr) {
0027     previousAlignable = nullptr;
0028     return;
0029   }
0030
0031   bool transpose = (ap2 > ap1);
0032   if (transpose)
0033     std::swap(ap1, ap2);
0034
0035   if (ap1 == previousAlignable) {
0036     ExtendedCorrelationsTable::const_iterator itC2 = previousCorrelations->find(ap2);
0037     if (itC2 != previousCorrelations->end()) {
0038       transpose ? fillCovarianceT(ap1, ap2, (*itC2).second, cov, row, col)
0039                 : fillCovariance(ap1, ap2, (*itC2).second, cov, row, col);
0040     }
0041   } else {
0042     ExtendedCorrelations::const_iterator itC1 = theCorrelations.find(ap1);
0043     if (itC1 != theCorrelations.end()) {
0044       previousAlignable = ap1;
0045       previousCorrelations = (*itC1).second;
0046
0047       ExtendedCorrelationsTable::const_iterator itC2 = (*itC1).second->find(ap2);
0048       if (itC2 != (*itC1).second->end()) {
0049         transpose ? fillCovarianceT(ap1, ap2, (*itC2).second, cov, row, col)
0050                   : fillCovariance(ap1, ap2, (*itC2).second, cov, row, col);
0051       }
0052     }
0053   }
0054
0055   // don't fill anything into the covariance if there's no entry
0056   return;
0057 }
0058
0059 void AlignmentExtendedCorrelationsStore::setCorrelations(
0060     Alignable* ap1, Alignable* ap2, const AlgebraicSymMatrix& cov, int row, int col) {
0061   static Alignable* previousAlignable = nullptr;
0062   static ExtendedCorrelationsTable* previousCorrelations;
0063
0064   // Needed by 'resetCorrelations()' to reset the static pointer:
0065   if (ap1 == nullptr) {
0066     previousAlignable = nullptr;
0067     return;
0068   }
0069
0070   bool transpose = (ap2 > ap1);
0071   if (transpose)
0072     std::swap(ap1, ap2);
0073
0074   if (ap1 == previousAlignable) {
0075     fillCorrelationsTable(ap1, ap2, previousCorrelations, cov, row, col, transpose);
0076   } else {
0077     ExtendedCorrelations::iterator itC = theCorrelations.find(ap1);
0078     if (itC != theCorrelations.end()) {
0079       fillCorrelationsTable(ap1, ap2, itC->second, cov, row, col, transpose);
0080       previousAlignable = ap1;
0081       previousCorrelations = itC->second;
0082     } else {
0083       // make new entry
0084       ExtendedCorrelationsTable* newTable = new ExtendedCorrelationsTable;
0085       fillCorrelationsTable(ap1, ap2, newTable, cov, row, col, transpose);
0086
0087       theCorrelations[ap1] = newTable;
0088
0089       previousAlignable = ap1;
0090       previousCorrelations = newTable;
0091     }
0092   }
0093 }
0094
0095 void AlignmentExtendedCorrelationsStore::setCorrelations(Alignable* ap1, Alignable* ap2, AlgebraicMatrix& mat) {
0096   bool transpose = (ap2 > ap1);
0097   if (transpose)
0098     std::swap(ap1, ap2);
0099
0100   ExtendedCorrelations::iterator itC1 = theCorrelations.find(ap1);
0101   if (itC1 != theCorrelations.end()) {
0102     ExtendedCorrelationsTable::iterator itC2 = itC1->second->find(ap1);
0103     if (itC2 != itC1->second->end()) {
0104       itC2->second = transpose ? ExtendedCorrelationsEntry(mat.T()) : ExtendedCorrelationsEntry(mat);
0105     } else {
0106       (*itC1->second)[ap2] = transpose ? ExtendedCorrelationsEntry(mat.T()) : ExtendedCorrelationsEntry(mat);
0107     }
0108   } else {
0109     ExtendedCorrelationsTable* newTable = new ExtendedCorrelationsTable;
0110     (*newTable)[ap2] = transpose ? ExtendedCorrelationsEntry(mat.T()) : ExtendedCorrelationsEntry(mat);
0111     theCorrelations[ap1] = newTable;
0112   }
0113 }
0114
0115 void AlignmentExtendedCorrelationsStore::getCorrelations(Alignable* ap1, Alignable* ap2, AlgebraicMatrix& mat) const {
0116   bool transpose = (ap2 > ap1);
0117   if (transpose)
0118     std::swap(ap1, ap2);
0119
0120   ExtendedCorrelations::const_iterator itC1 = theCorrelations.find(ap1);
0121   if (itC1 != theCorrelations.end()) {
0122     ExtendedCorrelationsTable::const_iterator itC2 = itC1->second->find(ap2);
0123     if (itC2 != itC1->second->end()) {
0124       mat = transpose ? itC2->second.matrix().T() : itC2->second.matrix();
0125       return;
0126     }
0127   }
0128
0129   mat = AlgebraicMatrix();
0130 }
0131
0132 bool AlignmentExtendedCorrelationsStore::correlationsAvailable(Alignable* ap1, Alignable* ap2) const {
0133   bool transpose = (ap2 > ap1);
0134   if (transpose)
0135     std::swap(ap1, ap2);
0136
0137   ExtendedCorrelations::const_iterator itC1 = theCorrelations.find(ap1);
0138   if (itC1 != theCorrelations.end()) {
0139     ExtendedCorrelationsTable::const_iterator itC2 = itC1->second->find(ap2);
0140     if (itC2 != itC1->second->end())
0141       return true;
0142   }
0143   return false;
0144 }
0145
0146 void AlignmentExtendedCorrelationsStore::resetCorrelations(void) {
0147   ExtendedCorrelations::iterator itC;
0148   for (itC = theCorrelations.begin(); itC != theCorrelations.end(); ++itC)
0149     delete (*itC).second;
0150   theCorrelations.erase(theCorrelations.begin(), theCorrelations.end());
0151
0152   // Reset the static pointers to the 'previous alignables'
0153   AlgebraicSymMatrix dummy;
0154   correlations(nullptr, nullptr, dummy, 0, 0);
0155   setCorrelations(nullptr, nullptr, dummy, 0, 0);
0156 }
0157
0158 unsigned int AlignmentExtendedCorrelationsStore::size(void) const {
0159   unsigned int size = 0;
0160   ExtendedCorrelations::const_iterator itC;
0161   for (itC = theCorrelations.begin(); itC != theCorrelations.end(); ++itC)
0162     size += itC->second->size();
0163
0164   return size;
0165 }
0166
0167 void AlignmentExtendedCorrelationsStore::fillCorrelationsTable(Alignable* ap1,
0168                                                                Alignable* ap2,
0169                                                                ExtendedCorrelationsTable* table,
0170                                                                const AlgebraicSymMatrix& cov,
0171                                                                int row,
0172                                                                int col,
0173                                                                bool transpose) {
0174   ExtendedCorrelationsTable::iterator itC = table->find(ap2);
0175
0176   if (itC != table->end()) {
0177     //if ( itC->second.counter() > theMaxUpdates ) return;
0178
0179     transpose ? readFromCovarianceT(ap1, ap2, itC->second, cov, row, col)
0180               : readFromCovariance(ap1, ap2, itC->second, cov, row, col);
0181
0182     //itC->second.incrementCounter();
0183   } else {
0184     int nRow = ap1->alignmentParameters()->numSelected();
0185     int nCol = ap2->alignmentParameters()->numSelected();
0186     ExtendedCorrelationsEntry newEntry(nRow, nCol);
0187
0188     transpose ? readFromCovarianceT(ap1, ap2, newEntry, cov, row, col)
0189               : readFromCovariance(ap1, ap2, newEntry, cov, row, col);
0190
0191     (*table)[ap2] = newEntry;
0192   }
0193 }
0194
0195 void AlignmentExtendedCorrelationsStore::fillCovariance(Alignable* ap1,
0196                                                         Alignable* ap2,
0197                                                         const ExtendedCorrelationsEntry& entry,
0198                                                         AlgebraicSymMatrix& cov,
0199                                                         int row,
0200                                                         int col) const {
0201   int nRow = entry.numRow();
0202   int nCol = entry.numCol();
0203
0204   for (int iRow = 0; iRow < nRow; ++iRow) {
0205     double factor = sqrt(cov[row + iRow][row + iRow]);
0206     if (edm::isNotFinite(factor))
0207       throw cms::Exception("LogicError") << "[AlignmentExtendedCorrelationsStore::fillCovariance] "
0208                                          << "NaN-factor: sqrt(" << cov[row + iRow][row + iRow] << ")";
0209
0210     for (int jCol = 0; jCol < nCol; ++jCol)
0211       cov[row + iRow][col + jCol] = entry(iRow, jCol) * factor;
0212   }
0213
0214   for (int jCol = 0; jCol < nCol; ++jCol) {
0215     double factor = sqrt(cov[col + jCol][col + jCol]);
0216     if (edm::isNotFinite(factor))
0217       throw cms::Exception("LogicError") << "[AlignmentExtendedCorrelationsStore::fillCovariance] "
0218                                          << "NaN-factor: sqrt(" << cov[col + jCol][col + jCol] << ")";
0219
0220     for (int iRow = 0; iRow < nRow; ++iRow)
0221       cov[row + iRow][col + jCol] *= factor;
0222   }
0223 }
0224
0225 void AlignmentExtendedCorrelationsStore::fillCovarianceT(Alignable* ap1,
0226                                                          Alignable* ap2,
0227                                                          const ExtendedCorrelationsEntry& entry,
0228                                                          AlgebraicSymMatrix& cov,
0229                                                          int row,
0230                                                          int col) const {
0231   int nRow = entry.numRow();
0232   int nCol = entry.numCol();
0233
0234   for (int iRow = 0; iRow < nRow; ++iRow) {
0235     double factor = sqrt(cov[col + iRow][col + iRow]);
0236     if (edm::isNotFinite(factor))
0237       throw cms::Exception("LogicError") << "[AlignmentExtendedCorrelationsStore::fillCovarianceT] "
0238                                          << "NaN-factor: sqrt(" << cov[col + iRow][col + iRow] << ")";
0239     for (int jCol = 0; jCol < nCol; ++jCol)
0240       cov[row + jCol][col + iRow] = entry(iRow, jCol) * factor;
0241   }
0242
0243   for (int jCol = 0; jCol < nCol; ++jCol) {
0244     double factor = sqrt(cov[row + jCol][row + jCol]);
0245     if (edm::isNotFinite(factor))
0246       throw cms::Exception("LogicError") << "[AlignmentExtendedCorrelationsStore::fillCovarianceT] "
0247                                          << "NaN-factor: sqrt(" << cov[row + jCol][row + jCol] << ")";
0248     for (int iRow = 0; iRow < nRow; ++iRow)
0249       cov[row + jCol][col + iRow] *= factor;
0250   }
0251 }
0252
0253 void AlignmentExtendedCorrelationsStore::readFromCovariance(
0254     Alignable* ap1, Alignable* ap2, ExtendedCorrelationsEntry& entry, const AlgebraicSymMatrix& cov, int row, int col) {
0255   int nRow = entry.numRow();
0256   int nCol = entry.numCol();
0257
0258   for (int iRow = 0; iRow < nRow; ++iRow) {
0259     double factor = sqrt(cov[row + iRow][row + iRow]);
0260     for (int jCol = 0; jCol < nCol; ++jCol)
0261       entry(iRow, jCol) = cov[row + iRow][col + jCol] / factor;
0262   }
0263
0264   double maxCorr = 0;
0265
0266   for (int jCol = 0; jCol < nCol; ++jCol) {
0267     double factor = sqrt(cov[col + jCol][col + jCol]);
0268     for (int iRow = 0; iRow < nRow; ++iRow) {
0269       entry(iRow, jCol) /= factor;
0270       if (fabs(entry(iRow, jCol)) > maxCorr)
0271         maxCorr = fabs(entry(iRow, jCol));
0272     }
0273   }
0274
0275   resizeCorruptCorrelations(entry, maxCorr);
0276 }
0277
0278 void AlignmentExtendedCorrelationsStore::readFromCovarianceT(
0279     Alignable* ap1, Alignable* ap2, ExtendedCorrelationsEntry& entry, const AlgebraicSymMatrix& cov, int row, int col) {
0280   int nRow = entry.numRow();
0281   int nCol = entry.numCol();
0282
0283   for (int iRow = 0; iRow < nRow; ++iRow) {
0284     double factor = sqrt(cov[col + iRow][col + iRow]);
0285     for (int jCol = 0; jCol < nCol; ++jCol)
0286       entry(iRow, jCol) = cov[row + jCol][col + iRow] / factor;
0287   }
0288
0289   double maxCorr = 0;
0290
0291   for (int jCol = 0; jCol < nCol; ++jCol) {
0292     double factor = sqrt(cov[row + jCol][row + jCol]);
0293     for (int iRow = 0; iRow < nRow; ++iRow) {
0294       entry(iRow, jCol) /= factor;
0295       if (fabs(entry(iRow, jCol)) > maxCorr)
0296         maxCorr = fabs(entry(iRow, jCol));
0297     }
0298   }
0299
0300   resizeCorruptCorrelations(entry, maxCorr);
0301 }
0302
0303 void AlignmentExtendedCorrelationsStore::resizeCorruptCorrelations(ExtendedCorrelationsEntry& entry, double maxCorr) {
0304   if (maxCorr > 1.) {
0305     entry *= theWeight / maxCorr;
0306   } else if (maxCorr > theCut) {
0307     entry *= 1. - (maxCorr - theCut) / (1. - theCut) * (1. - theWeight);
0308   }
0309 }