Page Menu
Home
HEPForge
Search
Configure Global Search
Log In
Files
F19250555
No One
Temporary
Actions
View File
Edit File
Delete File
View Transforms
Subscribe
Award Token
Flag For Later
Size
18 KB
Referenced Files
None
Subscribers
None
View Options
Index: trunk/npstat/stat/SmoothedDAgostiniUnfold1D.cc
===================================================================
--- trunk/npstat/stat/SmoothedDAgostiniUnfold1D.cc (revision 187)
+++ trunk/npstat/stat/SmoothedDAgostiniUnfold1D.cc (revision 188)
@@ -1,310 +1,315 @@
#include <cassert>
#include <utility>
#include <algorithm>
#include <sstream>
#include <numeric>
#include "npstat/nm/allocators.hh"
#include "npstat/stat/SmoothedDAgostiniUnfold1D.hh"
#include "npstat/stat/StatUtils.hh"
static double probDelta(const double* prev, const double* next,
const unsigned len)
{
long double del = 0.0L;
long double sum = 0.0L;
for (unsigned i=0; i<len; ++i)
{
del += fabs(*prev - *next);
sum += fabs(*prev++);
sum += fabs(*next++);
}
sum /= 2.0L;
double ratio = 0.0;
if (sum)
ratio = del/sum;
return ratio;
}
namespace npstat {
SmoothedDAgostiniUnfold1D::SmoothedDAgostiniUnfold1D(
const Matrix<double>& responseMatrix,
const LocalPolyFilter1D& f,
+ const bool i_useConvolutions,
const bool i_useMultinomialCovariance,
const bool i_smoothLastIteration,
const double i_convergenceEpsilon,
const unsigned i_maxIterations)
: AbsUnfold1D(responseMatrix),
filter_(0),
convergenceEpsilon_(0.0),
maxIterations_(i_maxIterations),
lastIterations_(0),
lastEPIterations_(0),
+ useConvolutions_(i_useConvolutions),
useMultinomialCovariance_(i_useMultinomialCovariance),
smoothLast_(i_smoothLastIteration),
v1_(responseMatrix.nColumns()),
v2_(responseMatrix.nColumns()),
v3_(std::max(responseMatrix.nRows(), responseMatrix.nColumns()))
{
setConvergenceEpsilon(i_convergenceEpsilon);
setFilter(f);
}
void SmoothedDAgostiniUnfold1D::setFilter(const LocalPolyFilter1D& f)
{
validateDimensions(responseMatrix().nRows(), f.dataLen());
filter_ = &f;
}
void SmoothedDAgostiniUnfold1D::setConvergenceEpsilon(const double eps)
{
convergenceEpsilon_ = eps;
if (convergenceEpsilon_ < 0.0)
convergenceEpsilon_ = 0.0;
}
bool SmoothedDAgostiniUnfold1D::unfold(
const double* observed, const unsigned lenObserved,
double* unfolded, const unsigned lenUnfolded,
Matrix<double>* unfoldedCovarianceFunction)
{
validateDimensions(lenObserved, lenUnfolded);
validateDensity(observed, lenObserved);
assert(unfolded);
// Initialize the "previous" approximation
double* prev = &v1_[0];
clearBuffer(prev, lenUnfolded);
// Set the "next" approximation to the initial one
if (getInitialApproximation().empty())
buildUniformInitialApproximation(observed, lenObserved, &v2_);
else
v2_ = getInitialApproximation();
double* next = &v2_[0];
bool converged = false;
for (lastIterations_=0U; lastIterations_<maxIterations_ && !converged;
++lastIterations_)
{
std::swap(prev, next);
update(observed, lenObserved, prev, next, lenUnfolded, true);
converged = probDelta(prev, next, lenUnfolded)<=convergenceEpsilon_;
}
// One more cycle in case we do not want to smooth the last iteration
if (!smoothLast_)
{
std::swap(prev, next);
update(observed, lenObserved, prev, next, lenUnfolded, false);
}
// Fill out the unfolded result
for (unsigned i=0; i<lenUnfolded; ++i)
unfolded[i] = next[i];
if (unfoldedCovarianceFunction)
{
// We will also need to calculate the result covariance matrix
unfoldedCovarianceFunction->uninitialize();
// Estimate yhat for the smoothed spectrum. This yhat
// has to be used to get the the error propagation matrix.
const double *smoothed = smoothLast_ ? next : prev;
const Matrix<double>& K = responseMatrix();
const Matrix<double>& yS = K.timesVector(smoothed, lenUnfolded);
bool conv = false;
const Matrix<double>& J = errorPropagationMatrix(
observed, lenObserved, smoothed, lenUnfolded, yS.data(), &conv);
converged = converged && conv;
// To construct the covariance matrix for the
// observations, we need the final yhat
if (smoothLast_)
{
// We have the final yhat already
const Matrix<double>& mcov = observationCovariance(
yS.data(), lenObserved);
*unfoldedCovarianceFunction = J*mcov*J.T();
}
else
{
// Must recalculate yhat
const Matrix<double>& yFinal = K.timesVector(next, lenUnfolded);
const Matrix<double>& mcov = observationCovariance(
yFinal.data(), lenObserved);
*unfoldedCovarianceFunction = J*mcov*J.T();
}
}
return converged;
}
Matrix<double> SmoothedDAgostiniUnfold1D::unfoldingMatrix0(
const double* unfolded, const double* yhat) const
{
Matrix<double> umat(transposedResponseMatrix());
const unsigned nRows = umat.nRows();
const unsigned nCols = umat.nColumns();
const double* eff = &efficiency()[0];
for (unsigned row=0; row<nRows; ++row)
{
double* rowData = umat[row];
const double factor = unfolded[row]/eff[row];
for (unsigned col=0; col<nCols; ++col)
rowData[col] *= factor/yhat[col];
}
return umat;
}
Matrix<double> SmoothedDAgostiniUnfold1D::observationCovariance(
const double* yhat, const unsigned lenObserved) const
{
Matrix<double> mcov(lenObserved, lenObserved, 0);
if (useMultinomialCovariance_)
{
const double N = std::accumulate(yhat, yhat+lenObserved, 0.0L);
for (unsigned i=0; i<lenObserved; ++i)
{
const double pi = yhat[i]/N;
double* mat = mcov[i];
mat[i] = N*pi*(1.0 - pi);
for (unsigned j=0; j<i; ++j)
{
const double c = -pi*yhat[j];
mat[j] = c;
mcov[j][i] = c;
}
}
}
else
{
for (unsigned i=0; i<lenObserved; ++i)
mcov[i][i] = yhat[i];
mcov.tagAsDiagonal();
}
return mcov;
}
Matrix<double> SmoothedDAgostiniUnfold1D::errorPropagationMatrix(
const double* observed, const unsigned lenObserved,
const double* unfolded, const unsigned lenUnfolded,
const double* yHat, bool* convergedPtr) const
{
const Matrix<double>& K = responseMatrix();
const Matrix<double>& S = filter_->getFilterMatrix();
const Matrix<double>& m0 = unfoldingMatrix0(unfolded, yHat);
const double* eff = &efficiency()[0];
Matrix<double> diag(lenUnfolded, lenUnfolded, 0);
double* tmp = const_cast<double*>(&v3_[0]);
for (unsigned iobs=0; iobs<lenObserved; ++iobs)
tmp[iobs] = observed[iobs]/yHat[iobs];
{
const Matrix<double>& numMat = K.rowMultiply(tmp, lenObserved);
const double* num = numMat.data();
for (unsigned i=0; i<lenUnfolded; ++i)
diag[i][i] = num[i]/eff[i];
}
for (unsigned iobs=0; iobs<lenObserved; ++iobs)
tmp[iobs] /= yHat[iobs];
Matrix<double> mat(lenUnfolded, lenUnfolded);
for (unsigned m=0; m<lenUnfolded; ++m)
{
const double factor = unfolded[m]/eff[m];
for (unsigned r=0; r<lenUnfolded; ++r)
{
long double sum = 0.0L;
for (unsigned i=0; i<lenObserved; ++i)
sum += K[i][m]*tmp[i]*K[i][r];
mat[m][r] = factor*sum;
}
}
const Matrix<double>& rmat = diag - mat;
Matrix<double> buf1(lenUnfolded, lenObserved, 0);
Matrix<double> buf2(buf1);
Matrix<double>* prev = &buf1;
Matrix<double>* next = &buf2;
// Error propagation matrix for iterations that run until convergence
bool converged = false;
double oldNorm = 0.0;
for (lastEPIterations_ = 0U;
lastEPIterations_<maxIterations_ && !converged;
++lastEPIterations_)
{
std::swap(prev, next);
*next = S*(m0 + rmat * *prev);
const double del = (*next - *prev).frobeniusNorm();
const double nn = next->frobeniusNorm();
const double pn = oldNorm;
oldNorm = nn;
const double rat = del*2.0/(nn + pn);
converged = rat < convergenceEpsilon_;
}
*convergedPtr = converged;
if (smoothLast_)
return *next;
else
return m0 + rmat * *next;
}
void SmoothedDAgostiniUnfold1D::update(
const double* observed, const unsigned lenObserved,
const double* prev, double* next, const unsigned lenUnfolded,
const bool performSmoothing) const
{
// Perform one D'Agostini iteration
const Matrix<double>& K = responseMatrix();
const Matrix<double>& yhatMat = K.timesVector(prev, lenUnfolded);
const double* yhat = yhatMat.data();
double* tmp = &v3_[0];
for (unsigned iobs=0; iobs<lenObserved; ++iobs)
{
if (yhat[iobs] <= 0.0)
{
if (observed[iobs] == 0.0)
tmp[iobs] = 0.0;
else
{
std::ostringstream os;
os << "In npstat::SmoothedDAgostiniUnfold1D::update: "
<< yhat[iobs] << " entries predicted, "
<< observed[iobs] << " observed for bin " << iobs
<< ". You need to change either the response "
<< "matrix or the initial approximation.";
throw std::runtime_error(os.str());
}
}
else
tmp[iobs] = observed[iobs]/yhat[iobs];
}
const Matrix<double>& numMat = K.rowMultiply(tmp, lenObserved);
const double* num = numMat.data();
const double* eff = &efficiency()[0];
double* out = performSmoothing ? tmp : next;
for (unsigned iunf=0; iunf<lenUnfolded; ++iunf)
out[iunf] = prev[iunf]*num[iunf]/eff[iunf];
if (performSmoothing)
{
// Apply the smoothing procedure
- filter_->filter(out, lenUnfolded, next);
+ if (useConvolutions_)
+ filter_->convolve(out, lenUnfolded, next);
+ else
+ filter_->filter(out, lenUnfolded, next);
// Make sure that smoothing result is non-negative
// and restore normalization
const double sum = std::accumulate(out, out+lenUnfolded, 0.0L);
normalizeArrayAsDensity(next, lenUnfolded, 1.0/sum);
}
}
}
Index: trunk/npstat/stat/SmoothedDAgostiniUnfold1D.hh
===================================================================
--- trunk/npstat/stat/SmoothedDAgostiniUnfold1D.hh (revision 187)
+++ trunk/npstat/stat/SmoothedDAgostiniUnfold1D.hh (revision 188)
@@ -1,150 +1,160 @@
#ifndef NPSTAT_SMOOTHEDDAGOSTINIUNFOLD1D_HH_
#define NPSTAT_SMOOTHEDDAGOSTINIUNFOLD1D_HH_
/*!
// \file SmoothedDAgostiniUnfold1D.hh
//
// \brief Implements D'Agostini unfolding with smoothing
//
// Author: I. Volobouev
//
// May 2014
*/
#include "npstat/stat/AbsUnfold1D.hh"
#include "npstat/stat/LocalPolyFilter1D.hh"
namespace npstat {
class SmoothedDAgostiniUnfold1D : public AbsUnfold1D
{
public:
/**
// The constructor arguments are:
//
// responseMatrix -- Naturally, the problem response matrix.
//
// filter -- The filter to use for smoothing the
// unfolded values. This object will not
// make a copy of the filter. It is
// a responsibility of the caller to ensure
// that the argument filter exists while
// this object is in use.
//
+ // useConvolutions -- If "true", the code will call the
+ // "convolve" method of the filter rather
+ // than its "filter" method.
+ //
// useMultinomialCovariance -- Specifies whether we should use
// multinomial distribution to estimate
// covariance of fitted observations
// (otherwise Poisson assumption is used).
//
// smoothLastIter -- If "false", smoothing will not be
// applied after the last iteration.
// Setting this parameter to "false" is
// not recommended for production results
// because it becomes unclear how to
// compare such results with models.
//
// convergenceEpsilon -- Convergence criterion for various
// iterations.
//
// maxIterations -- Maximum number of iterations allowed
// (both for D'Agostini iterations and for
// the code estimating the error propagation
// matrix).
*/
SmoothedDAgostiniUnfold1D(const Matrix<double>& responseMatrix,
const LocalPolyFilter1D& filter,
+ bool useConvolutions,
bool useMultinomialCovariance,
bool smoothLastIter = true,
double convergenceEpsilon = 1.0e-10,
unsigned maxIterations = 100000U);
inline virtual ~SmoothedDAgostiniUnfold1D() {}
/** Change maximum number of iterations */
inline void setMaxIterations(const unsigned u) {maxIterations_ = u;}
/** Switch between multinomial/Poisson covariance for observed space */
inline void useMultinomialCovariance(const bool b)
{useMultinomialCovariance_ = b;}
/** Switch between smoothing/not smoothing last iteration */
inline void smoothLastIteration(const bool b) {smoothLast_ = b;}
+ /** Switch between using filtering/convolution */
+ inline void useConvolutions(const bool b) {useConvolutions_ = b;}
+
/** Change the filter used. For use in bandwidth scans, etc. */
void setFilter(const LocalPolyFilter1D& f);
/** Change convergence criterion */
void setConvergenceEpsilon(double eps);
//@{
/** Simple inspector of the object properties */
inline const LocalPolyFilter1D& getFilter() const {return *filter_;}
inline double convergenceEpsilon() const {return convergenceEpsilon_;}
inline unsigned maxIterations() const {return maxIterations_;}
inline bool usingMultinomialCovariance() const
{return useMultinomialCovariance_;}
inline bool smoothingLastIteration() const {return smoothLast_;}
+ inline bool usingConvolutions() const {return useConvolutions_;}
//@}
/**
// Last number of iterations used to calculate the unfolded results.
// This number will be filled after each "unfold" call.
*/
inline unsigned lastNIterations() const {return lastIterations_;}
/**
// Last number of iterations used to calculate the error
// propagation matrix
*/
inline unsigned lastEPIterations() const {return lastEPIterations_;}
/**
// The main unfolding method. The "unfoldedCovarianceFunction"
// parameter can be NULL in which case the corresponding calculation
// is not performed.
*/
virtual bool unfold(const double* observed, unsigned lenObserved,
double* unfolded, unsigned lenUnfolded,
Matrix<double>* unfoldedCovarianceFunction);
private:
// D'Agostini iteration
void update(const double* observed, unsigned lenObserved,
const double* prev, double* next, unsigned lenUnfolded,
bool performSmoothing) const;
// The matrix which D'Agostini calls "unfolding matrix" in his
// original paper. This matrix takes into account neither the fixed
// point iteration nor the smoothing, so it can not be used for
// error propagation. It was wrong to use it for error propagation
// even in D'Agostini's case -- because changes in the data affect
// every iteration, not just the last one.
Matrix<double> unfoldingMatrix0(const double* unfolded,
const double* yHat) const;
// Error propagation matrix which does take into account
// the fixed point iteration and the smoothing
Matrix<double> errorPropagationMatrix(const double* observed,
unsigned lenObserved,
const double* unfolded,
unsigned lenUnfolded,
const double* yHat,
bool *converged) const;
// Covariance matrix of the observations
Matrix<double> observationCovariance(const double* yhat,
unsigned lenObserved) const;
const LocalPolyFilter1D* filter_;
double convergenceEpsilon_;
unsigned maxIterations_;
unsigned lastIterations_;
mutable unsigned lastEPIterations_;
+ bool useConvolutions_;
bool useMultinomialCovariance_;
bool smoothLast_;
// Memory buffers
std::vector<double> v1_, v2_;
mutable std::vector<double> v3_;
};
}
#endif // NPSTAT_SMOOTHEDDAGOSTINIUNFOLD1D_HH_
File Metadata
Details
Attached
Mime Type
text/x-diff
Expires
Tue, Sep 30, 5:42 AM (1 h, 8 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6566265
Default Alt Text
(18 KB)
Attached To
Mode
rNPSTATSVN npstatsvn
Attached
Detach File
Event Timeline
Log In to Comment