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
18Class RooUnbinnedL implements a -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 <RooNLLVar.h> // RooNLLVar::ComputeScalar
33#include <RooChangeTracker.h>
34#include <RooNaNPacker.h>
35#include <RooFit/Evaluator.h>
36
37#include "../RooFit/BatchModeDataHelpers.h"
38
39namespace RooFit {
40namespace TestStatistics {
41
42namespace {
43
44RooAbsL::ClonePdfData clonePdfData(RooAbsPdf &pdf, RooAbsData &data, RooFit::EvalBackend evalBackend)
45{
46 if (evalBackend.value() == RooFit::EvalBackend::Value::Legacy) {
47 return {&pdf, &data};
48 }
49 // For the evaluation with the BatchMode, the pdf needs to be "compiled" for
50 // a given normalization set.
51 return {RooFit::Detail::compileForNormSet(pdf, *data.get()), &data};
52}
53
54} // namespace
55
57 RooFit::EvalBackend evalBackend)
58 : RooAbsL(clonePdfData(*pdf, *data, evalBackend), data->numEntries(), 1, extended)
59{
60 std::unique_ptr<RooArgSet> params(pdf->getParameters(data));
61 paramTracker_ = std::make_unique<RooChangeTracker>("chtracker", "change tracker", *params, true);
62
63 if (evalBackend.value() != RooFit::EvalBackend::Value::Legacy) {
64 evaluator_ = std::make_unique<RooFit::Evaluator>(*pdf_, evalBackend.value() == RooFit::EvalBackend::Value::Cuda);
65 std::stack<std::vector<double>>{}.swap(_vectorBuffers);
66 auto dataSpans =
67 RooFit::BatchModeDataHelpers::getDataSpans(*data, "", nullptr, /*skipZeroWeights=*/true,
68 /*takeGlobalObservablesFromData=*/false, _vectorBuffers);
69 for (auto const &item : dataSpans) {
70 evaluator_->setInput(item.first->GetName(), item.second, false);
71 }
72 }
73}
74
76 : RooAbsL(other),
77 apply_weight_squared(other.apply_weight_squared),
78 _first(other._first),
79 lastSection_(other.lastSection_),
80 cachedResult_(other.cachedResult_),
81 evaluator_(other.evaluator_)
82{
83 paramTracker_ = std::make_unique<RooChangeTracker>(*other.paramTracker_);
84}
85
87
88//////////////////////////////////////////////////////////////////////////////////
89
90/// Returns true if value was changed, false otherwise.
92{
93 if (apply_weight_squared != flag) {
95 return true;
96 }
97 // setValueDirty();
98 return false;
99}
100
101namespace {
102
103// For now, almost exact copy of RooNLLVar::computeScalarFunc.
104RooNLLVar::ComputeResult computeBatchFunc(std::span<const double> probas, RooAbsData *dataClone, bool weightSq,
105 std::size_t stepSize, std::size_t firstEvent, std::size_t lastEvent)
106{
109 RooNaNPacker packedNaN(0.f);
110
111 for (auto i = firstEvent; i < lastEvent; i += stepSize) {
112 dataClone->get(i);
113
114 double weight = dataClone->weight();
115
116 if (0. == weight * weight)
117 continue;
118 if (weightSq)
119 weight = dataClone->weightSquared();
120
121 double logProba = std::log(probas[i]);
122 const double term = -weight * logProba;
123
124 kahanWeight.Add(weight);
125 kahanProb.Add(term);
126 packedNaN.accumulate(term);
127 }
128
129 if (packedNaN.getPayload() != 0.) {
130 // Some events with evaluation errors. Return "badness" of errors.
131 return {ROOT::Math::KahanSum<double>{packedNaN.getNaNWithPayload()}, kahanWeight.Sum()};
132 }
133
134 return {kahanProb, kahanWeight.Sum()};
135}
136
137} // namespace
138
139//////////////////////////////////////////////////////////////////////////////////
140/// Calculate and return likelihood on subset of data from firstEvent to lastEvent
141/// processed with a step size of 'stepSize'. If this an extended likelihood and
142/// and the zero event is processed the extended term is added to the return
143/// likelihood.
144///
146RooUnbinnedL::evaluatePartition(Section events, std::size_t /*components_begin*/, std::size_t /*components_end*/)
147{
148 // Throughout the calculation, we use Kahan's algorithm for summing to
149 // prevent loss of precision - this is a factor four more expensive than
150 // straight addition, but since evaluating the PDF is usually much more
151 // expensive than that, we tolerate the additional cost...
153 double sumWeight;
154
155 // Do not reevaluate likelihood if parameters nor event range have changed
156 if (!paramTracker_->hasChanged(true) && events == lastSection_ &&
157 (cachedResult_.Sum() != 0 || cachedResult_.Carry() != 0))
158 return cachedResult_;
159
160 if (evaluator_) {
161 // Here, we have a memory allocation that should be avoided when this
162 // code needs to be optimized.
163 std::span<const double> probas = evaluator_->run();
164 std::tie(result, sumWeight) =
165 computeBatchFunc(probas, data_.get(), apply_weight_squared, 1, events.begin(N_events_), events.end(N_events_));
166 } else {
167 data_->store()->recalculateCache(nullptr, events.begin(N_events_), events.end(N_events_), 1, true);
168 std::tie(result, sumWeight) =
170 events.begin(N_events_), events.end(N_events_));
171 }
172
173 // include the extended maximum likelihood term, if requested
174 if (extended_ && events.begin_fraction == 0) {
175 result += pdf_->extendedTerm(*data_, apply_weight_squared);
176 }
177
178 // If part of simultaneous PDF normalize probability over
179 // number of simultaneous PDFs: -sum(log(p/n)) = -sum(log(p)) + N*log(n)
180 if (sim_count_ > 1) {
181 result += sumWeight * log(1.0 * sim_count_);
182 }
183
184 // At the end of the first full calculation, wire the caches. This doesn't
185 // need to be done in BatchMode with the RooFit driver.
186 if (_first && !evaluator_) {
187 _first = false;
188 pdf_->wireAllCaches();
189 }
190
192 lastSection_ = events;
193 return result;
194}
195
196} // namespace TestStatistics
197} // namespace RooFit
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...
RooAbsData is the common 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
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
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)
std::pair< ROOT::Math::KahanSum< double >, double > ComputeResult
Definition RooNLLVar.h:56
std::map< RooFit::Detail::DataKey, std::span< const double > > getDataSpans(RooAbsData const &data, std::string const &rangeName, RooSimultaneous const *simPdf, bool skipZeroWeights, bool takeGlobalObservablesFromData, std::stack< std::vector< double > > &buffers)
Extract all content from a RooFit datasets as a map of spans.
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...
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.