File indexing completed on 2023-03-17 11:18:42
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) {
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
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;
0082 _pulsemat.resize(Eigen::NoChange, npulse + nPedestals);
0083 _pulsemat.col(npulse + nPedestals - 1) = pedestal;
0084 }
0085 }
0086
0087
0088 for (int isample = 0; isample < SampleVector::RowsAtCompileTime; ++isample) {
0089 if (badSamples.coeff(isample) > 0) {
0090 SampleVector step = SampleVector::Zero();
0091
0092 step[isample] = -1.;
0093
0094 ++nPedestals;
0095 _bxs.resize(npulse + nPedestals);
0096 _bxs[npulse + nPedestals - 1] =
0097 -100 - isample;
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
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
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
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
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
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
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
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
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;
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
0284
0285
0286 return _covdecomp.matrixL().solve(_pulsemat * _ampvec - _sampvec).squaredNorm();
0287 }
0288
0289 double PulseChiSqSNNLS::ComputeApproxUncertainty(unsigned int ipulse) {
0290
0291
0292
0293
0294 return 1. / _covdecomp.matrixL().solve(_pulsemat.col(ipulse)).norm();
0295 }
0296
0297 bool PulseChiSqSNNLS::NNLS() {
0298
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
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
0325 if (wmax < threshold || (idxwmax == idxwmaxprev && wmax == wmaxprev))
0326 break;
0327
0328
0329 if (iter >= 500) {
0330 LogDebug("PulseChiSqSNNLS::NNLS()") << "Max Iterations reached at iter " << iter;
0331 break;
0332 }
0333
0334
0335 Index idxp = _nP + idxwmax;
0336 NNLSUnconstrainParameter(idxp);
0337 }
0338
0339 while (true) {
0340
0341
0342 if (_nP == 0)
0343 break;
0344
0345 ampvecpermtest = _ampvec;
0346
0347
0348
0349
0350 eigen_solve_submatrix(aTamat, aTbvec, ampvecpermtest, _nP);
0351
0352
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
0362 Index minratioidx = 0;
0363
0364
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
0380 _ampvec.coeffRef(minratioidx) = 0.;
0381
0382
0383 NNLSConstrainParameter(minratioidx);
0384 }
0385 ++iter;
0386
0387
0388
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
0419
0420
0421
0422 invcovp = _covdecomp.matrixL().solve(_pulsemat);
0423
0424
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 }