From 6d280ad011f7878ab47f72dcae7507477e66935f Mon Sep 17 00:00:00 2001 From: Giovanni Date: Mon, 20 Jan 2014 11:03:39 +0100 Subject: [PATCH 1/4] Expose who's the base and who's the exponent --- interface/HGGRooPdfs.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/interface/HGGRooPdfs.h b/interface/HGGRooPdfs.h index 1185737bd5d7c..cd168c17f012f 100644 --- a/interface/HGGRooPdfs.h +++ b/interface/HGGRooPdfs.h @@ -33,6 +33,8 @@ class RooPower : public RooAbsPdf { Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const ; Double_t analyticalIntegral(Int_t code, const char* rangeName=0) const ; + const RooAbsReal & base() const { return x.arg(); } + const RooAbsReal & exponent() const { return c.arg(); } protected: RooRealProxy x; From 31056390a79106cc73bc516817d0a1660b950117 Mon Sep 17 00:00:00 2001 From: Giovanni Date: Mon, 20 Jan 2014 14:11:00 +0100 Subject: [PATCH 2/4] vectorized MultiPdf, nexted AddPdfs, RooExponential and RooPower --- interface/CachingMultiPdf.h | 39 +++++++++++ interface/CachingNLL.h | 21 ++++-- interface/VectorizedSimplePdfs.h | 31 +++++++++ src/CachingMultiPdf.cc | 116 +++++++++++++++++++++++++++++++ src/CachingNLL.cc | 108 +++++++++++++++++++++------- src/VectorizedSimplePdfs.cc | 55 +++++++++++++++ src/vectorized.cc | 21 ++++++ src/vectorized.h | 6 ++ 8 files changed, 367 insertions(+), 30 deletions(-) create mode 100644 interface/CachingMultiPdf.h create mode 100644 interface/VectorizedSimplePdfs.h create mode 100644 src/CachingMultiPdf.cc create mode 100644 src/VectorizedSimplePdfs.cc diff --git a/interface/CachingMultiPdf.h b/interface/CachingMultiPdf.h new file mode 100644 index 0000000000000..c78a621fa9462 --- /dev/null +++ b/interface/CachingMultiPdf.h @@ -0,0 +1,39 @@ +#ifndef CachingMultiPdf_h +#define CachingMultiPdf_h + +#include "../interface/RooMultiPdf.h" +#include "../interface/CachingNLL.h" +#include +#include +#include + +namespace cacheutils { + class CachingMultiPdf : public CachingPdfBase { + public: + CachingMultiPdf(const RooMultiPdf &pdf, const RooArgSet &obs) ; + ~CachingMultiPdf() ; + virtual const std::vector & eval(const RooAbsData &data) ; + const RooAbsReal *pdf() const { return pdf_; } + virtual void setDataDirty() ; + protected: + const RooMultiPdf * pdf_; + boost::ptr_vector cachingPdfs_; + }; + + class CachingAddPdf : public CachingPdfBase { + public: + CachingAddPdf(const RooAddPdf &pdf, const RooArgSet &obs) ; + ~CachingAddPdf() ; + virtual const std::vector & eval(const RooAbsData &data) ; + const RooAbsReal *pdf() const { return pdf_; } + virtual void setDataDirty() ; + protected: + const RooAddPdf * pdf_; + std::vector coeffs_; + boost::ptr_vector cachingPdfs_; + std::vector work_; + }; + +} // namespace + +#endif diff --git a/interface/CachingNLL.h b/interface/CachingNLL.h index 60e6c71870da0..a0a6b9c410b82 100644 --- a/interface/CachingNLL.h +++ b/interface/CachingNLL.h @@ -29,6 +29,8 @@ namespace cacheutils { private: std::vector vars_; std::vector vals_; + std::vector cats_; + std::vector states_; }; // Part zero point five: Cache of pdf values for different parameters @@ -56,14 +58,22 @@ namespace cacheutils { Item *items[MaxItems_]; }; // Part one: cache all values of a pdf -class CachingPdf { +class CachingPdfBase { + public: + CachingPdfBase() {} + virtual ~CachingPdfBase() {} + virtual const std::vector & eval(const RooAbsData &data) = 0; + virtual const RooAbsReal *pdf() const = 0; + virtual void setDataDirty() = 0; +}; +class CachingPdf : public CachingPdfBase { public: CachingPdf(RooAbsReal *pdf, const RooArgSet *obs) ; CachingPdf(const CachingPdf &other) ; virtual ~CachingPdf() ; - const std::vector & eval(const RooAbsData &data) ; + virtual const std::vector & eval(const RooAbsData &data) ; const RooAbsReal *pdf() const { return pdf_; } - void setDataDirty() { lastData_ = 0; } + virtual void setDataDirty() { lastData_ = 0; } protected: const RooArgSet *obs_; RooAbsReal *pdfOriginal_; @@ -91,6 +101,7 @@ class OptimizedCachingPdfT : public CachingPdf { VPdfT *vpdf_; }; +CachingPdfBase * makeCachingPdf(RooAbsReal *pdf, const RooArgSet *obs) ; class CachingAddNLL : public RooAbsReal { public: @@ -118,10 +129,10 @@ class CachingAddNLL : public RooAbsReal { std::vector weights_; double sumWeights_; mutable std::vector coeffs_; - mutable boost::ptr_vector pdfs_; + mutable boost::ptr_vector pdfs_; mutable boost::ptr_vector prods_; mutable std::vector integrals_; - mutable std::vector > multiPdfs_; + mutable std::vector > multiPdfs_; mutable std::vector partialSum_; mutable std::vector workingArea_; mutable bool isRooRealSum_, fastExit_; diff --git a/interface/VectorizedSimplePdfs.h b/interface/VectorizedSimplePdfs.h new file mode 100644 index 0000000000000..64b41a271a1cb --- /dev/null +++ b/interface/VectorizedSimplePdfs.h @@ -0,0 +1,31 @@ +#ifndef VectorizedExponential_h +#define VectorizedExponential_h + +#include +#include +#include "../interface/HGGRooPdfs.h" +#include + +class VectorizedExponential { + public: + VectorizedExponential(const RooExponential &pdf, const RooAbsData &data) ; + void fill(std::vector &out) const ; + private: + const RooRealVar * x_; + const RooAbsReal * lambda_; + std::vector xvals_; + mutable std::vector work_; +}; + +class VectorizedPower { + public: + VectorizedPower(const RooPower &pdf, const RooAbsData &data) ; + void fill(std::vector &out) const ; + private: + const RooRealVar * x_; + const RooAbsReal * exponent_; + std::vector xvals_; + mutable std::vector work_; +}; + +#endif diff --git a/src/CachingMultiPdf.cc b/src/CachingMultiPdf.cc new file mode 100644 index 0000000000000..ee25217062370 --- /dev/null +++ b/src/CachingMultiPdf.cc @@ -0,0 +1,116 @@ +#include "../interface/CachingMultiPdf.h" +#include "vectorized.h" + +// Uncomment do do regression testing wrt uncached multipdf +//#define CachingMultiPdf_VALIDATE + +cacheutils::CachingMultiPdf::CachingMultiPdf(const RooMultiPdf &pdf, const RooArgSet &obs) : + pdf_(&pdf) +{ + //std::cout << "Making a CachingMultiPdf for " << pdf.GetName() << " with " << pdf_->getNumPdfs() << " pdfs." << std::endl; + for (int i = 0, n = pdf_->getNumPdfs(); i < n; ++i) { + cachingPdfs_.push_back(makeCachingPdf(pdf_->getPdf(i), &obs)); + //cachingPdfs_.push_back(new CachingPdf(pdf_->getPdf(i), &obs)); + //std::cout << " MultiPdfAdding " << pdf.GetName() << "[" << i << "]: " << pdf_->getPdf(i)->ClassName() << " " << pdf_->getPdf(i)->GetName() << " using " << typeid(cachingPdfs_.back()).name() << std::endl; + } +#ifdef CachingMultiPdf_VALIDATE + cachingPdfs_.push_back(new CachingPdf((RooAbsPdf*)&pdf,&obs)); +#endif + +} + +cacheutils::CachingMultiPdf::~CachingMultiPdf() +{ +} + +const std::vector & cacheutils::CachingMultiPdf::eval(const RooAbsData &data) +{ + const std::vector & ret = cachingPdfs_[pdf_->getCurrentIndex()].eval(data); +#ifdef CachingMultiPdf_VALIDATE + cachingPdfs_.back().setDataDirty(); + const std::vector & chk = cachingPdfs_.back().eval(data); + double normratio = 0.0; + for (unsigned int i = 0, n = ret.size(); i < n; ++i) { + if (chk[i] > 1e-6) normratio = max(normratio, abs(ret[i]/chk[i]-1)); + } + printf("%-70s: normalization match %+13.10f, selected pdf %s\n", pdf_->GetName(), normratio, pdf_->getCurrentPdf()->ClassName()); fflush(stdout); + if (normratio > 1e-5 && typeid(*pdf_->getCurrentPdf()) == typeid(RooAddPdf) ) { + RooAddPdf *add = (RooAddPdf*)pdf_->getCurrentPdf(); + for (unsigned int i = 0, n = add->pdfList().getSize(); i < n; ++i) { + printf("%-70s: RooAdd[%d] is a %s\n", pdf_->GetName(), i, add->pdfList().at(i)->ClassName()); + } + } + if (normratio > 1e-5) { + double sum = 0.0, sumc = 0; + std::cout << "Before normalization " << std::endl; + for (unsigned int i = 0, n = ret.size(); i < n; ++i) { + printf("%s[%6d] %18.10f %18.10f %+13.10f %+13.10f\n", pdf_->GetName(), i, ret[i], chk[i], ret[i]-chk[i], ret[i]/chk[i]-1); + sum += ret[i]; sumc += chk[i]; + } + printf("%s[%6s] %18.10f %18.10f %+13.10f %+13.10f\n", pdf_->GetName(), "TOTAL", sum, sumc, sum-sumc, sum/sumc-1); + } + return chk; +#endif + return ret; +} + +void cacheutils::CachingMultiPdf::setDataDirty() +{ + for (CachingPdfBase &pdf : cachingPdfs_) { + pdf.setDataDirty(); + } +} + +cacheutils::CachingAddPdf::CachingAddPdf(const RooAddPdf &pdf, const RooArgSet &obs) : + pdf_(&pdf) +{ + const RooArgList & pdfs = pdf.pdfList(); + const RooArgList & coeffs = pdf.coefList(); + //std::cout << "Making a CachingAddPdf for " << pdf.GetName() << " with " << pdfs.getSize() << " pdfs, " << coeffs.getSize() << " coeffs." << std::endl; + for (int i = 0, n = pdfs.getSize(); i < n; ++i) { + RooAbsPdf *pdfi = (RooAbsPdf*) pdfs.at(i); + cachingPdfs_.push_back(makeCachingPdf(pdfi, &obs)); + //std::cout << " AddPdfAdding " << pdf.GetName() << "[" << i << "]: " << pdfi->ClassName() << " " << pdfi->GetName() << " using " << typeid(cachingPdfs_.back()).name() << std::endl; + } + for (int i = 0, n = coeffs.getSize(); i < n; ++i) { + coeffs_.push_back((RooAbsReal*)coeffs.at(i)); + } +} + +cacheutils::CachingAddPdf::~CachingAddPdf() +{ +} + +const std::vector & cacheutils::CachingAddPdf::eval(const RooAbsData &data) +{ + double coefSum = 0.0; + std::vector coefVals(cachingPdfs_.size()); + for (int i = 0, n = coeffs_.size(); i < n; ++i) { + coefVals[i] = coeffs_[i]->getVal(); + coefSum += coefVals[i]; + } + if (coeffs_.size() != coefVals.size()) { + coefVals.back() = 1.0 - coefSum; + } else if (coefSum != 1.0) { + for (int i = 0, n = coeffs_.size(); i < n; ++i) { + coefVals[i] /= coefSum; + } + } + const std::vector & one = cachingPdfs_.front().eval(data); + unsigned int size = one.size(); + work_.resize(size); + std::fill(work_.begin(), work_.end(), 0.0); + for (int i = 0, n = coefVals.size(); i < n; ++i) { + vectorized::mul_add(size, coefVals[i], &(i ? cachingPdfs_[i].eval(data) : one)[0], &work_[0]); + } + //std::cout << "Evaluated " << pdf_->GetName() << " of type CachingAddPdf" << std::endl; + return work_; +} + +void cacheutils::CachingAddPdf::setDataDirty() +{ + for (CachingPdfBase &pdf : cachingPdfs_) { + pdf.setDataDirty(); + } +} + diff --git a/src/CachingNLL.cc b/src/CachingNLL.cc index bd4ddaf55663c..3ada57ce2c8b7 100644 --- a/src/CachingNLL.cc +++ b/src/CachingNLL.cc @@ -8,12 +8,16 @@ #include <../interface/RooMultiPdf.h> #include <../interface/VerticalInterpHistPdf.h> #include <../interface/VectorizedGaussian.h> +#include <../interface/VectorizedSimplePdfs.h> +#include <../interface/CachingMultiPdf.h> #include "vectorized.h" namespace cacheutils { typedef OptimizedCachingPdfT CachingHistPdf; typedef OptimizedCachingPdfT CachingHistPdf2; typedef OptimizedCachingPdfT CachingGaussPdf; + typedef OptimizedCachingPdfT CachingExpoPdf; + typedef OptimizedCachingPdfT CachingPowerPdf; class ReminderSum : public RooAbsReal { public: @@ -35,6 +39,9 @@ namespace cacheutils { //---- Uncomment this to get a '.' printed every some evals //#define TRACE_NLL_EVALS +//---- Uncomment this to get a total of the evals done +//#define TRACE_NLL_EVAL_COUNT + //---- Uncomment this and run with --perfCounters to get cache statistics //#define DEBUG_CACHE @@ -86,6 +93,9 @@ namespace { #define TRACE_NLL(x) #endif +#ifdef TRACE_NLL_EVAL_COUNT +namespace { unsigned long CachingSimNLLEvalCount = 0; } +#endif cacheutils::ArgSetChecker::ArgSetChecker(const RooAbsCollection &set) { @@ -97,6 +107,12 @@ cacheutils::ArgSetChecker::ArgSetChecker(const RooAbsCollection &set) if (rrv) { // && !rrv->isConstant()) { vars_.push_back(rrv); vals_.push_back(rrv->getVal()); + continue; + } + RooCategory *cat = dynamic_cast(a); + if (cat) { + cats_.push_back(cat); + states_.push_back(cat->getIndex()); } } } @@ -104,18 +120,30 @@ cacheutils::ArgSetChecker::ArgSetChecker(const RooAbsCollection &set) bool cacheutils::ArgSetChecker::changed(bool updateIfChanged) { + bool changed = false; std::vector::const_iterator it = vars_.begin(), ed = vars_.end(); std::vector::iterator itv = vals_.begin(); - bool changed = false; for ( ; it != ed; ++it, ++itv) { double val = (*it)->getVal(); if (val != *itv) { - //std::cerr << "var::CachingPdfable " << (*it)->GetName() << " changed: " << *itv << " -> " << val << std::endl; + //std::cerr << "var::CachingPdf " << (*it)->GetName() << " changed: " << *itv << " -> " << val << std::endl; changed = true; if (updateIfChanged) { *itv = val; } else break; } } + std::vector::const_iterator itc = cats_.begin(), edc = cats_.end(); + std::vector::iterator itvc = states_.begin(); + for ( ; itc != edc; ++itc, ++itvc) { + int val = (*itc)->getIndex(); + if (val != *itvc) { + //std::cerr << "cat::CachingPdf " << (*itc)->GetName() << " changed: " << *itvc << " -> " << val << std::endl; + changed = true; + if (updateIfChanged) { *itvc = val; } + else break; + } + } + return changed; } @@ -132,6 +160,7 @@ cacheutils::ValuesCache::ValuesCache(const RooAbsReal &pdf, const RooArgSet &obs { assert(size <= MaxItems_); std::auto_ptr params(pdf.getParameters(obs)); + //std::cout << "Parameters for pdf " << pdf.GetName() << " (" << pdf.ClassName() << "):"; params->Print(""); items[0] = new Item(*params); } @@ -231,7 +260,7 @@ cacheutils::CachingPdf::eval(const RooAbsData &data) bool newdata = (lastData_ != &data); if (newdata) newData_(data); std::pair *, bool> hit = cache_.get(); - if (!hit.second) { + if (!hit.second) { realFill_(data, *hit.first); } return *hit.first; @@ -265,6 +294,9 @@ cacheutils::CachingPdf::realFill_(const RooAbsData &data, std::vector #ifdef DEBUG_CACHE PerfCounter::add("CachingPdf::realFill_ called"); #endif + //std::cout << "CachingPdf::realFill_ called for " << pdf_->GetName() << " (" << pdf_->ClassName() << ")\n"; + //utils::printPdf((RooAbsPdf*)pdf_); + int n = data.numEntries(); vals.resize(nonZeroWEntries_); // should be a no-op if size is already >= n. std::vector::iterator itv = vals.begin(); @@ -350,8 +382,6 @@ cacheutils::CachingAddNLL::clone(const char *name) const void cacheutils::CachingAddNLL::addPdfs_(RooAddPdf *addpdf, bool recursive, const RooArgList & basecoeffs) { - bool histNll = runtimedef::get("ADDNLL_HISTNLL"); - bool gaussNll = runtimedef::get("ADDNLL_GAUSSNLL"); int npdf = addpdf->pdfList().getSize(); //std::cout << "Unpacking RooAddPdf " << addpdf->GetName() << " with " << npdf << " components:" << std::endl; RooAbsReal *lastcoeff = 0; @@ -368,7 +398,7 @@ void cacheutils::CachingAddNLL::addPdfs_(RooAddPdf *addpdf, bool recursive, cons RooAddPdf *apdfi = static_cast(pdfi); RooArgList list(*coeff); if (basecoeffs.getSize()) list.add(basecoeffs); - //std::cout << " Invoking recursive unpack on " << i << ": RooAddPdf " << apdfi->GetName() << std::endl; + //std::cout << " Invoking recursive unpack on " << addpdf->GetName() << "[" << i << "]: RooAddPdf " << apdfi->GetName() << std::endl; addPdfs_(apdfi, recursive, list); continue; } @@ -381,17 +411,35 @@ void cacheutils::CachingAddNLL::addPdfs_(RooAddPdf *addpdf, bool recursive, cons coeffs_.push_back(&prods_.back()); } const RooArgSet *obs = data_->get(); - //std::cout << " Adding " << i << ": " << pdfi->ClassName() << " " << pdfi->GetName() << std::endl; - if (histNll && typeid(*pdfi) == typeid(FastVerticalInterpHistPdf)) { - pdfs_.push_back(new CachingHistPdf(pdfi, obs)); - } else if (histNll && typeid(*pdfi) == typeid(FastVerticalInterpHistPdf2)) { - pdfs_.push_back(new CachingHistPdf2(pdfi, obs)); - } else if (gaussNll && typeid(*pdfi) == typeid(RooGaussian)) { - pdfs_.push_back(new CachingGaussPdf(pdfi, obs)); - } else { - pdfs_.push_back(new CachingPdf(pdfi, obs)); - } + pdfs_.push_back(makeCachingPdf(pdfi,obs)); + //std::cout << " Adding " << addpdf->GetName() << "[" << i << "]: " << pdfi->ClassName() << " " << pdfi->GetName() << " using " << typeid(pdfs_.back()).name() << std::endl; + } +} + +cacheutils::CachingPdfBase * +cacheutils::makeCachingPdf(RooAbsReal *pdf, const RooArgSet *obs) { + static bool histNll = runtimedef::get("ADDNLL_HISTNLL"); + static bool gaussNll = runtimedef::get("ADDNLL_GAUSSNLL"); + static bool multiNll = runtimedef::get("ADDNLL_MULTINLL"); + + if (histNll && typeid(*pdf) == typeid(FastVerticalInterpHistPdf)) { + return new CachingHistPdf(pdf, obs); + } else if (histNll && typeid(*pdf) == typeid(FastVerticalInterpHistPdf2)) { + return new CachingHistPdf2(pdf, obs); + } else if (gaussNll && typeid(*pdf) == typeid(RooGaussian)) { + return new CachingGaussPdf(pdf, obs); + } else if (gaussNll && typeid(*pdf) == typeid(RooExponential)) { + return new CachingExpoPdf(pdf, obs); + } else if (gaussNll && typeid(*pdf) == typeid(RooPower)) { + return new CachingPowerPdf(pdf, obs); + } else if (multiNll && typeid(*pdf) == typeid(RooMultiPdf)) { + return new CachingMultiPdf(static_cast(*pdf), *obs); + } else if (multiNll && typeid(*pdf) == typeid(RooAddPdf)) { + return new CachingAddPdf(static_cast(*pdf), *obs); + } else { + return new CachingPdf(pdf, obs); } + } void @@ -436,13 +484,17 @@ cacheutils::CachingAddNLL::setup_() } // For multi pdf's need to reset the cache if index changed before evaluations + // unless they're being properly treated in the CachingPdf + static bool multiNll = runtimedef::get("ADDNLL_MULTINLL"); multiPdfs_.clear(); - for (auto itp = pdfs_.begin(), edp = pdfs_.end(); itp != edp; ++itp) { - bool isMultiPdf = itp->pdf()->IsA()->InheritsFrom(RooMultiPdf::Class()); - if (isMultiPdf) { - const RooMultiPdf *mpdf = dynamic_cast((*itp).pdf()); - multiPdfs_.push_back(std::make_pair(mpdf, &*itp)); - } + if (!multiNll) { + for (auto itp = pdfs_.begin(), edp = pdfs_.end(); itp != edp; ++itp) { + bool isMultiPdf = itp->pdf()->IsA()->InheritsFrom(RooMultiPdf::Class()); + if (isMultiPdf) { + const RooMultiPdf *mpdf = dynamic_cast((*itp).pdf()); + multiPdfs_.push_back(std::make_pair(mpdf, &*itp)); + } + } } } @@ -456,7 +508,7 @@ cacheutils::CachingAddNLL::evaluate() const // For multi pdf's need to reset the cache if index changed before evaluations if (!multiPdfs_.empty()) { - for (std::vector >::iterator itp = multiPdfs_.begin(), edp = multiPdfs_.end(); itp != edp; ++itp) { + for (std::vector >::iterator itp = multiPdfs_.begin(), edp = multiPdfs_.end(); itp != edp; ++itp) { bool hasChangedPdf = itp->first->checkIndexDirty(); if (hasChangedPdf) itp->second->setDataDirty(); } @@ -465,7 +517,7 @@ cacheutils::CachingAddNLL::evaluate() const std::fill( partialSum_.begin(), partialSum_.end(), 0.0 ); std::vector::iterator itc = coeffs_.begin(), edc = coeffs_.end(); - boost::ptr_vector::iterator itp = pdfs_.begin();//, edp = pdfs_.end(); + boost::ptr_vector::iterator itp = pdfs_.begin();//, edp = pdfs_.end(); std::vector::const_iterator itw, bgw = weights_.begin();//, edw = weights_.end(); std::vector::iterator its, bgs = partialSum_.begin(), eds = partialSum_.end(); double sumCoeff = 0; @@ -538,7 +590,7 @@ cacheutils::CachingAddNLL::evaluate() const // multipdfs want to add a correction factor to the NLL if (!multiPdfs_.empty()) { double correctionFactor = 0; - for (std::vector >::iterator itp = multiPdfs_.begin(), edp = multiPdfs_.end(); itp != edp; ++itp) { + for (std::vector >::iterator itp = multiPdfs_.begin(), edp = multiPdfs_.end(); itp != edp; ++itp) { correctionFactor += itp->first->getCorrection(); } // Add correction @@ -630,6 +682,9 @@ cacheutils::CachingSimNLL::~CachingSimNLL() for (std::vector::iterator it = constrainPdfsFast_.begin(), ed = constrainPdfsFast_.end(); it != ed; ++it, ++ito) { if (*ito) { delete *it; } } + #ifdef TRACE_NLL_EVAL_COUNT + std::cout << "CachingSimNLLEvalCount: " << ::CachingSimNLLEvalCount << std::endl; + #endif } void @@ -712,6 +767,9 @@ Double_t cacheutils::CachingSimNLL::evaluate() const { TRACE_POINT(params_) +#ifdef TRACE_NLL_EVAL_COUNT + ::CachingSimNLLEvalCount++; +#endif #ifdef DEBUG_CACHE PerfCounter::add("CachingSimNLL::evaluate called"); #endif diff --git a/src/VectorizedSimplePdfs.cc b/src/VectorizedSimplePdfs.cc new file mode 100644 index 0000000000000..5bc4c21c1f23e --- /dev/null +++ b/src/VectorizedSimplePdfs.cc @@ -0,0 +1,55 @@ +#include "../interface/VectorizedSimplePdfs.h" +#include "RooMath.h" +#include "vectorized.h" +#include +#include +#include + +VectorizedExponential::VectorizedExponential(const RooExponential &pdf, const RooAbsData &data) +{ + RooArgSet obs(*data.get()); + std::auto_ptr params(pdf.getParameters(data)); + if (params->getSize() != 1) throw std::invalid_argument("Can't resolve which is the parameter of the exponential"); + + x_ = dynamic_cast(obs.first()); + lambda_ = dynamic_cast(params->first()); + + xvals_.reserve(data.numEntries()); + for (unsigned int i = 0, n = data.numEntries(); i < n; ++i) { + obs.assignValueOnly(*data.get(i), true); + if (data.weight()) xvals_.push_back(x_->getVal()); + } + work_.resize(xvals_.size()); +} + +void VectorizedExponential::fill(std::vector &out) const { + Double_t xmax = x_->getMax(), xmin = x_->getMin(), lambda = lambda_->getVal(); + Double_t norm = (lambda != 0 ? (std::exp(lambda*xmax)-std::exp(lambda*xmin))/lambda : xmax-xmin); + out.resize(xvals_.size()); + vectorized::exponentials(xvals_.size(), lambda, norm, &xvals_[0], &out[0], &work_[0]); +} + +VectorizedPower::VectorizedPower(const RooPower &pdf, const RooAbsData &data) +{ + RooArgSet obs(*data.get()); + if (obs.find(pdf.base()) == 0) throw std::invalid_argument("Dataset does not depend on the base of the power"); + x_ = dynamic_cast(obs.first()); + exponent_ = &pdf.exponent(); + + xvals_.reserve(data.numEntries()); + for (unsigned int i = 0, n = data.numEntries(); i < n; ++i) { + obs.assignValueOnly(*data.get(i), true); + if (data.weight()) xvals_.push_back(x_->getVal()); + } + work_.resize(xvals_.size()); +} + +void VectorizedPower::fill(std::vector &out) const { + Double_t xmax = x_->getMax(), xmin = x_->getMin(), exponent = exponent_->getVal(); + Double_t norm; + if (exponent == -1) norm = log(xmax) - log(xmin); + else if (exponent == 0) norm = xmax - xmin; + else norm = (std::pow(xmax,exponent+1) - std::pow(xmin,exponent+1))/(exponent + 1); + out.resize(xvals_.size()); + vectorized::powers(xvals_.size(), exponent, norm, &xvals_[0], &out[0], &work_[0]); +} diff --git a/src/vectorized.cc b/src/vectorized.cc index 87c06f4f09923..f7becd0395a79 100644 --- a/src/vectorized.cc +++ b/src/vectorized.cc @@ -38,3 +38,24 @@ void vectorized::gaussians(const uint32_t size, double mean, double sigma, doubl out[i] = inorm*workingArea2[i]; } } + +void vectorized::exponentials(const uint32_t size, double lambda, double norm, const double* __restrict__ xvals, double * __restrict__ out, double * __restrict__ workingArea) +{ + //out[i] = std::exp(xvals[i]*lambda) * nfact; nfact = 1.0/norm + double lognfact = -std::log(norm); + for (uint32_t i = 0; i < size; ++i) { + workingArea[i] = xvals[i] * lambda + lognfact; + } + vdt::fast_expv(size, workingArea, out); +} + +void vectorized::powers(const uint32_t size, double exponent, double norm, const double* __restrict__ xvals, double * __restrict__ out, double * __restrict__ workingArea) +{ + //out[i] = std::pow(xvals[i],exponent) * nfact; // nfact = 1.0/norm + double lognfact = -std::log(norm); + vdt::fast_logv(size, xvals, workingArea); + for (uint32_t i = 0; i < size; ++i) { + workingArea[i] = workingArea[i]*exponent + lognfact; + } + vdt::fast_expv(size, workingArea, out); +} diff --git a/src/vectorized.h b/src/vectorized.h index 4a76111cc984d..8f2a76ded187e 100644 --- a/src/vectorized.h +++ b/src/vectorized.h @@ -9,5 +9,11 @@ namespace vectorized { // gaussians void gaussians(const uint32_t size, double mean, double sigma, double norm, const double* __restrict__ xvals, double * __restrict__ out, double * __restrict__ workingArea, double * __restrict__ workingArea2) ; + + // exponentials + void exponentials(const uint32_t size, double lambda, double norm, const double* __restrict__ xvals, double * __restrict__ out, double * __restrict__ workingArea) ; + + // powers + void powers(const uint32_t size, double lambda, double norm, const double* __restrict__ xvals, double * __restrict__ out, double * __restrict__ workingArea) ; } From e7efb1cf00e9673bf38435d9daf84a8ce96ec13c Mon Sep 17 00:00:00 2001 From: Giovanni Date: Mon, 20 Jan 2014 18:38:09 +0100 Subject: [PATCH 3/4] fix bug in recursive addpdf splitting (the default constructor of RooAbsReal doesn't seem to work) --- src/CachingNLL.cc | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/CachingNLL.cc b/src/CachingNLL.cc index 3ada57ce2c8b7..b32d3e8af5b62 100644 --- a/src/CachingNLL.cc +++ b/src/CachingNLL.cc @@ -333,6 +333,7 @@ cacheutils::OptimizedCachingPdfT::realFill_(const RooAbsData &data, cacheutils::ReminderSum::ReminderSum(const char *name, const char *title, const RooArgList& sumSet) : + RooAbsReal(name,title), list_("deps","",this) { RooLinkedListIter iter(sumSet.iterator()); @@ -388,9 +389,10 @@ void cacheutils::CachingAddNLL::addPdfs_(RooAddPdf *addpdf, bool recursive, cons if (npdf == addpdf->coefList().getSize()) { lastcoeff = dynamic_cast(addpdf->coefList().at(npdf-1)); } else { - prods_.push_back(new ReminderSum("","", addpdf->coefList())); + prods_.push_back(new ReminderSum((std::string("reminder_of_")+addpdf->GetName()).c_str(),"", addpdf->coefList())); lastcoeff = & prods_.back(); } + //std::cout << " Last coefficient is a " << lastcoeff->ClassName() << " aka " << typeid(*lastcoeff).name() << ": "; lastcoeff->Print(""); for (int i = 0; i < npdf; ++i) { RooAbsReal * coeff = (i < npdf-1 ? dynamic_cast(addpdf->coefList().at(i)) : lastcoeff); RooAbsPdf * pdfi = dynamic_cast(addpdf->pdfList().at(i)); @@ -409,6 +411,7 @@ void cacheutils::CachingAddNLL::addPdfs_(RooAddPdf *addpdf, bool recursive, cons list.add(basecoeffs); prods_.push_back(new RooProduct("","",list)); coeffs_.push_back(&prods_.back()); + //std::cout << "Coefficient of " << pdfi->GetName() << std::endl; prods_.back().Print(""); } const RooArgSet *obs = data_->get(); pdfs_.push_back(makeCachingPdf(pdfi,obs)); From 205c5ee9369a303c0d101258480c6c0ea8bdd470 Mon Sep 17 00:00:00 2001 From: Giovanni Date: Mon, 20 Jan 2014 18:40:07 +0100 Subject: [PATCH 4/4] MultiPDFs: no need of setDataDirty if unpacking them, but still need to remember them to get the penalty term --- src/CachingNLL.cc | 29 ++++++++++++++++++----------- 1 file changed, 18 insertions(+), 11 deletions(-) diff --git a/src/CachingNLL.cc b/src/CachingNLL.cc index b32d3e8af5b62..f2f1d3d14e022 100644 --- a/src/CachingNLL.cc +++ b/src/CachingNLL.cc @@ -45,6 +45,9 @@ namespace cacheutils { //---- Uncomment this and run with --perfCounters to get cache statistics //#define DEBUG_CACHE +//---- Uncomment to dump PDF values inside CachingAddNLL +//#define LOG_ADDPDFS + //---- Uncomment to enable Kahan's summation (if enabled at runtime with --X-rtd = ... // http://en.wikipedia.org/wiki/Kahan_summation_algorithm //#define ADDNLL_KAHAN_SUM @@ -486,17 +489,12 @@ cacheutils::CachingAddNLL::setup_() if (rrv != 0) params_.add(*rrv); } - // For multi pdf's need to reset the cache if index changed before evaluations - // unless they're being properly treated in the CachingPdf - static bool multiNll = runtimedef::get("ADDNLL_MULTINLL"); multiPdfs_.clear(); - if (!multiNll) { - for (auto itp = pdfs_.begin(), edp = pdfs_.end(); itp != edp; ++itp) { - bool isMultiPdf = itp->pdf()->IsA()->InheritsFrom(RooMultiPdf::Class()); - if (isMultiPdf) { - const RooMultiPdf *mpdf = dynamic_cast((*itp).pdf()); - multiPdfs_.push_back(std::make_pair(mpdf, &*itp)); - } + for (auto itp = pdfs_.begin(), edp = pdfs_.end(); itp != edp; ++itp) { + bool isMultiPdf = itp->pdf()->IsA()->InheritsFrom(RooMultiPdf::Class()); + if (isMultiPdf) { + const RooMultiPdf *mpdf = dynamic_cast((*itp).pdf()); + multiPdfs_.push_back(std::make_pair(mpdf, &*itp)); } } @@ -510,7 +508,9 @@ cacheutils::CachingAddNLL::evaluate() const #endif // For multi pdf's need to reset the cache if index changed before evaluations - if (!multiPdfs_.empty()) { + // unless they're being properly treated in the CachingPdf + static bool multiNll = runtimedef::get("ADDNLL_MULTINLL"); + if (!multiNll && !multiPdfs_.empty()) { for (std::vector >::iterator itp = multiPdfs_.begin(), edp = multiPdfs_.end(); itp != edp; ++itp) { bool hasChangedPdf = itp->first->checkIndexDirty(); if (hasChangedPdf) itp->second->setDataDirty(); @@ -536,6 +536,13 @@ cacheutils::CachingAddNLL::evaluate() const } // get vals const std::vector &pdfvals = itp->eval(*data_); +#ifdef LOG_ADDPDFS + printf("%s coefficient %s (%s) = %20.15f\n", itp->pdf()->GetName(), (*itc)->GetName(), (*itc)->ClassName(), coeff); + //(*itc)->Print(""); + for (unsigned int i = 0, n = pdfvals.size(); i < n; ++i) { + if (i%84==0) printf("%-80s[%3d] = %20.15f\n", itp->pdf()->GetName(), i, pdfvals[i]); + } +#endif // update running sum // std::vector::const_iterator itv = pdfvals.begin(); // for (its = bgs; its != eds; ++its, ++itv) {