PulseChiSqSNNLS Class Reference

#include <PulseChiSqSNNLS.h>

Public Types
typedef BXVector::Index	Index

Public Member Functions
const BXVector &	BXs () const
double	ChiSq () const
void	disableErrorCalculation ()
bool	DoFit (const SampleVector &samples, const SampleMatrix &samplecor, double pederr, const BXVector &bxs, const FullSampleVector &fullpulse, const FullSampleMatrix &fullpulsecov)
const PulseVector &	Errors () const
const SampleMatrix &	invcov () const
	PulseChiSqSNNLS ()
const SamplePulseMatrix &	pulsemat () const
const PulseVector &	X () const
	~PulseChiSqSNNLS ()

Protected Member Functions
double	ComputeApproxUncertainty (unsigned int ipulse)
double	ComputeChiSq ()
bool	Minimize (const SampleMatrix &samplecor, double pederr, const FullSampleMatrix &fullpulsecov)
bool	NNLS ()
bool	updateCov (const SampleMatrix &samplecor, double pederr, const FullSampleMatrix &fullpulsecov)

Protected Attributes
PulseVector	_ampvec
PulseVector	_ampvecmin
BXVector	_bxs
BXVector	_bxsmin
double	_chisq
bool	_computeErrors
SampleDecompLLT	_covdecomp
SampleMatrix	_covdecompLinv
PulseVector	_errvec
SampleMatrix	_invcov
unsigned int	_nP
unsigned int	_npulsetot
PulseDecompLDLT	_pulsedecomp
SamplePulseMatrix	_pulsemat
SampleVector	_sampvec
PulseMatrix	_topleft_work
PulseVector	ampvecpermtest
PulseMatrix	aTamat
PulseVector	aTbvec
SamplePulseMatrix	invcovp
PulseVector	updatework
PulseVector	wvec

Detailed Description

Definition at line 10 of file PulseChiSqSNNLS.h.

Member Typedef Documentation

typedef BXVector::Index PulseChiSqSNNLS::Index

Definition at line 13 of file PulseChiSqSNNLS.h.

Constructor & Destructor Documentation

PulseChiSqSNNLS::PulseChiSqSNNLS

(

)

Definition at line 75 of file PulseChiSqSNNLS.cc.

                                 :
  _chisq(0.),
  _computeErrors(true)
{
  
  Eigen::initParallel();
      
}  

PulseChiSqSNNLS::~PulseChiSqSNNLS

(

)

Definition at line 84 of file PulseChiSqSNNLS.cc.

                                  {
  
}

Member Function Documentation

const BXVector& PulseChiSqSNNLS::BXs ( void  ) const

inline

Definition at line 26 of file PulseChiSqSNNLS.h.

References _bxsmin.

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

{ return _bxsmin; }

double PulseChiSqSNNLS::ChiSq ( ) const

inline

Definition at line 28 of file PulseChiSqSNNLS.h.

References _chisq.

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

{ return _chisq; }

double PulseChiSqSNNLS::ComputeApproxUncertainty ( unsigned int  ipulse )

protected

Definition at line 268 of file PulseChiSqSNNLS.cc.

References _covdecomp, and _pulsemat.

Referenced by DoFit().

                                                                    {
  //compute approximate uncertainties
  //(using 1/second derivative since full Hessian is not meaningful in
  //presence of positive amplitude boundaries.)
      
  return 1./_covdecomp.matrixL().solve(_pulsemat.col(ipulse)).norm();
  
}

double PulseChiSqSNNLS::ComputeChiSq ( )

protected

Definition at line 259 of file PulseChiSqSNNLS.cc.

References _ampvec, _covdecomp, _pulsemat, and _sampvec.

Referenced by DoFit(), and Minimize().

                                     {
  
//   SampleVector resvec = _pulsemat*_ampvec - _sampvec;
//   return resvec.transpose()*_covdecomp.solve(resvec);
  
  return _covdecomp.matrixL().solve(_pulsemat*_ampvec - _sampvec).squaredNorm();
  
}

void PulseChiSqSNNLS::disableErrorCalculation ( )

inline

Definition at line 29 of file PulseChiSqSNNLS.h.

References _computeErrors.

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

{ _computeErrors = false; }

bool PulseChiSqSNNLS::DoFit	(	const SampleVector &	samples,
		const SampleMatrix &	samplecor,
		double	pederr,
		const BXVector &	bxs,
		const FullSampleVector &	fullpulse,
		const FullSampleMatrix &	fullpulsecov
	)

Definition at line 88 of file PulseChiSqSNNLS.cc.

References _ampvec, _ampvecmin, _bxs, _bxsmin, _chisq, _computeErrors, _errvec, _nP, _npulsetot, _pulsemat, _sampvec, funct::abs(), aTamat, ComputeApproxUncertainty(), ComputeChiSq(), max(), Minimize(), evf::evtn::offset(), mathSSE::sqrt(), ntuplemaker::status, std::swap(), and wvec.

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

                                                                                                                                                                                                   {
 
  //const unsigned int nsample = SampleVector::RowsAtCompileTime;
  _npulsetot = bxs.rows();
  //const unsigned int npulse = bxs.rows();

  _sampvec = samples;
  _bxs = bxs;
  
  //_pulsemat = SamplePulseMatrix::Zero(nsample,npulse);
  _pulsemat.resize(Eigen::NoChange,_npulsetot);
  _ampvec = PulseVector::Zero(_npulsetot);
  _errvec = PulseVector::Zero(_npulsetot);  
  _nP = 0;
  _chisq = 0.;
  
  aTamat.resize(_npulsetot,_npulsetot);
  wvec.resize(_npulsetot);

  //initialize pulse template matrix
  for (unsigned int ipulse=0; ipulse<_npulsetot; ++ipulse) {
    int bx = _bxs.coeff(ipulse);
    //int firstsamplet = std::max(0,bx + 3);
    int offset = 7-3-bx;
    
    //const unsigned int nsamplepulse = nsample-firstsamplet;
    //_pulsemat.col(ipulse).segment(firstsamplet,nsamplepulse) = fullpulse.segment(firstsamplet+offset,nsamplepulse);
    
    _pulsemat.col(ipulse) = fullpulse.segment<SampleVector::RowsAtCompileTime>(offset);
  }

  //do the actual fit
  bool status = Minimize(samplecor,pederr,fullpulsecov);
  _ampvecmin = _ampvec;
  _bxsmin = _bxs;
  
  if (!status) return status;
  
  if(!_computeErrors) return status;

  //compute MINOS-like uncertainties for in-time amplitude
  bool foundintime = false;
  unsigned int ipulseintime = 0;
  for (unsigned int ipulse=0; ipulse<_npulsetot; ++ipulse) {
    if (_bxs.coeff(ipulse)==0) {
      ipulseintime = ipulse;
      foundintime = true;
      break;
    }
  }
  if (!foundintime) return status;
  
  const unsigned int ipulseintimemin = ipulseintime;
  
  double approxerr = ComputeApproxUncertainty(ipulseintime);
  double chisq0 = _chisq;
  double x0 = _ampvecmin[ipulseintime];  
  
  //move in time pulse first to active set if necessary
  if (ipulseintime<_nP) {
    _pulsemat.col(_nP-1).swap(_pulsemat.col(ipulseintime));
    std::swap(_ampvec.coeffRef(_nP-1),_ampvec.coeffRef(ipulseintime));
    std::swap(_bxs.coeffRef(_nP-1),_bxs.coeffRef(ipulseintime));
    ipulseintime = _nP - 1;
    --_nP;    
  }
  
  
  SampleVector pulseintime = _pulsemat.col(ipulseintime);
  _pulsemat.col(ipulseintime).setZero();
  
  //two point interpolation for upper uncertainty when amplitude is away from boundary
  double xplus100 = x0 + approxerr;
  _ampvec.coeffRef(ipulseintime) = xplus100;
  _sampvec = samples - _ampvec.coeff(ipulseintime)*pulseintime;  
  status &= Minimize(samplecor,pederr,fullpulsecov);
  if (!status) return status;
  double chisqplus100 = ComputeChiSq();
  
  double sigmaplus = std::abs(xplus100-x0)/sqrt(chisqplus100-chisq0);
  
  //if amplitude is sufficiently far from the boundary, compute also the lower uncertainty and average them
  if ( (x0/sigmaplus) > 0.5 ) {
    for (unsigned int ipulse=0; ipulse<_npulsetot; ++ipulse) {
      if (_bxs.coeff(ipulse)==0) {
        ipulseintime = ipulse;
        break;
      }
    }    
    double xminus100 = std::max(0.,x0-approxerr);
    _ampvec.coeffRef(ipulseintime) = xminus100;
   _sampvec = samples - _ampvec.coeff(ipulseintime)*pulseintime;
    status &= Minimize(samplecor,pederr,fullpulsecov);
    if (!status) return status;
    double chisqminus100 = ComputeChiSq();
    
    double sigmaminus = std::abs(xminus100-x0)/sqrt(chisqminus100-chisq0);
    _errvec[ipulseintimemin] = 0.5*(sigmaplus + sigmaminus);
    
  }
  else {
    _errvec[ipulseintimemin] = sigmaplus;
  }
            
  _chisq = chisq0;  

  return status;
  
}

const PulseVector& PulseChiSqSNNLS::Errors ( ) const

inline

Definition at line 25 of file PulseChiSqSNNLS.h.

References _errvec.

Referenced by EcalUncalibRecHitMultiFitAlgo::makeRecHit().

{ return _errvec; }

const SampleMatrix& PulseChiSqSNNLS::invcov ( ) const

inline

Definition at line 22 of file PulseChiSqSNNLS.h.

References _invcov.

{ return _invcov; }

bool PulseChiSqSNNLS::Minimize ( const SampleMatrix &  samplecor, double  pederr, const FullSampleMatrix &  fullpulsecov  )

protected

Definition at line 198 of file PulseChiSqSNNLS.cc.

References _chisq, funct::abs(), ComputeChiSq(), alignCSCRings::e, getDQMSummary::iter, maxiter, NNLS(), ntuplemaker::status, and updateCov().

Referenced by DoFit().

                                                                                                                 {

  
  const int maxiter = 50;
  int iter = 0;
  bool status = false;
  while (true) {    
    
    if (iter>=maxiter) {
      edm::LogWarning("PulseChiSqSNNLS::Minimize") << "Max Iterations reached at iter " << iter <<  std::endl;
      break;
    }    
    
    status = updateCov(samplecor,pederr,fullpulsecov);    
    if (!status) break;    
    status = NNLS();
    if (!status) break;
        
    double chisqnow = ComputeChiSq();
    double deltachisq = chisqnow-_chisq;
        
    _chisq = chisqnow;
    if (std::abs(deltachisq)<1e-3) {
      break;
    }
    ++iter;    
  }  
  
  return status;  
  
}

bool PulseChiSqSNNLS::NNLS ( )

protected

Definition at line 277 of file PulseChiSqSNNLS.cc.

References _ampvec, _bxs, _covdecomp, _nP, _pulsemat, _sampvec, ampvecpermtest, aTamat, aTbvec, alignCSCRings::e, eigen_solve_submatrix(), invcovp, getDQMSummary::iter, max(), std::swap(), updatework, and wvec.

Referenced by Minimize().

                           {
  
  //Fast NNLS (fnnls) algorithm as per http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.157.9203&rep=rep1&type=pdf
  
  const unsigned int npulse = _bxs.rows();

  invcovp = _covdecomp.matrixL().solve(_pulsemat);
  aTamat = invcovp.transpose()*invcovp; //.triangularView<Eigen::Lower>()
  //aTamat = aTamat.selfadjointView<Eigen::Lower>();  
  aTbvec = invcovp.transpose()*_covdecomp.matrixL().solve(_sampvec);  
  
  int iter = 0;
  Index idxwmax;
  double wmax = 0.0;
  //work = PulseVector::zeros();
  while (true) {    
    //can only perform this step if solution is guaranteed viable
    if (iter>0 || _nP==0) {
      if ( _nP==npulse ) break;                  
      
      const unsigned int nActive = npulse - _nP;
      
      updatework = aTbvec - aTamat*_ampvec;      
      wmax = updatework.tail(nActive).maxCoeff(&idxwmax);
      
      //convergence
      if (wmax<1e-11) break;
      
      //unconstrain parameter
      Index idxp = _nP + idxwmax;
      //printf("adding index %i, orig index %i\n",int(idxp),int(_bxs.coeff(idxp)));
      aTamat.col(_nP).swap(aTamat.col(idxp));
      aTamat.row(_nP).swap(aTamat.row(idxp));
      _pulsemat.col(_nP).swap(_pulsemat.col(idxp));
      std::swap(aTbvec.coeffRef(_nP),aTbvec.coeffRef(idxp));
      std::swap(_ampvec.coeffRef(_nP),_ampvec.coeffRef(idxp));
      std::swap(_bxs.coeffRef(_nP),_bxs.coeffRef(idxp));
      
      // update now that we are done doing work
      wvec.tail(nActive) = updatework.tail(nActive); 
      ++_nP;
    }

    
    while (true) {
      //printf("iter in, idxsP = %i\n",int(_idxsP.size()));
      
      if (_nP==0) break;     
      
      ampvecpermtest = _ampvec;
      
      //solve for unconstrained parameters 
      //need to have specialized function to call optimized versions
      // of matrix solver... this is truly amazing...
      eigen_solve_submatrix(aTamat,aTbvec,ampvecpermtest,_nP);
      
      //check solution
      auto ampvecpermhead = ampvecpermtest.head(_nP);
      if ( ampvecpermhead.minCoeff()>0. ) {
        _ampvec.head(_nP) = ampvecpermhead.head(_nP);
        break;
      }      

      //update parameter vector
      Index minratioidx=0;
      
      // no realizable optimization here (because it autovectorizes!)
      double minratio = std::numeric_limits<double>::max();
      for (unsigned int ipulse=0; ipulse<_nP; ++ipulse) {
        if (ampvecpermtest.coeff(ipulse)<=0.) {
      const double c_ampvec = _ampvec.coeff(ipulse);
          const double ratio = c_ampvec/(c_ampvec-ampvecpermtest.coeff(ipulse));
          if (ratio<minratio) {
            minratio = ratio;
            minratioidx = ipulse;
          }
        }
      }

      _ampvec.head(_nP) += minratio*(ampvecpermhead - _ampvec.head(_nP));
      
      //avoid numerical problems with later ==0. check
      _ampvec.coeffRef(minratioidx) = 0.;
            
      //printf("removing index %i, orig idx %i\n",int(minratioidx),int(_bxs.coeff(minratioidx)));
      aTamat.col(_nP-1).swap(aTamat.col(minratioidx));
      aTamat.row(_nP-1).swap(aTamat.row(minratioidx));
      _pulsemat.col(_nP-1).swap(_pulsemat.col(minratioidx));
      std::swap(aTbvec.coeffRef(_nP-1),aTbvec.coeffRef(minratioidx));
      std::swap(_ampvec.coeffRef(_nP-1),_ampvec.coeffRef(minratioidx));
      std::swap(_bxs.coeffRef(_nP-1),_bxs.coeffRef(minratioidx));
      --_nP;      
    }
    ++iter;
  }
  
  return true;
  
  
}

const SamplePulseMatrix& PulseChiSqSNNLS::pulsemat ( ) const

inline

Definition at line 21 of file PulseChiSqSNNLS.h.

References _pulsemat.

{ return _pulsemat; }

bool PulseChiSqSNNLS::updateCov ( const SampleMatrix &  samplecor, double  pederr, const FullSampleMatrix &  fullpulsecov  )

protected

Definition at line 230 of file PulseChiSqSNNLS.cc.

References _ampvec, _bxs, _covdecomp, _invcov, max(), evf::evtn::offset(), and ntuplemaker::status.

Referenced by Minimize().

                                                                                                                  {
 
  const unsigned int nsample = SampleVector::RowsAtCompileTime;
  const unsigned int npulse = _bxs.rows();

  const double pederr2 = pederr*pederr;
  _invcov = pederr2*samplecor; //
  
  for (unsigned int ipulse=0; ipulse<npulse; ++ipulse) {
    if (_ampvec.coeff(ipulse)==0.) continue;
    int bx = _bxs.coeff(ipulse);
    int firstsamplet = std::max(0,bx + 3);
    int offset = 7-3-bx;
    
    const double ampveccoef = _ampvec.coeff(ipulse);
    const double ampsq = ampveccoef*ampveccoef;
    
    const unsigned int nsamplepulse = nsample-firstsamplet;    
    _invcov.block(firstsamplet,firstsamplet,nsamplepulse,nsamplepulse) += 
      ampsq*fullpulsecov.block(firstsamplet+offset,firstsamplet+offset,nsamplepulse,nsamplepulse);   
  }
  
  _covdecomp.compute(_invcov);
  
  bool status = true;
  return status;
    
}

const PulseVector& PulseChiSqSNNLS::X ( ) const

inline

Definition at line 24 of file PulseChiSqSNNLS.h.

References _ampvecmin.

Referenced by svgfig.Curve.Sample::__repr__(), and EcalUncalibRecHitMultiFitAlgo::makeRecHit().

{ return _ampvecmin; }

Member Data Documentation

PulseVector PulseChiSqSNNLS::_ampvec

protected

Definition at line 43 of file PulseChiSqSNNLS.h.

Referenced by ComputeChiSq(), DoFit(), NNLS(), and updateCov().

PulseVector PulseChiSqSNNLS::_ampvecmin

protected

Definition at line 45 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), and X().

BXVector PulseChiSqSNNLS::_bxs

protected

Definition at line 52 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), NNLS(), and updateCov().

BXVector PulseChiSqSNNLS::_bxsmin

protected

Definition at line 53 of file PulseChiSqSNNLS.h.

Referenced by BXs(), and DoFit().

double PulseChiSqSNNLS::_chisq

protected

Definition at line 65 of file PulseChiSqSNNLS.h.

Referenced by ChiSq(), DoFit(), and Minimize().

bool PulseChiSqSNNLS::_computeErrors

protected

Definition at line 66 of file PulseChiSqSNNLS.h.

Referenced by disableErrorCalculation(), and DoFit().

SampleDecompLLT PulseChiSqSNNLS::_covdecomp

protected

Definition at line 47 of file PulseChiSqSNNLS.h.

Referenced by ComputeApproxUncertainty(), ComputeChiSq(), NNLS(), and updateCov().

SampleMatrix PulseChiSqSNNLS::_covdecompLinv

protected

Definition at line 48 of file PulseChiSqSNNLS.h.

PulseVector PulseChiSqSNNLS::_errvec

protected

Definition at line 44 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), and Errors().

SampleMatrix PulseChiSqSNNLS::_invcov

protected

Definition at line 41 of file PulseChiSqSNNLS.h.

Referenced by invcov(), and updateCov().

unsigned int PulseChiSqSNNLS::_nP

protected

Definition at line 55 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), and NNLS().

unsigned int PulseChiSqSNNLS::_npulsetot

protected

Definition at line 54 of file PulseChiSqSNNLS.h.

Referenced by DoFit().

PulseDecompLDLT PulseChiSqSNNLS::_pulsedecomp

protected

Definition at line 50 of file PulseChiSqSNNLS.h.

SamplePulseMatrix PulseChiSqSNNLS::_pulsemat

protected

Definition at line 42 of file PulseChiSqSNNLS.h.

Referenced by ComputeApproxUncertainty(), ComputeChiSq(), DoFit(), NNLS(), and pulsemat().

SampleVector PulseChiSqSNNLS::_sampvec

protected

Definition at line 40 of file PulseChiSqSNNLS.h.

Referenced by ComputeChiSq(), DoFit(), and NNLS().

PulseMatrix PulseChiSqSNNLS::_topleft_work

protected

Definition at line 49 of file PulseChiSqSNNLS.h.

PulseVector PulseChiSqSNNLS::ampvecpermtest

protected

Definition at line 63 of file PulseChiSqSNNLS.h.

Referenced by NNLS().

PulseMatrix PulseChiSqSNNLS::aTamat

protected

Definition at line 58 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), and NNLS().

PulseVector PulseChiSqSNNLS::aTbvec

protected

Definition at line 59 of file PulseChiSqSNNLS.h.

Referenced by NNLS().

SamplePulseMatrix PulseChiSqSNNLS::invcovp

protected

Definition at line 57 of file PulseChiSqSNNLS.h.

Referenced by NNLS().

PulseVector PulseChiSqSNNLS::updatework

protected

Definition at line 61 of file PulseChiSqSNNLS.h.

Referenced by NNLS().

PulseVector PulseChiSqSNNLS::wvec

protected

Definition at line 60 of file PulseChiSqSNNLS.h.

Referenced by DoFit(), and NNLS().