Tools/src/GenericHouseholder.cc

0001 /** \file GenericHouseholder.cc
0002  *
0003  *
0004  * \author Lorenzo Agostino, R.Ofierzynski, CERN
0005  */
0006
0007 #include "Calibration/Tools/interface/GenericHouseholder.h"
0008 #include <cfloat>
0009 #include <cmath>
0010
0011 GenericHouseholder::GenericHouseholder(bool normalise) : normaliseFlag(normalise) {}
0012
0013 GenericHouseholder::~GenericHouseholder() {}
0014
0015 std::vector<float> GenericHouseholder::iterate(const std::vector<std::vector<float> >& eventMatrix,
0016                                                const std::vector<float>& energyVector,
0017                                                const int nIter) {
0018   std::vector<float> solution;
0019   std::vector<float> theCalibVector(energyVector.size(), 1.);
0020   std::vector<std::vector<float> > myEventMatrix(eventMatrix);
0021   int Nevents = eventMatrix.size();       // Number of events to calibrate with
0022   int Nchannels = eventMatrix[0].size();  // Number of channel coefficients
0023
0024   // Iterate the correction
0025   for (int iter = 1; iter <= nIter; iter++) {
0026     // make one iteration
0027     solution = iterate(myEventMatrix, energyVector);
0028
0029     if (solution.empty())
0030       return solution;
0031     // R.O.: or throw an exception, what's the standard CMS way ?
0032
0033     // re-calibrate eventMatrix with solution
0034     for (int i = 0; i < Nchannels; i++) {
0035       for (int ievent = 0; ievent < Nevents; ievent++) {
0036         myEventMatrix[ievent][i] *= solution[i];
0037       }
0038       // save solution into theCalibVector
0039       theCalibVector[i] *= solution[i];
0040     }
0041
0042   }  // end iterate the correction
0043
0044   return theCalibVector;
0045 }
0046
0047 std::vector<float> GenericHouseholder::iterate(const std::vector<std::vector<float> >& eventMatrix,
0048                                                const std::vector<float>& energyVector) {
0049   // An implementation of the Householder in-situ calibration algorithm
0050   // (Least squares minimisation of residual R=b-Ax, with QR decomposition of A)
0051   // A: matrix of channel response for all calib events
0052   // x: vector of channel calibration coefficients
0053   // b: vector of energies
0054   // adapted from the original code by Matt Probert 9/08/01.
0055
0056   std::vector<float> solution;
0057
0058   unsigned int m = eventMatrix.size();     // Number of events to calibrate with
0059   unsigned int n = eventMatrix[0].size();  // Number of channel coefficients to optimize
0060
0061   std::cout << "Householder::runIter(): starting calibration optimization:" << std::endl;
0062   std::cout << "  Events:" << m << ", channels: " << n << std::endl;
0063
0064   // Sanity check
0065   if (m != energyVector.size()) {
0066     std::cout << "Householder::runIter(): matrix dimensions non-conformant. " << std::endl;
0067     std::cout << "  energyVector.size()=" << energyVector.size() << std::endl;
0068     std::cout << "  eventMatrix[0].size()=" << eventMatrix[0].size() << std::endl;
0069     std::cout << " ******************    ERROR   *********************" << std::endl;
0070     return solution;  // empty vector
0071   }
0072
0073   // Reserve workspace
0074   float e25p;
0075   unsigned int i, j;
0076   std::vector<std::vector<float> > A(eventMatrix);
0077   std::vector<float> energies(energyVector);
0078
0079   float normalisation = 0.;
0080
0081   // Normalise if normaliseFlag is set
0082   if (normaliseFlag) {
0083     std::cout << "Householder::iterate(): Normalising event data" << std::endl;
0084     std::cout << "  WARNING: assuming 5x5 filtering has already been done" << std::endl;
0085
0086     for (i = 0; i < m; i++) {
0087       e25p = 0.;
0088       for (j = 0; j < n; j++) {
0089         e25p += eventMatrix[i][j];  // lorenzo -> trying to use ESetup which already performs calibs on rechits
0090       }
0091       e25p /= energyVector[i];
0092       normalisation += e25p;  // SUM e25p for all events
0093     }
0094     normalisation /= m;
0095     std::cout << "  Normalisation = " << normalisation << std::endl;
0096
0097     for (i = 0; i < energies.size(); ++i)
0098       energies[i] *= normalisation;
0099   }
0100
0101   // This is where the work goes on...
0102   // matrix decomposition
0103   std::vector<std::vector<float> > Acopy(A);
0104   std::vector<float> alpha(n);
0105   std::vector<int> pivot(n);
0106   if (!decompose(m, n, A, alpha, pivot)) {
0107     std::cout << "Householder::runIter(): Failed: Singular condition in decomposition." << std::endl;
0108     std::cout << "***************** PROBLEM in DECOMPOSITION *************************" << std::endl;
0109     return solution;  // empty vector
0110   }
0111
0112   /* DBL_EPSILON: Difference between 1.0 and the minimum float greater than 1.0 */
0113   float etasqr = DBL_EPSILON * DBL_EPSILON;
0114   std::cout << "LOOK at DBL_EPSILON :" << DBL_EPSILON << std::endl;
0115
0116   std::vector<float> r(energies);  // copy energies vector
0117   std::vector<float> e(n);
0118
0119   // apply transformations to rhs - find solution vector
0120   solution.assign(n, 0.);
0121   solve(m, n, A, alpha, pivot, r, solution);
0122
0123   // compute residual vector r
0124   for (i = 0; i < m; i++) {
0125     r[i] = energies[i];
0126     for (j = 0; j < n; j++)
0127       r[i] -= Acopy[i][j] * solution[j];
0128   }
0129   // compute first correction vector e
0130   solve(m, n, A, alpha, pivot, r, e);
0131
0132   float normy0 = 0.;
0133   float norme1 = 0.;
0134   float norme0;
0135
0136   for (i = 0; i < n; i++) {
0137     normy0 += solution[i] * solution[i];
0138     norme1 += e[i] * e[i];
0139   }
0140
0141   std::cout << "Householder::runIter(): applying first correction" << std::endl;
0142   std::cout << " normy0 = " << normy0 << std::endl;
0143   std::cout << " norme1 = " << norme1 << std::endl;
0144
0145   // not attempt at obtaining the solution is made unless the norm of the first
0146   // correction  is significantly smaller than the norm of the initial solution
0147   if (norme1 > (0.0625 * normy0)) {
0148     std::cout << "Householder::runIter(): first correction is too large. Failed." << std::endl;
0149   }
0150
0151   // improve the solution
0152   for (i = 0; i < n; i++)
0153     solution[i] += e[i];
0154
0155   std::cout << "Householder::runIter(): improving solution...." << std::endl;
0156
0157   // only continue iteration if the correction was significant
0158   while (norme1 > (etasqr * normy0)) {
0159     std::cout << "Householder::runIter(): norme1 = " << norme1 << std::endl;
0160
0161     for (i = 0; i < m; i++) {
0162       r[i] = energies[i];
0163       for (j = 0; j < n; j++)
0164         r[i] -= Acopy[i][j] * solution[j];
0165     }
0166
0167     // compute next correction vector
0168     solve(m, n, A, alpha, pivot, r, e);
0169
0170     norme0 = norme1;
0171     norme1 = 0.;
0172     for (i = 0; i < n; i++)
0173       norme1 += e[i] * e[i];
0174
0175     // terminate iteration if the norm of the new correction failed to decrease
0176     // significantly compared to the norm of the previous correction
0177     if (norme1 > (0.0625 * norme0))
0178       break;
0179
0180     // apply correction vector
0181     for (i = 0; i < n; i++)
0182       solution[i] += e[i];
0183   }
0184
0185   return solution;
0186 }
0187
0188 bool GenericHouseholder::decompose(const int m,
0189                                    const int n,
0190                                    std::vector<std::vector<float> >& qr,
0191                                    std::vector<float>& alpha,
0192                                    std::vector<int>& pivot) {
0193   int i, j, jbar, k;
0194   float beta, sigma, alphak, qrkk;
0195   std::vector<float> y(n);
0196   std::vector<float> sum(n);
0197
0198   std::cout << "Householder::decompose() started" << std::endl;
0199
0200   for (j = 0; j < n; j++) {
0201     // jth column sum
0202
0203     sum[j] = 0.;
0204     for (i = 0; i < m; i++)
0205       //      std::cout << "0: qr[i][j]" << qr[i][j] << " i = " << i << " j = " << j << std::endl;
0206       sum[j] += qr[i][j] * qr[i][j];
0207
0208     pivot[j] = j;
0209   }
0210
0211   for (k = 0; k < n; k++) {
0212     // kth Householder transformation
0213
0214     sigma = sum[k];
0215     jbar = k;
0216
0217     for (j = k + 1; j < n; j++) {
0218       if (sigma < sum[j]) {
0219         sigma = sum[j];
0220         jbar = j;
0221       }
0222     }
0223
0224     if (jbar != k) {
0225       // column interchange
0226       i = pivot[k];
0227       pivot[k] = pivot[jbar];
0228       pivot[jbar] = i;
0229       sum[jbar] = sum[k];
0230       sum[k] = sigma;
0231
0232       for (i = 0; i < m; i++) {
0233         sigma = qr[i][k];
0234         qr[i][k] = qr[i][jbar];
0235         //      std::cout << "A: qr[i][k]" << qr[i][k] << " i = " << i << " k = " << k << std::endl;
0236         qr[i][jbar] = sigma;
0237         //      std::cout << "B: qr[i][jbar]" << qr[i][k] << " i = " << i << " jbar = " << jbar << std::endl;
0238       }
0239     }  // end column interchange
0240
0241     sigma = 0.;
0242     for (i = k; i < m; i++) {
0243       sigma += qr[i][k] * qr[i][k];
0244       //      std::cout << "C: qr[i][k]" << qr[i][k] << " i = " << i << " k = " << k << std::endl;
0245     }
0246
0247     if (sigma == 0.) {
0248       std::cout << "Householder::decompose() failed" << std::endl;
0249       return false;
0250     }
0251
0252     qrkk = qr[k][k];
0253
0254     if (qrkk < 0.)
0255       alpha[k] = sqrt(sigma);
0256     else
0257       alpha[k] = sqrt(sigma) * (-1.);
0258     alphak = alpha[k];
0259
0260     beta = 1 / (sigma - qrkk * alphak);
0261     qr[k][k] = qrkk - alphak;
0262
0263     for (j = k + 1; j < n; j++) {
0264       y[j] = 0.;
0265       for (i = k; i < m; i++)
0266         y[j] += qr[i][k] * qr[i][j];
0267       y[j] *= beta;
0268     }
0269
0270     for (j = k + 1; j < n; j++) {
0271       for (i = k; i < m; i++) {
0272         qr[i][j] -= qr[i][k] * y[j];
0273         sum[j] -= qr[k][j] * qr[k][j];
0274       }
0275     }
0276   }  // end of kth householder transformation
0277
0278   std::cout << "Householder::decompose() finished" << std::endl;
0279
0280   return true;
0281 }
0282
0283 void GenericHouseholder::solve(int m,
0284                                int n,
0285                                const std::vector<std::vector<float> >& qr,
0286                                const std::vector<float>& alpha,
0287                                const std::vector<int>& pivot,
0288                                std::vector<float>& r,
0289                                std::vector<float>& y) {
0290   std::vector<float> z(n, 0.);
0291
0292   float gamma;
0293   int i, j;
0294
0295   std::cout << "Householder::solve() begin" << std::endl;
0296
0297   for (j = 0; j < n; j++) {
0298     // apply jth transformation to the right hand side
0299     gamma = 0.;
0300     for (i = j; i < m; i++)
0301       gamma += qr[i][j] * r[i];
0302     gamma /= (alpha[j] * qr[j][j]);
0303
0304     for (i = j; i < m; i++)
0305       r[i] += gamma * qr[i][j];
0306   }
0307
0308   //  std::cout<<"OK1:"<<std::endl;
0309   z[n - 1] = r[n - 1] / alpha[n - 1];
0310   //  std::cout<<"OK2:"<<std::endl;
0311
0312   for (i = n - 2; i >= 0; i--) {
0313     z[i] = r[i];
0314     for (j = i + 1; j < n; j++)
0315       z[i] -= qr[i][j] * z[j];
0316     z[i] /= alpha[i];
0317   }
0318   //  std::cout<<"OK3:"<<std::endl;
0319
0320   for (i = 0; i < n; i++)
0321     y[pivot[i]] = z[i];
0322
0323   std::cout << "Householder::solve() finished." << std::endl;
0324 }