Logo ROOT  
Reference Guide
RooConstraintSum.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooFitCore *
4 * @(#)root/roofitcore:$Id$
5 * Authors: *
6 * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7 * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/**
18\file RooConstraintSum.cxx
19\class RooConstraintSum
20\ingroup Roofitcore
21
22RooConstraintSum calculates the sum of the -(log) likelihoods of
23a set of RooAbsPfs that represent constraint functions. This class
24is used to calculate the composite -log(L) of constraints to be
25added to the regular -log(L) in RooAbsPdf::fitTo() with Constrain(..)
26arguments.
27**/
28
29
30#include "RooConstraintSum.h"
31#include "RooAbsData.h"
32#include "RooAbsReal.h"
33#include "RooAbsPdf.h"
34#include "RooErrorHandler.h"
35#include "RooArgSet.h"
36#include "RooMsgService.h"
37#include "RooHelpers.h"
38#include "RooWorkspace.h"
40
41#include <memory>
42
44
45
46////////////////////////////////////////////////////////////////////////////////
47/// Constructor with set of constraint p.d.f.s. All elements in constraintSet must inherit from RooAbsPdf.
48
49RooConstraintSum::RooConstraintSum(const char* name, const char* title, const RooArgSet& constraintSet, const RooArgSet& normSet, bool takeGlobalObservablesFromData) :
50 RooAbsReal(name, title),
51 _set1("set1","First set of components",this),
52 _takeGlobalObservablesFromData{takeGlobalObservablesFromData}
53{
54 for (const auto comp : constraintSet) {
55 if (!dynamic_cast<RooAbsPdf*>(comp)) {
56 coutE(InputArguments) << "RooConstraintSum::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
57 << " is not of type RooAbsPdf" << std::endl ;
59 }
60 _set1.add(*comp) ;
61 }
62
63 _paramSet.add(normSet) ;
64}
65
66
67////////////////////////////////////////////////////////////////////////////////
68/// Copy constructor.
69
71 RooAbsReal(other, name),
72 _set1("set1",this,other._set1),
73 _paramSet(other._paramSet),
74 _takeGlobalObservablesFromData{other._takeGlobalObservablesFromData}
75{
76}
77
78
79////////////////////////////////////////////////////////////////////////////////
80/// Return sum of -log of constraint p.d.f.s.
81
83{
84 double sum(0);
85
86 for (const auto comp : _set1) {
87 sum -= static_cast<RooAbsPdf*>(comp)->getLogVal(&_paramSet);
88 }
89
90 return sum;
91}
92
93
94void RooConstraintSum::computeBatch(cudaStream_t *, double *output, size_t /*size*/,
95 RooFit::Detail::DataMap const &dataMap) const
96{
97 double sum(0);
98
99 for (const auto comp : _set1) {
100 sum -= std::log(dataMap.at(comp)[0]);
101 }
102
103 output[0] = sum;
104}
105
106
107std::unique_ptr<RooArgSet> RooConstraintSum::fillNormSetForServer(RooArgSet const& /*normSet*/, RooAbsArg const& /*server*/) const {
108 return std::make_unique<RooArgSet>(_paramSet);
109}
110
111
112////////////////////////////////////////////////////////////////////////////////
113/// Replace the variables in this RooConstraintSum with the global observables
114/// in the dataset if they match by name. This function will do nothing if this
115/// RooConstraintSum is configured to not use the global observables stored in
116/// datasets.
117bool RooConstraintSum::setData(RooAbsData const& data, bool /*cloneData=true*/) {
118 if(_takeGlobalObservablesFromData && data.getGlobalObservables()) {
120 }
121 return true;
122}
123
124
125namespace {
126
127std::unique_ptr<RooArgSet> getGlobalObservables(
128 RooAbsPdf const& pdf, RooArgSet const* globalObservables, const char* globalObservablesTag)
129{
130
131 if(globalObservables && globalObservablesTag) {
132 // error!
133 std::string errMsg = "RooAbsPdf::fitTo: GlobalObservables and GlobalObservablesTag options mutually exclusive!";
134 oocoutE(&pdf, Minimization) << errMsg << std::endl;
135 throw std::invalid_argument(errMsg);
136 }
137 if(globalObservables) {
138 // pass-throught of global observables
139 return std::make_unique<RooArgSet>(*globalObservables);
140 }
141
142 if(globalObservablesTag) {
143 oocoutI(&pdf, Minimization) << "User-defined specification of global observables definition with tag named '"
144 << globalObservablesTag << "'" << std::endl;
145 } else {
146 // Neither GlobalObservables nor GlobalObservablesTag has been processed -
147 // try if a default tag is defined in the head node Check if head not
148 // specifies default global observable tag
149 if(auto defaultGlobalObservablesTag = pdf.getStringAttribute("DefaultGlobalObservablesTag")) {
150 oocoutI(&pdf, Minimization) << "p.d.f. provides built-in specification of global observables definition "
151 << "with tag named '" << defaultGlobalObservablesTag << "'" << std::endl;
152 globalObservablesTag = defaultGlobalObservablesTag;
153 }
154 }
155
156 if(globalObservablesTag) {
157 std::unique_ptr<RooArgSet> allVars{pdf.getVariables()} ;
158 return std::unique_ptr<RooArgSet>{static_cast<RooArgSet*>(allVars->selectByAttrib(globalObservablesTag, true))};
159 }
160
161 // no global observables specified
162 return nullptr;
163}
164
165
166RooArgSet const* tryToGetConstraintSetFromWorkspace(
167 RooAbsPdf const& pdf, RooWorkspace * workspace, std::string const& constraintSetCacheName) {
168 if(!workspace) return nullptr;
169
170 if(workspace->set(constraintSetCacheName.c_str())) {
171 // retrieve from cache
172 const RooArgSet *constr = workspace->set(constraintSetCacheName.c_str());
173 oocoutI(&pdf, Minimization)
174 << "createConstraintTerm picked up cached constraints from workspace with " << constr->size()
175 << " entries" << std::endl;
176 return constr;
177 }
178 return nullptr;
179}
180
181
182} // namespace
183
184
185////////////////////////////////////////////////////////////////////////////////
186/// Create the parameter constraint sum to add to the negative log-likelihood.
187/// \return a `nullptr` if the parameters are unconstrained.
188/// \param[in] name Name of the created RooConstraintSum object.
189/// \param[in] pdf The pdf model whose parameters should be constrained.
190/// Constraint terms will be extracted from RooProdPdf instances
191/// that are servers of the pdf (internal constraints).
192/// \param[in] data Dataset used in the fit with the constraint sum. It is
193/// used to figure out which are the observables and also to get the
194/// global observables definition and values if they are stored in
195/// the dataset.
196/// \param[in] constrainedParameters Set of parameters to constrain. If `nullptr`, all
197/// parameters will be considered.
198/// \param[in] externalConstraints Set of constraint terms that are not
199/// embedded in the pdf (external constraints).
200/// \param[in] globalObservables The normalization set for the constraint terms.
201/// If it is `nullptr`, the set of all constrained parameters will
202/// be used as the normalization set.
203/// \param[in] globalObservablesTag Alternative to define the normalization set
204/// for the constraint terms. All constrained parameters that have
205/// the attribute with the tag defined by `globalObservablesTag` are
206/// used. The `globalObservables` and `globalObservablesTag`
207/// parameters are mutually exclusive, meaning at least one of them
208/// has to be `nullptr`.
209/// \param[in] takeGlobalObservablesFromData If the dataset should be used to automatically
210/// define the set of global observables. If this is the case and the
211/// set of global observables is still defined manually with the
212/// `globalObservables` or `globalObservablesTag` parameters, the
213/// values of all global observables that are not stored in the
214/// dataset are taken from the model.
215/// \param[in] cloneConstraints true to clone constraints
216/// \param[in] workspace RooWorkspace to cache the set of constraints.
217std::unique_ptr<RooAbsReal> RooConstraintSum::createConstraintTerm(
218 std::string const& name,
219 RooAbsPdf const& pdf,
220 RooAbsData const& data,
221 RooArgSet const* constrainedParameters,
222 RooArgSet const* externalConstraints,
223 RooArgSet const* globalObservables,
224 const char* globalObservablesTag,
225 bool takeGlobalObservablesFromData,
226 bool cloneConstraints,
227 RooWorkspace * workspace)
228{
229 RooArgSet const& observables = *data.get();
230
231 bool doStripDisconnected = false ;
232
233 // If no explicit list of parameters to be constrained is specified apply default algorithm
234 // All terms of RooProdPdfs that do not contain observables and share a parameters with one or more
235 // terms that do contain observables are added as constrainedParameters.
236 RooArgSet cPars;
237 if(constrainedParameters) {
238 cPars.add(*constrainedParameters);
239 } else {
240 pdf.getParameters(&observables,cPars,false);
241 doStripDisconnected = true;
242 }
243
244 // Collect internal and external constraint specifications
245 RooArgSet allConstraints ;
246
247 auto observableNames = RooHelpers::getColonSeparatedNameString(observables);
248 auto constraintSetCacheName = std::string("CACHE_CONSTR_OF_PDF_") + pdf.GetName() + "_FOR_OBS_" + observableNames;
249
250 if (RooArgSet const* constr = tryToGetConstraintSetFromWorkspace(pdf, workspace, constraintSetCacheName)) {
251 allConstraints.add(*constr);
252 } else {
253
254 if (!cPars.empty()) {
255 std::unique_ptr<RooArgSet> internalConstraints{pdf.getAllConstraints(observables, cPars, doStripDisconnected)};
256 allConstraints.add(*internalConstraints);
257 }
258 if (externalConstraints) {
259 allConstraints.add(*externalConstraints);
260 }
261
262 // write to cache
263 if (workspace) {
264 oocoutI(&pdf, Minimization)
265 << "createConstraintTerm: caching constraint set under name "
266 << constraintSetCacheName << " with " << allConstraints.size() << " entries" << std::endl;
267 workspace->defineSetInternal(constraintSetCacheName.c_str(), allConstraints);
268 }
269 }
270
271 if (!allConstraints.empty()) {
272
273 oocoutI(&pdf, Minimization) << " Including the following constraint terms in minimization: "
274 << allConstraints << std::endl ;
275
276 // Identify global observables in the model.
277 auto glObs = getGlobalObservables(pdf, globalObservables, globalObservablesTag);
278 if(data.getGlobalObservables() && takeGlobalObservablesFromData) {
279 if(!glObs) {
280 // There were no global observables specified, but there are some in the
281 // dataset. We will just take them from the dataset.
282 oocoutI(&pdf, Minimization)
283 << "The following global observables have been automatically defined according to the dataset "
284 << "which also provides their values: " << *data.getGlobalObservables() << std::endl;
285 glObs = std::make_unique<RooArgSet>(*data.getGlobalObservables());
286 } else {
287 // There are global observables specified by the user and also some in
288 // the dataset.
289 RooArgSet globalsFromDataset;
290 data.getGlobalObservables()->selectCommon(*glObs, globalsFromDataset);
291 oocoutI(&pdf, Minimization)
292 << "The following global observables have been defined: " << *glObs
293 << "," << " with the values of " << globalsFromDataset
294 << " obtained from the dataset and the other values from the model." << std::endl;
295 }
296 } else if(glObs) {
297 oocoutI(&pdf, Minimization)
298 << "The following global observables have been defined and their values are taken from the model: "
299 << *glObs << std::endl;
300 // in this case we don;t take global observables from data
301 takeGlobalObservablesFromData = false;
302 } else {
303 if (!glObs)
304 oocoutI(&pdf, Minimization)
305 << "The global observables are not defined , normalize constraints with respect to the parameters " << cPars
306 << std::endl;
307 takeGlobalObservablesFromData = false;
308 }
309
310 // The constraint terms need to be cloned, because the global observables
311 // might be changed to have the same values as stored in data.
312 if (cloneConstraints) {
313 RooConstraintSum constraintTerm{name.c_str(), "nllCons", allConstraints, glObs ? *glObs : cPars,
314 takeGlobalObservablesFromData};
315 std::unique_ptr<RooAbsReal> constraintTermClone{static_cast<RooAbsReal *>(constraintTerm.cloneTree())};
316
317 // The parameters that are not connected to global observables from data
318 // need to be redirected to the original args to get the changes made by
319 // the minimizer. This excludes the global observables, where we take the
320 // clones with the values set to the values from the dataset if available.
321 RooArgSet allOriginalParams;
322 constraintTerm.getParameters(nullptr, allOriginalParams);
323 constraintTermClone->recursiveRedirectServers(allOriginalParams);
324
325 // Redirect the global observables to the ones from the dataset if applicable.
326 static_cast<RooConstraintSum *>(constraintTermClone.get())->setData(data, false);
327
328 // The computation graph for the constraints is very small, no need to do
329 // the tracking of clean and dirty nodes here.
330 constraintTermClone->setOperMode(RooAbsArg::ADirty);
331
332 return constraintTermClone;
333 }
334 // case we do not clone constraints (e.g. when using new Driver)
335 else {
336 // when we don't clone we need that global observables are not from data
337 if (takeGlobalObservablesFromData) {
338 oocoutE(&pdf, InputArguments) << "RooAbsPdf::Fit: Batch mode does not support yet GlobalObservable from data, use model as GlobalObservableSource " << std::endl;
339 throw std::invalid_argument("Invalid arguments for GlobalObservables in batch mode");
340 }
341 std::unique_ptr<RooAbsReal> constraintTerm(
342 new RooConstraintSum(name.c_str(), "nllCons", allConstraints, glObs ? *glObs : cPars, false));
343 return constraintTerm;
344 }
345 }
346
347 // no constraints
348 return nullptr;
349}
#define oocoutE(o, a)
Definition: RooMsgService.h:52
#define oocoutI(o, a)
Definition: RooMsgService.h:49
#define coutE(a)
Definition: RooMsgService.h:37
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
char name[80]
Definition: TGX11.cxx:110
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:71
bool recursiveRedirectServers(const RooAbsCollection &newServerList, bool mustReplaceAll=false, bool nameChange=false, bool recurseInNewSet=true)
Recursively replace all servers with the new servers in newSet.
Definition: RooAbsArg.cxx:1165
const Text_t * getStringAttribute(const Text_t *key) const
Get string attribute mapped under key 'key'.
Definition: RooAbsArg.cxx:299
RooArgSet * getVariables(bool stripDisconnected=true) const
Return RooArgSet with all variables (tree leaf nodes of expresssion tree)
Definition: RooAbsArg.cxx:2035
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...
Definition: RooAbsArg.cxx:541
bool empty() const
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:62
RooArgSet const * getGlobalObservables() const
Returns snapshot of global observables stored in this data.
Definition: RooAbsData.h:304
virtual RooArgSet * getAllConstraints(const RooArgSet &observables, RooArgSet &constrainedParams, bool stripDisconnected=true) const
This helper function finds and collects all constraints terms of all component p.d....
Definition: RooAbsPdf.cxx:3374
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:56
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
RooConstraintSum calculates the sum of the -(log) likelihoods of a set of RooAbsPfs that represent co...
bool setData(RooAbsData const &data, bool cloneData=true)
Replace the variables in this RooConstraintSum with the global observables in the dataset if they mat...
RooListProxy _set1
Set of constraint terms.
double evaluate() const override
Return sum of -log of constraint p.d.f.s.
void computeBatch(cudaStream_t *, double *output, size_t size, RooFit::Detail::DataMap const &) const override
Base function for computing multiple values of a RooAbsReal.
RooArgSet _paramSet
Set of parameters to which constraints apply.
static std::unique_ptr< RooAbsReal > createConstraintTerm(std::string const &name, RooAbsPdf const &pdf, RooAbsData const &data, RooArgSet const *constrainedParameters, RooArgSet const *externalConstraints, RooArgSet const *globalObservables, const char *globalObservablesTag, bool takeGlobalObservablesFromData, bool cloneConstraints, RooWorkspace *workspace)
Create the parameter constraint sum to add to the negative log-likelihood.
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...
const bool _takeGlobalObservablesFromData
If the global observable values are taken from data.
static void softAbort()
Soft abort function that interrupts macro execution but doesn't kill ROOT.
auto & at(RooAbsArg const *arg, RooAbsArg const *=nullptr)
Definition: DataMap.h:88
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
bool defineSetInternal(const char *name, const RooArgSet &aset)
const RooArgSet * set(const char *name)
Return pointer to previously defined named set with given nmame If no such set is found a null pointe...
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1748
@ Minimization
Definition: RooGlobalFunc.h:63
@ InputArguments
Definition: RooGlobalFunc.h:64
std::string getColonSeparatedNameString(RooArgSet const &argSet)
Create a string with all sorted names of RooArgSet elements separated by colons.
Definition: RooHelpers.cxx:252
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345
static void output()