Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooUnbinnedL.cxx
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 * Authors:
4 * PB, Patrick Bos, Netherlands eScience Center, p.bos@esciencecenter.nl
5 *
6 * Copyright (c) 2021, CERN
7 *
8 * Redistribution and use in source and binary forms,
9 * with or without modification, are permitted according to the terms
10 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
11 */
12
13/**
14\file RooUnbinnedL.cxx
15\class RooUnbinnedL
16\ingroup Roofitcore
17
18A -log(likelihood) calculation from a dataset
19(assumed to be unbinned) and a PDF. The NLL is calculated as
20\f[
21 \sum_\mathrm{data} -\log( \mathrm{pdf}(x_\mathrm{data}))
22\f]
23In extended mode, a
24\f$ N_\mathrm{expect} - N_\mathrm{observed}*log(N_\mathrm{expect}) \f$ term is added.
25**/
26
28
29#include <RooAbsData.h>
30#include <RooAbsPdf.h>
31#include <RooAbsDataStore.h>
32#include <RooChangeTracker.h>
33#include <RooNaNPacker.h>
34#include <RooFit/Evaluator.h>
35
37
38namespace RooFit {
39namespace TestStatistics {
40
41namespace {
42
43RooAbsL::ClonePdfData clonePdfData(RooAbsPdf &pdf, RooAbsData &data, RooFit::EvalBackend evalBackend)
44{
45 if (evalBackend.value() == RooFit::EvalBackend::Value::Legacy) {
46 return {&pdf, &data};
47 }
48 // For the evaluation with the BatchMode, the pdf needs to be "compiled" for
49 // a given normalization set.
50 return {RooFit::Detail::compileForNormSet(pdf, *data.get()), &data};
51}
52
53} // namespace
54
56 RooFit::EvalBackend evalBackend)
57 : RooAbsL(clonePdfData(*pdf, *data, evalBackend), data->numEntries(), 1, extended)
58{
59 std::unique_ptr<RooArgSet> params(pdf->getParameters(data));
60 paramTracker_ = std::make_unique<RooChangeTracker>("chtracker", "change tracker", *params, true);
61
62 if (evalBackend.value() != RooFit::EvalBackend::Value::Legacy) {
63 evaluator_ = std::make_unique<RooFit::Evaluator>(*pdf_, evalBackend.value() == RooFit::EvalBackend::Value::Cuda);
64 std::stack<std::vector<double>>{}.swap(_vectorBuffers);
65 auto dataSpans =
66 RooFit::Detail::BatchModeDataHelpers::getDataSpans(*data, "", nullptr, /*skipZeroWeights=*/true,
67 /*takeGlobalObservablesFromData=*/false, _vectorBuffers);
68 for (auto const &item : dataSpans) {
69 evaluator_->setInput(item.first->GetName(), item.second, false);
70 }
71 }
72}
73
75 : RooAbsL(other),
76 apply_weight_squared(other.apply_weight_squared),
77 _first(other._first),
78 lastSection_(other.lastSection_),
79 cachedResult_(other.cachedResult_),
80 evaluator_(other.evaluator_)
81{
82 paramTracker_ = std::make_unique<RooChangeTracker>(*other.paramTracker_);
83}
84
86
87//////////////////////////////////////////////////////////////////////////////////
88
89/// Returns true if value was changed, false otherwise.
91{
92 if (apply_weight_squared != flag) {
94 return true;
95 }
96 // setValueDirty();
97 return false;
98}
99
100namespace {
101
102using ComputeResult = std::pair<ROOT::Math::KahanSum<double>, double>;
103
104// Copy of RooNLLVar::computeScalarFunc.
105ComputeResult computeScalarFunc(const RooAbsPdf *pdfClone, RooAbsData *dataClone, RooArgSet *normSet, bool weightSq,
106 std::size_t stepSize, std::size_t firstEvent, std::size_t lastEvent,
107 RooAbsPdf const *offsetPdf = nullptr)
108{
111 RooNaNPacker packedNaN(0.f);
112
113 for (auto i = firstEvent; i < lastEvent; i += stepSize) {
114 dataClone->get(i);
115
116 double weight = dataClone->weight(); // FIXME
117
118 if (0. == weight * weight)
119 continue;
120 if (weightSq)
121 weight = dataClone->weightSquared();
122
123 double logProba = pdfClone->getLogVal(normSet);
124
125 if (offsetPdf) {
126 logProba -= offsetPdf->getLogVal(normSet);
127 }
128
129 const double term = -weight * logProba;
130
131 kahanWeight.Add(weight);
132 kahanProb.Add(term);
133 packedNaN.accumulate(term);
134 }
135
136 if (packedNaN.getPayload() != 0.) {
137 // Some events with evaluation errors. Return "badness" of errors.
138 return {ROOT::Math::KahanSum<double>{packedNaN.getNaNWithPayload()}, kahanWeight.Sum()};
139 }
140
141 return {kahanProb, kahanWeight.Sum()};
142}
143
144// For now, almost exact copy of computeScalarFunc.
145ComputeResult computeBatchFunc(std::span<const double> probas, RooAbsData *dataClone, bool weightSq,
146 std::size_t stepSize, std::size_t firstEvent, std::size_t lastEvent)
147{
150 RooNaNPacker packedNaN(0.f);
151
152 for (auto i = firstEvent; i < lastEvent; i += stepSize) {
153 dataClone->get(i);
154
155 double weight = dataClone->weight();
156
157 if (0. == weight * weight)
158 continue;
159 if (weightSq)
160 weight = dataClone->weightSquared();
161
162 double logProba = std::log(probas[i]);
163 const double term = -weight * logProba;
164
165 kahanWeight.Add(weight);
166 kahanProb.Add(term);
167 packedNaN.accumulate(term);
168 }
169
170 if (packedNaN.getPayload() != 0.) {
171 // Some events with evaluation errors. Return "badness" of errors.
172 return {ROOT::Math::KahanSum<double>{packedNaN.getNaNWithPayload()}, kahanWeight.Sum()};
173 }
174
175 return {kahanProb, kahanWeight.Sum()};
176}
177
178} // namespace
179
180//////////////////////////////////////////////////////////////////////////////////
181/// Calculate and return likelihood on subset of data from firstEvent to lastEvent
182/// processed with a step size of 'stepSize'. If this an extended likelihood and
183/// and the zero event is processed the extended term is added to the return
184/// likelihood.
185///
187RooUnbinnedL::evaluatePartition(Section events, std::size_t /*components_begin*/, std::size_t /*components_end*/)
188{
189 // Throughout the calculation, we use Kahan's algorithm for summing to
190 // prevent loss of precision - this is a factor four more expensive than
191 // straight addition, but since evaluating the PDF is usually much more
192 // expensive than that, we tolerate the additional cost...
194 double sumWeight;
195 auto numEvalErrorsBefore = RooAbsReal::numEvalErrors();
196
197 // Do not reevaluate likelihood if parameters nor event range have changed
198 if (!paramTracker_->hasChanged(true) && events == lastSection_ &&
199 (cachedResult_.Sum() != 0 || cachedResult_.Carry() != 0))
200 return cachedResult_;
201
202 if (evaluator_) {
203 // Here, we have a memory allocation that should be avoided when this
204 // code needs to be optimized.
205 std::span<const double> probas = evaluator_->run();
206 std::tie(result, sumWeight) =
207 computeBatchFunc(probas, data_.get(), apply_weight_squared, 1, events.begin(N_events_), events.end(N_events_));
208 } else {
209 data_->store()->recalculateCache(nullptr, events.begin(N_events_), events.end(N_events_), 1, true);
210 std::tie(result, sumWeight) = computeScalarFunc(pdf_.get(), data_.get(), normSet_.get(), apply_weight_squared, 1,
211 events.begin(N_events_), events.end(N_events_));
212 }
213
214 // include the extended maximum likelihood term, if requested
215 if (extended_ && events.begin_fraction == 0) {
216 result += pdf_->extendedTerm(*data_, apply_weight_squared);
217 }
218
219 // If part of simultaneous PDF normalize probability over
220 // number of simultaneous PDFs: -sum(log(p/n)) = -sum(log(p)) + N*log(n)
221 if (sim_count_ > 1) {
222 result += sumWeight * log(1.0 * sim_count_);
223 }
224
225 // At the end of the first full calculation, wire the caches. This doesn't
226 // need to be done in BatchMode with the RooFit driver.
227 if (_first && !evaluator_) {
228 _first = false;
229 pdf_->wireAllCaches();
230 }
231
234 numEvalErrorsBefore == RooAbsReal::numEvalErrors()) {
236 lastSection_ = events;
237 }
238 return result;
239}
240
241} // namespace TestStatistics
242} // namespace RooFit
static RooNLLVar::ComputeResult computeScalarFunc(const RooAbsPdf *pdfClone, RooAbsData *dataClone, RooArgSet *normSet, bool weightSq, std::size_t stepSize, std::size_t firstEvent, std::size_t lastEvent, RooAbsPdf const *offsetPdf=nullptr)
bool _first
!
Definition RooNLLVar.h:46
std::pair< ROOT::Math::KahanSum< double >, double > ComputeResult
Definition RooNLLVar.h:25
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
The Kahan summation is a compensated summation algorithm, which significantly reduces numerical error...
Definition Util.h:122
T Sum() const
Definition Util.h:240
T Carry() const
Definition Util.h:250
void Add(T x)
Single-element accumulation. Will not vectorise.
Definition Util.h:165
RooFit::OwningPtr< RooArgSet > getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
Abstract base class for binned and unbinned datasets.
Definition RooAbsData.h:57
virtual double weight() const =0
virtual const RooArgSet * get() const
Definition RooAbsData.h:101
virtual double weightSquared() const =0
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
virtual double getLogVal(const RooArgSet *set=nullptr) const
Return the log of the current value with given normalization An error message is printed if the argum...
static ErrorLoggingMode evalErrorLoggingMode()
Return current evaluation error logging mode.
static Int_t numEvalErrors()
Return the number of logged evaluation errors since the last clearing.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Value value() const
std::shared_ptr< RooAbsData > data_
Definition RooAbsL.h:143
std::unique_ptr< RooArgSet > normSet_
Pointer to set with observables used for normalization.
Definition RooAbsL.h:144
std::shared_ptr< RooAbsPdf > pdf_
Definition RooAbsL.h:142
ROOT::Math::KahanSum< double > cachedResult_
bool setApplyWeightSquared(bool flag)
Returns true if value was changed, false otherwise.
ROOT::Math::KahanSum< double > evaluatePartition(Section events, std::size_t components_begin, std::size_t components_end) override
Calculate and return likelihood on subset of data from firstEvent to lastEvent processed with a step ...
std::unique_ptr< RooChangeTracker > paramTracker_
bool apply_weight_squared
Apply weights squared?
RooUnbinnedL(RooAbsPdf *pdf, RooAbsData *data, RooAbsL::Extended extended=RooAbsL::Extended::Auto, RooFit::EvalBackend evalBackend=RooFit::EvalBackend::Legacy())
std::stack< std::vector< double > > _vectorBuffers
std::shared_ptr< RooFit::Evaluator > evaluator_
! For batched evaluation
std::unique_ptr< T > compileForNormSet(T const &arg, RooArgSet const &normSet)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Definition JSONIO.h:26
A part of some range delimited by two fractional points between 0 and 1 (inclusive).
Definition RooAbsL.h:65
std::size_t begin(std::size_t N_total) const
Definition RooAbsL.h:73
std::size_t end(std::size_t N_total) const
Definition RooAbsL.h:75
Little struct that can pack a float into the unused bits of the mantissa of a NaN double.