Logo ROOT  
Reference Guide
RooNLLVarNew.cxx
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 * Authors:
4 * Jonas Rembser, CERN 2021
5 * Emmanouil Michalainas, CERN 2021
6 *
7 * Copyright (c) 2021, CERN
8 *
9 * Redistribution and use in source and binary forms,
10 * with or without modification, are permitted according to the terms
11 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
12 */
13
14/**
15\file RooNLLVarNew.cxx
16\class RooNLLVarNew
17\ingroup Roofitcore
18
19This is a simple class designed to produce the nll values needed by the fitter.
20In contrast to the `RooNLLVar` class, any logic except the bare minimum has been
21transfered away to other classes, like the `RooFitDriver`. This class also calls
22functions from `RooBatchCompute` library to provide faster computation times.
23**/
24
25#include <RooNLLVarNew.h>
26
27#include <RooAddition.h>
28#include <RooFormulaVar.h>
29#include <RooNaNPacker.h>
30#include <RooRealVar.h>
32
33#include <ROOT/StringUtils.hxx>
34
35#include <TClass.h>
36#include <TMath.h>
37#include <Math/Util.h>
38#include <TMath.h>
39
40#include <numeric>
41#include <stdexcept>
42#include <vector>
43
44using namespace ROOT::Experimental;
45
46// Declare constexpr static members to make them available if odr-used in C++14.
47constexpr const char *RooNLLVarNew::weightVarName;
48constexpr const char *RooNLLVarNew::weightVarNameSumW2;
49
50namespace {
51
52template <class Input>
53double kahanSum(Input const &input)
54{
55 return ROOT::Math::KahanSum<double, 4u>::Accumulate(input.begin(), input.end()).Sum();
56}
57
58RooArgSet getObs(RooAbsArg const &arg, RooArgSet const &observables)
59{
60 RooArgSet out;
61 arg.getObservables(&observables, out);
62 return out;
63}
64
65} // namespace
66
67/** Construct a RooNLLVarNew
68\param name the name
69\param title the title
70\param pdf The pdf for which the nll is computed for
71\param observables The observabes of the pdf
72\param isExtended Set to true if this is an extended fit
73**/
74RooNLLVarNew::RooNLLVarNew(const char *name, const char *title, RooAbsPdf &pdf, RooArgSet const &observables,
75 bool isExtended, bool doOffset, bool binnedL)
76 : RooAbsReal(name, title), _pdf{"pdf", "pdf", this, pdf}, _observables{getObs(pdf, observables)},
77 _isExtended{isExtended}, _binnedL{binnedL},
78 _weightVar{"weightVar", "weightVar", this, *new RooRealVar(weightVarName, weightVarName, 1.0), true, false, true},
79 _weightSquaredVar{weightVarNameSumW2,
80 weightVarNameSumW2,
81 this,
82 *new RooRealVar("weightSquardVar", "weightSquaredVar", 1.0),
83 true,
84 false,
85 true}
86{
87 if (_binnedL) {
88 if (_observables.size() != 1) {
89 throw std::runtime_error("BinnedPdf optimization only works with a 1D pdf.");
90 } else {
91 auto *var = static_cast<RooRealVar *>(_observables.first());
92 std::list<double> *boundaries = pdf.binBoundaries(*var, var->getMin(), var->getMax());
93 std::list<double>::iterator biter = boundaries->begin();
94 _binw.resize(boundaries->size() - 1);
95 double lastBound = (*biter);
96 ++biter;
97 int ibin = 0;
98 while (biter != boundaries->end()) {
99 _binw[ibin] = (*biter) - lastBound;
100 lastBound = (*biter);
101 ibin++;
102 ++biter;
103 }
104 }
105 }
106
108 enableOffsetting(doOffset);
109}
110
112 : RooAbsReal(other, name), _pdf{"pdf", this, other._pdf}, _observables{other._observables},
113 _isExtended{other._isExtended}, _weightSquared{other._weightSquared}, _binnedL{other._binnedL},
114 _doOffset{other._doOffset}, _simCount{other._simCount}, _prefix{other._prefix},
115 _weightVar{"weightVar", this, other._weightVar}, _weightSquaredVar{"weightSquaredVar", this,
116 other._weightSquaredVar}
117{
118}
119
120/** Compute multiple negative logs of propabilities
121
122\param output An array of doubles where the computation results will be stored
123\param nOut not used
124\note nEvents is the number of events to be processed (the dataMap size)
125\param dataMap A map containing spans with the input data for the computation
126**/
127void RooNLLVarNew::computeBatch(cudaStream_t * /*stream*/, double *output, size_t /*nOut*/,
128 RooFit::Detail::DataMap const &dataMap) const
129{
130 std::size_t nEvents = dataMap.at(_pdf).size();
131
132 auto weights = dataMap.at(_weightVar);
133 auto weightsSumW2 = dataMap.at(_weightSquaredVar);
134 auto weightSpan = _weightSquared ? weightsSumW2 : weights;
135
136 if (_binnedL) {
138 ROOT::Math::KahanSum<double> sumWeightKahanSum{0.0};
139 auto preds = dataMap.at(&*_pdf);
140
141 for (std::size_t i = 0; i < nEvents; ++i) {
142
143 double eventWeight = weightSpan[i];
144
145 // Calculate log(Poisson(N|mu) for this bin
146 double N = eventWeight;
147 double mu = preds[i] * _binw[i];
148
149 if (mu <= 0 && N > 0) {
150
151 // Catch error condition: data present where zero events are predicted
152 logEvalError(Form("Observed %f events in bin %lu with zero event yield", N, (unsigned long)i));
153
154 } else if (std::abs(mu) < 1e-10 && std::abs(N) < 1e-10) {
155
156 // Special handling of this case since log(Poisson(0,0)=0 but can't be calculated with usual log-formula
157 // since log(mu)=0. No update of result is required since term=0.
158
159 } else {
160
161 result += -1 * (-mu + N * log(mu) - TMath::LnGamma(N + 1));
162 sumWeightKahanSum += eventWeight;
163 }
164 }
165
166 output[0] = finalizeResult(std::move(result), sumWeightKahanSum.Sum());
167
168 return;
169 }
170
171 auto probas = dataMap.at(_pdf);
172
173 _logProbasBuffer.resize(nEvents);
174 (*_pdf).getLogProbabilities(probas, _logProbasBuffer.data());
175
176 _sumWeight = weights.size() == 1 ? weights[0] * nEvents : kahanSum(weights);
177
178 if (_isExtended && _weightSquared && _sumWeight2 == 0.0) {
179 _sumWeight2 = weights.size() == 1 ? weightsSumW2[0] * nEvents : kahanSum(weightsSumW2);
180 }
181
183 RooNaNPacker packedNaN(0.f);
184
185 for (std::size_t i = 0; i < nEvents; ++i) {
186
187 double eventWeight = weightSpan.size() > 1 ? weightSpan[i] : weightSpan[0];
188 if (0. == eventWeight * eventWeight)
189 continue;
190
191 const double term = -eventWeight * _logProbasBuffer[i];
192
193 kahanProb.Add(term);
194 packedNaN.accumulate(term);
195 }
196
197 if (packedNaN.getPayload() != 0.) {
198 // Some events with evaluation errors. Return "badness" of errors.
199 kahanProb = packedNaN.getNaNWithPayload();
200 }
201
202 if (_isExtended) {
203 double expected = _pdf->expectedEvents(&_observables);
204 kahanProb += _pdf->extendedTerm(_sumWeight, expected, _weightSquared ? _sumWeight2 : 0.0);
205 }
206
207 output[0] = finalizeResult(std::move(kahanProb), _sumWeight);
208}
209
210void RooNLLVarNew::getParametersHook(const RooArgSet * /*nset*/, RooArgSet *params, bool /*stripDisconnected*/) const
211{
212 // strip away the observables and weights
213 params->remove(_observables, true, true);
214 params->remove(RooArgList{*_weightVar, *_weightSquaredVar}, true, true);
215}
216
217////////////////////////////////////////////////////////////////////////////////
218/// Clones the PDF recursively and prefixes the names of all nodes, except for
219/// parameter nodes. Used for simultaneous fits.
220/// \return A RooArgSet with the new observable args.
221/// \param[in] prefix The prefix to add to the observables and weight names.
222RooArgSet RooNLLVarNew::prefixArgNames(std::string const &prefix)
223{
224 _prefix = prefix;
225
226 std::unique_ptr<RooAbsReal> pdfClone = RooHelpers::cloneTreeWithSameParameters(*_pdf, &_observables);
227
228 redirectServers(RooArgList{*pdfClone});
229
230 RooArgSet parameters;
231 pdfClone->getParameters(&_observables, parameters);
232
234
235 RooArgSet nodes;
236 pdfClone->treeNodeServerList(&nodes);
237 for (RooAbsArg *arg : nodes) {
238 if (!parameters.find(*arg)) {
239 arg->SetName((prefix + arg->GetName()).c_str());
240 if (dynamic_cast<RooRealVar *>(arg)) {
241 // It's an observable
242 static_cast<RooRealVar *>(arg)->setConstant();
243 _observables.add(*arg);
244 arg->setAttribute("__obs__");
245 }
246 }
247 }
248
249 addOwnedComponents(std::move(pdfClone));
250
252
253 return _observables;
254}
255
257{
260}
261
262////////////////////////////////////////////////////////////////////////////////
263/// Toggles the weight square correction.
265{
266 _weightSquared = flag;
267}
268
269std::unique_ptr<RooArgSet>
270RooNLLVarNew::fillNormSetForServer(RooArgSet const & /*normSet*/, RooAbsArg const & /*server*/) const
271{
272 if (_binnedL) {
273 return std::make_unique<RooArgSet>();
274 }
275 return nullptr;
276}
277
279{
280 _doOffset = flag;
281 _offset = {};
282}
283
285{
286 // If part of simultaneous PDF normalize probability over
287 // number of simultaneous PDFs: -sum(log(p/n)) = -sum(log(p)) + N*log(n)
288 if (_simCount > 1) {
289 result += weightSum * std::log(static_cast<double>(_simCount));
290 }
291
292 // Check if value offset flag is set.
293 if (_doOffset) {
294
295 // If no offset is stored enable this feature now
296 if (_offset == 0 && result != 0) {
297 _offset = result;
298 }
299
300 // Subtract offset
301 if (!RooAbsReal::hideOffset()) {
302 result -= _offset;
303 }
304 }
305 return result.Sum();
306}
#define e(i)
Definition: RSha256.hxx:103
#define N
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void input
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
char name[80]
Definition: TGX11.cxx:110
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2456
RooTemplateProxy< RooAbsReal > _weightVar
Definition: RooNLLVarNew.h:75
double finalizeResult(ROOT::Math::KahanSum< double > &&result, double weightSum) const
void enableOffsetting(bool) override
void computeBatch(cudaStream_t *, double *output, size_t nOut, RooFit::Detail::DataMap const &) const override
Compute multiple negative logs of propabilities.
void applyWeightSquared(bool flag) override
Toggles the weight square correction.
void getParametersHook(const RooArgSet *nset, RooArgSet *list, bool stripDisconnected) const override
ROOT::Math::KahanSum< double > _offset
! Offset as KahanSum to avoid loss of precision
Definition: RooNLLVarNew.h:79
RooTemplateProxy< RooAbsPdf > _pdf
Definition: RooNLLVarNew.h:65
std::vector< double > _logProbasBuffer
!
Definition: RooNLLVarNew.h:78
RooTemplateProxy< RooAbsReal > _weightSquaredVar
Definition: RooNLLVarNew.h:76
RooArgSet prefixArgNames(std::string const &prefix)
Clones the PDF recursively and prefixes the names of all nodes, except for parameter nodes.
std::vector< double > _binw
!
Definition: RooNLLVarNew.h:77
static constexpr const char * weightVarName
Definition: RooNLLVarNew.h:32
std::unique_ptr< RooArgSet > fillNormSetForServer(RooArgSet const &normSet, RooAbsArg const &server) const override
Fills a RooArgSet to be used as the normalization set for a server, given a normalization set for thi...
static constexpr const char * weightVarNameSumW2
Definition: RooNLLVarNew.h:33
static KahanSum< T, N > Accumulate(Iterator begin, Iterator end, T initialValue=T{})
Iterate over a range and return an instance of a KahanSum.
Definition: Util.h:211
void Add(T x)
Single-element accumulation. Will not vectorise.
Definition: Util.h:165
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:71
RooArgSet * getObservables(const RooArgSet &set, bool valueOnly=true) const
Given a set of possible observables, return the observables that this PDF depends on.
Definition: RooAbsArg.h:293
void SetName(const char *name) override
Set the name of the TNamed.
Definition: RooAbsArg.cxx:2314
bool addOwnedComponents(const RooAbsCollection &comps)
Take ownership of the contents of 'comps'.
Definition: RooAbsArg.cxx:2185
bool redirectServers(const RooAbsCollection &newServerList, bool mustReplaceAll=false, bool nameChange=false, bool isRecursionStep=false)
Replace all direct servers of this object with the new servers in newServerList.
Definition: RooAbsArg.cxx:999
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * first() const
void clear()
Clear contents. If the collection is owning, it will also delete the contents.
virtual double expectedEvents(const RooArgSet *nset) const
Return expected number of events to be used in calculation of extended likelihood.
Definition: RooAbsPdf.cxx:3232
double extendedTerm(double sumEntries, double expected, double sumEntriesW2=0.0) const
Definition: RooAbsPdf.cxx:826
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
virtual std::list< double > * binBoundaries(RooAbsRealLValue &obs, double xlo, double xhi) const
Retrieve bin boundaries if this distribution is binned in obs.
static bool hideOffset()
Definition: RooAbsReal.cxx:117
void logEvalError(const char *message, const char *serverValueString=nullptr) const
Log evaluation error message.
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
auto & at(RooAbsArg const *arg, RooAbsArg const *=nullptr)
Definition: DataMap.h:88
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
RVec< PromoteType< T > > abs(const RVec< T > &v)
Definition: RVec.hxx:1778
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1787
std::unique_ptr< T > cloneTreeWithSameParameters(T const &arg, RooArgSet const *observables=nullptr)
Clone RooAbsArg object and reattach to original parameters.
Definition: RooHelpers.h:145
void probas(TString dataset, TString fin="TMVA.root", Bool_t useTMVAStyle=kTRUE)
Double_t LnGamma(Double_t z)
Computation of ln[gamma(z)] for all z.
Definition: TMath.cxx:509
Little struct that can pack a float into the unused bits of the mantissa of a NaN double.
Definition: RooNaNPacker.h:28
float getPayload() const
Retrieve packed float.
Definition: RooNaNPacker.h:85
double getNaNWithPayload() const
Retrieve a NaN with the current float payload packed into the mantissa.
Definition: RooNaNPacker.h:90
void accumulate(double val)
Accumulate a packed float from another NaN into this.
Definition: RooNaNPacker.h:57
static void output()