EcalRecAlgos/src/PulseChiSqSNNLS.cc

0001 #include "RecoLocalCalo/EcalRecAlgos/interface/PulseChiSqSNNLS.h"
0002 #include <cmath>
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include <iostream>
0005
0006 void eigen_solve_submatrix(PulseMatrix &mat, PulseVector &invec, PulseVector &outvec, unsigned NP) {
0007   using namespace Eigen;
0008   switch (NP) {  // pulse matrix is always square.
0009     case 10: {
0010       Matrix<double, 10, 10> temp = mat.topLeftCorner<10, 10>();
0011       outvec.head<10>() = temp.ldlt().solve(invec.head<10>());
0012     } break;
0013     case 9: {
0014       Matrix<double, 9, 9> temp = mat.topLeftCorner<9, 9>();
0015       outvec.head<9>() = temp.ldlt().solve(invec.head<9>());
0016     } break;
0017     case 8: {
0018       Matrix<double, 8, 8> temp = mat.topLeftCorner<8, 8>();
0019       outvec.head<8>() = temp.ldlt().solve(invec.head<8>());
0020     } break;
0021     case 7: {
0022       Matrix<double, 7, 7> temp = mat.topLeftCorner<7, 7>();
0023       outvec.head<7>() = temp.ldlt().solve(invec.head<7>());
0024     } break;
0025     case 6: {
0026       Matrix<double, 6, 6> temp = mat.topLeftCorner<6, 6>();
0027       outvec.head<6>() = temp.ldlt().solve(invec.head<6>());
0028     } break;
0029     case 5: {
0030       Matrix<double, 5, 5> temp = mat.topLeftCorner<5, 5>();
0031       outvec.head<5>() = temp.ldlt().solve(invec.head<5>());
0032     } break;
0033     case 4: {
0034       Matrix<double, 4, 4> temp = mat.topLeftCorner<4, 4>();
0035       outvec.head<4>() = temp.ldlt().solve(invec.head<4>());
0036     } break;
0037     case 3: {
0038       Matrix<double, 3, 3> temp = mat.topLeftCorner<3, 3>();
0039       outvec.head<3>() = temp.ldlt().solve(invec.head<3>());
0040     } break;
0041     case 2: {
0042       Matrix<double, 2, 2> temp = mat.topLeftCorner<2, 2>();
0043       outvec.head<2>() = temp.ldlt().solve(invec.head<2>());
0044     } break;
0045     case 1: {
0046       Matrix<double, 1, 1> temp = mat.topLeftCorner<1, 1>();
0047       outvec.head<1>() = temp.ldlt().solve(invec.head<1>());
0048     } break;
0049     default:
0050       throw cms::Exception("MultFitWeirdState")
0051           << "Weird number of pulses encountered in multifit, module is configured incorrectly!";
0052   }
0053 }
0054
0055 PulseChiSqSNNLS::PulseChiSqSNNLS() : _chisq(0.), _computeErrors(true), _maxiters(50), _maxiterwarnings(true) {}
0056
0057 PulseChiSqSNNLS::~PulseChiSqSNNLS() {}
0058
0059 bool PulseChiSqSNNLS::DoFit(const SampleVector &samples,
0060                             const SampleMatrix &samplecov,
0061                             const BXVector &bxs,
0062                             const FullSampleVector &fullpulse,
0063                             const FullSampleMatrix &fullpulsecov,
0064                             const SampleGainVector &gains,
0065                             const SampleGainVector &badSamples) {
0066   int npulse = bxs.rows();
0067
0068   _sampvec = samples;
0069   _bxs = bxs;
0070   _pulsemat.resize(Eigen::NoChange, npulse);
0071
0072   //construct dynamic pedestals if applicable
0073   int ngains = gains.maxCoeff() + 1;
0074   int nPedestals = 0;
0075   for (int gainidx = 0; gainidx < ngains; ++gainidx) {
0076     SampleGainVector mask = gainidx * SampleGainVector::Ones();
0077     SampleVector pedestal = (gains.array() == mask.array()).cast<SampleVector::value_type>();
0078     if (pedestal.maxCoeff() > 0.) {
0079       ++nPedestals;
0080       _bxs.resize(npulse + nPedestals);
0081       _bxs[npulse + nPedestals - 1] = 100 + gainidx;  //bx values >=100 indicate dynamic pedestals
0082       _pulsemat.resize(Eigen::NoChange, npulse + nPedestals);
0083       _pulsemat.col(npulse + nPedestals - 1) = pedestal;
0084     }
0085   }
0086
0087   //construct negative step functions for saturated or potentially slew-rate-limited samples
0088   for (int isample = 0; isample < SampleVector::RowsAtCompileTime; ++isample) {
0089     if (badSamples.coeff(isample) > 0) {
0090       SampleVector step = SampleVector::Zero();
0091       //step correction has negative sign for saturated or slew-limited samples which have been forced to zero
0092       step[isample] = -1.;
0093
0094       ++nPedestals;
0095       _bxs.resize(npulse + nPedestals);
0096       _bxs[npulse + nPedestals - 1] =
0097           -100 - isample;  //bx values <=-100 indicate step corrections for saturated or slew-limited samples
0098       _pulsemat.resize(Eigen::NoChange, npulse + nPedestals);
0099       _pulsemat.col(npulse + nPedestals - 1) = step;
0100     }
0101   }
0102
0103   _npulsetot = npulse + nPedestals;
0104
0105   _ampvec = PulseVector::Zero(_npulsetot);
0106   _errvec = PulseVector::Zero(_npulsetot);
0107   _nP = 0;
0108   _chisq = 0.;
0109
0110   if (_bxs.rows() == 1 && std::abs(_bxs.coeff(0)) < 100) {
0111     _ampvec.coeffRef(0) = _sampvec.coeff(_bxs.coeff(0) + 5);
0112   }
0113
0114   aTamat.resize(_npulsetot, _npulsetot);
0115
0116   //initialize pulse template matrix
0117   for (int ipulse = 0; ipulse < npulse; ++ipulse) {
0118     int bx = _bxs.coeff(ipulse);
0119     int offset = 7 - 3 - bx;
0120     _pulsemat.col(ipulse) = fullpulse.segment<SampleVector::RowsAtCompileTime>(offset);
0121   }
0122
0123   //unconstrain pedestals already for first iteration since they should always be non-zero
0124   if (nPedestals > 0) {
0125     for (int i = 0; i < _bxs.rows(); ++i) {
0126       int bx = _bxs.coeff(i);
0127       if (bx >= 100) {
0128         NNLSUnconstrainParameter(i);
0129       }
0130     }
0131   }
0132
0133   //do the actual fit
0134   bool status = Minimize(samplecov, fullpulsecov);
0135   _ampvecmin = _ampvec;
0136   _bxsmin = _bxs;
0137
0138   if (!status)
0139     return status;
0140
0141   if (!_computeErrors)
0142     return status;
0143
0144   //compute MINOS-like uncertainties for in-time amplitude
0145   bool foundintime = false;
0146   unsigned int ipulseintime = 0;
0147   for (unsigned int ipulse = 0; ipulse < _npulsetot; ++ipulse) {
0148     if (_bxs.coeff(ipulse) == 0) {
0149       ipulseintime = ipulse;
0150       foundintime = true;
0151       break;
0152     }
0153   }
0154   if (!foundintime)
0155     return status;
0156
0157   const unsigned int ipulseintimemin = ipulseintime;
0158
0159   double approxerr = ComputeApproxUncertainty(ipulseintime);
0160   double chisq0 = _chisq;
0161   double x0 = _ampvecmin[ipulseintime];
0162
0163   //move in time pulse first to active set if necessary
0164   if (ipulseintime < _nP) {
0165     _pulsemat.col(_nP - 1).swap(_pulsemat.col(ipulseintime));
0166     std::swap(_ampvec.coeffRef(_nP - 1), _ampvec.coeffRef(ipulseintime));
0167     std::swap(_bxs.coeffRef(_nP - 1), _bxs.coeffRef(ipulseintime));
0168     ipulseintime = _nP - 1;
0169     --_nP;
0170   }
0171
0172   SampleVector pulseintime = _pulsemat.col(ipulseintime);
0173   _pulsemat.col(ipulseintime).setZero();
0174
0175   //two point interpolation for upper uncertainty when amplitude is away from boundary
0176   double xplus100 = x0 + approxerr;
0177   _ampvec.coeffRef(ipulseintime) = xplus100;
0178   _sampvec = samples - _ampvec.coeff(ipulseintime) * pulseintime;
0179   status &= Minimize(samplecov, fullpulsecov);
0180   if (!status)
0181     return status;
0182   double chisqplus100 = ComputeChiSq();
0183
0184   double sigmaplus = std::abs(xplus100 - x0) / sqrt(chisqplus100 - chisq0);
0185
0186   //if amplitude is sufficiently far from the boundary, compute also the lower uncertainty and average them
0187   if ((x0 / sigmaplus) > 0.5) {
0188     for (unsigned int ipulse = 0; ipulse < _npulsetot; ++ipulse) {
0189       if (_bxs.coeff(ipulse) == 0) {
0190         ipulseintime = ipulse;
0191         break;
0192       }
0193     }
0194     double xminus100 = std::max(0., x0 - approxerr);
0195     _ampvec.coeffRef(ipulseintime) = xminus100;
0196     _sampvec = samples - _ampvec.coeff(ipulseintime) * pulseintime;
0197     status &= Minimize(samplecov, fullpulsecov);
0198     if (!status)
0199       return status;
0200     double chisqminus100 = ComputeChiSq();
0201
0202     double sigmaminus = std::abs(xminus100 - x0) / sqrt(chisqminus100 - chisq0);
0203     _errvec[ipulseintimemin] = 0.5 * (sigmaplus + sigmaminus);
0204
0205   } else {
0206     _errvec[ipulseintimemin] = sigmaplus;
0207   }
0208
0209   _chisq = chisq0;
0210
0211   return status;
0212 }
0213
0214 bool PulseChiSqSNNLS::Minimize(const SampleMatrix &samplecov, const FullSampleMatrix &fullpulsecov) {
0215   const unsigned int npulse = _bxs.rows();
0216
0217   int iter = 0;
0218   bool status = false;
0219   while (true) {
0220     if (iter >= _maxiters) {
0221       if (_maxiterwarnings) {
0222         LogDebug("PulseChiSqSNNLS::Minimize") << "Max Iterations reached at iter " << iter;
0223       }
0224       break;
0225     }
0226
0227     status = updateCov(samplecov, fullpulsecov);
0228     if (!status)
0229       break;
0230     if (npulse > 1) {
0231       status = NNLS();
0232     } else {
0233       //special case for one pulse fit (performance optimized)
0234       status = OnePulseMinimize();
0235     }
0236     if (!status)
0237       break;
0238
0239     double chisqnow = ComputeChiSq();
0240     double deltachisq = chisqnow - _chisq;
0241
0242     _chisq = chisqnow;
0243     if (std::abs(deltachisq) < 1e-3) {
0244       break;
0245     }
0246     ++iter;
0247   }
0248
0249   return status;
0250 }
0251
0252 bool PulseChiSqSNNLS::updateCov(const SampleMatrix &samplecov, const FullSampleMatrix &fullpulsecov) {
0253   const unsigned int nsample = SampleVector::RowsAtCompileTime;
0254   const unsigned int npulse = _bxs.rows();
0255
0256   _invcov = samplecov;  //
0257
0258   for (unsigned int ipulse = 0; ipulse < npulse; ++ipulse) {
0259     if (_ampvec.coeff(ipulse) == 0.)
0260       continue;
0261     int bx = _bxs.coeff(ipulse);
0262     if (std::abs(bx) >= 100)
0263       continue;  //no contribution to covariance from pedestal or saturation/slew step correction
0264
0265     int firstsamplet = std::max(0, bx + 3);
0266     int offset = 7 - 3 - bx;
0267
0268     const double ampveccoef = _ampvec.coeff(ipulse);
0269     const double ampsq = ampveccoef * ampveccoef;
0270
0271     const unsigned int nsamplepulse = nsample - firstsamplet;
0272     _invcov.block(firstsamplet, firstsamplet, nsamplepulse, nsamplepulse) +=
0273         ampsq * fullpulsecov.block(firstsamplet + offset, firstsamplet + offset, nsamplepulse, nsamplepulse);
0274   }
0275
0276   _covdecomp.compute(_invcov);
0277
0278   bool status = true;
0279   return status;
0280 }
0281
0282 double PulseChiSqSNNLS::ComputeChiSq() {
0283   //   SampleVector resvec = _pulsemat*_ampvec - _sampvec;
0284   //   return resvec.transpose()*_covdecomp.solve(resvec);
0285
0286   return _covdecomp.matrixL().solve(_pulsemat * _ampvec - _sampvec).squaredNorm();
0287 }
0288
0289 double PulseChiSqSNNLS::ComputeApproxUncertainty(unsigned int ipulse) {
0290   //compute approximate uncertainties
0291   //(using 1/second derivative since full Hessian is not meaningful in
0292   //presence of positive amplitude boundaries.)
0293
0294   return 1. / _covdecomp.matrixL().solve(_pulsemat.col(ipulse)).norm();
0295 }
0296
0297 bool PulseChiSqSNNLS::NNLS() {
0298   //Fast NNLS (fnnls) algorithm as per http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.9203&rep=rep1&type=pdf
0299
0300   const unsigned int npulse = _bxs.rows();
0301   constexpr unsigned int nsamples = SampleVector::RowsAtCompileTime;
0302
0303   invcovp = _covdecomp.matrixL().solve(_pulsemat);
0304   aTamat.noalias() = invcovp.transpose().lazyProduct(invcovp);
0305   aTbvec.noalias() = invcovp.transpose().lazyProduct(_covdecomp.matrixL().solve(_sampvec));
0306
0307   int iter = 0;
0308   Index idxwmax = 0;
0309   double wmax = 0.0;
0310   double threshold = 1e-11;
0311   while (true) {
0312     //can only perform this step if solution is guaranteed viable
0313     if (iter > 0 || _nP == 0) {
0314       if (_nP == std::min(npulse, nsamples))
0315         break;
0316
0317       const unsigned int nActive = npulse - _nP;
0318
0319       updatework = aTbvec - aTamat * _ampvec;
0320       Index idxwmaxprev = idxwmax;
0321       double wmaxprev = wmax;
0322       wmax = updatework.tail(nActive).maxCoeff(&idxwmax);
0323
0324       //convergence
0325       if (wmax < threshold || (idxwmax == idxwmaxprev && wmax == wmaxprev))
0326         break;
0327
0328       //worst case protection
0329       if (iter >= 500) {
0330         LogDebug("PulseChiSqSNNLS::NNLS()") << "Max Iterations reached at iter " << iter;
0331         break;
0332       }
0333
0334       //unconstrain parameter
0335       Index idxp = _nP + idxwmax;
0336       NNLSUnconstrainParameter(idxp);
0337     }
0338
0339     while (true) {
0340       //printf("iter in, idxsP = %i\n",int(_idxsP.size()));
0341
0342       if (_nP == 0)
0343         break;
0344
0345       ampvecpermtest = _ampvec;
0346
0347       //solve for unconstrained parameters
0348       //need to have specialized function to call optimized versions
0349       // of matrix solver... this is truly amazing...
0350       eigen_solve_submatrix(aTamat, aTbvec, ampvecpermtest, _nP);
0351
0352       //check solution
0353       bool positive = true;
0354       for (unsigned int i = 0; i < _nP; ++i)
0355         positive &= (ampvecpermtest(i) > 0);
0356       if (positive) {
0357         _ampvec.head(_nP) = ampvecpermtest.head(_nP);
0358         break;
0359       }
0360
0361       //update parameter vector
0362       Index minratioidx = 0;
0363
0364       // no realizable optimization here (because it autovectorizes!)
0365       double minratio = std::numeric_limits<double>::max();
0366       for (unsigned int ipulse = 0; ipulse < _nP; ++ipulse) {
0367         if (ampvecpermtest.coeff(ipulse) <= 0.) {
0368           const double c_ampvec = _ampvec.coeff(ipulse);
0369           const double ratio = c_ampvec / (c_ampvec - ampvecpermtest.coeff(ipulse));
0370           if (ratio < minratio) {
0371             minratio = ratio;
0372             minratioidx = ipulse;
0373           }
0374         }
0375       }
0376
0377       _ampvec.head(_nP) += minratio * (ampvecpermtest.head(_nP) - _ampvec.head(_nP));
0378
0379       //avoid numerical problems with later ==0. check
0380       _ampvec.coeffRef(minratioidx) = 0.;
0381
0382       //printf("removing index %i, orig idx %i\n",int(minratioidx),int(_bxs.coeff(minratioidx)));
0383       NNLSConstrainParameter(minratioidx);
0384     }
0385     ++iter;
0386
0387     //adaptive convergence threshold to avoid infinite loops but still
0388     //ensure best value is used
0389     if (iter % 16 == 0) {
0390       threshold *= 2;
0391     }
0392   }
0393
0394   return true;
0395 }
0396
0397 void PulseChiSqSNNLS::NNLSUnconstrainParameter(Index idxp) {
0398   aTamat.col(_nP).swap(aTamat.col(idxp));
0399   aTamat.row(_nP).swap(aTamat.row(idxp));
0400   _pulsemat.col(_nP).swap(_pulsemat.col(idxp));
0401   std::swap(aTbvec.coeffRef(_nP), aTbvec.coeffRef(idxp));
0402   std::swap(_ampvec.coeffRef(_nP), _ampvec.coeffRef(idxp));
0403   std::swap(_bxs.coeffRef(_nP), _bxs.coeffRef(idxp));
0404   ++_nP;
0405 }
0406
0407 void PulseChiSqSNNLS::NNLSConstrainParameter(Index minratioidx) {
0408   aTamat.col(_nP - 1).swap(aTamat.col(minratioidx));
0409   aTamat.row(_nP - 1).swap(aTamat.row(minratioidx));
0410   _pulsemat.col(_nP - 1).swap(_pulsemat.col(minratioidx));
0411   std::swap(aTbvec.coeffRef(_nP - 1), aTbvec.coeffRef(minratioidx));
0412   std::swap(_ampvec.coeffRef(_nP - 1), _ampvec.coeffRef(minratioidx));
0413   std::swap(_bxs.coeffRef(_nP - 1), _bxs.coeffRef(minratioidx));
0414   --_nP;
0415 }
0416
0417 bool PulseChiSqSNNLS::OnePulseMinimize() {
0418   //Fast NNLS (fnnls) algorithm as per http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.9203&rep=rep1&type=pdf
0419
0420   //   const unsigned int npulse = 1;
0421
0422   invcovp = _covdecomp.matrixL().solve(_pulsemat);
0423   //   aTamat = invcovp.transpose()*invcovp;
0424   //   aTbvec = invcovp.transpose()*_covdecomp.matrixL().solve(_sampvec);
0425
0426   SingleMatrix aTamatval = invcovp.transpose() * invcovp;
0427   SingleVector aTbvecval = invcovp.transpose() * _covdecomp.matrixL().solve(_sampvec);
0428   _ampvec.coeffRef(0) = std::max(0., aTbvecval.coeff(0) / aTamatval.coeff(0));
0429
0430   return true;
0431 }