[RF] Translate RooStats tutorial rs101_limitexample.C to Python #10853

busorgin · 2022-06-29T14:02:08Z

This Pull request:

Changes or fixes:

Add rs101_limitexample.py - Converted from rs101_limitexample.C

Checklist:

tested changes locally
updated the docs (if necessary)

This PR fixes # Part of 8758

phsft-bot · 2022-06-29T14:02:13Z

Can one of the admins verify this patch?

guitargeek · 2022-06-29T17:42:48Z

tutorials/roostats/rs101_limitexample.py

+## \macro_output
+## \macro_code
+##
+## \author Kyle Cranmer


Suggested change

## \author Kyle Cranmer

## \date June 2022

## \authors Artem Busorgin, Kyle Cranmer (C++ version)

Please add also the date, and you can also add your name here to get some credits :)

guitargeek · 2022-06-29T17:44:08Z

tutorials/roostats/rs101_limitexample.py

+t.Stop()
+t.Print()
+
+dataCanvas.SaveAs("rs101_limitexample.png")


Suggested change

dataCanvas.SaveAs("rs101_limitexample.png")

dataCanvas.SaveAs("rs101_limitexample.png")

# TODO: The MCMCCalculator has to be destructed first. Otherwise, we can get

# segmentation faults depending on the destruction order, which is random in

# Python. Probably the issue is that some object has a non-owning pointer to

# another object, which it uses in its destructor. This should be fixed either

# in the design of RooStats in C++, or with phythonizations.

del mc

Can you please add this at the bottom please? By deleting mc first, we are avoiding some destruction order fiasco that occasionally causes crashes in this tutorial.

guitargeek · 2022-06-29T17:45:30Z

tutorials/roostats/rs101_limitexample.py

+# also plot the points in the markov chain
+chainData = mcInt.GetChainAsDataSet()
+
+assert chainData


Suggested change

assert chainData

No asserts in tutorials please (I know they were in the C++ version too, but asserts are only to check if things went well in debug builds or unit tests. For user code like these tutorials, it should never be necessary to use them)

guitargeek · 2022-06-29T17:45:53Z

tutorials/roostats/rs101_limitexample.py

+assert chainData
+print("plotting the chain data - nentries = ", chainData.numEntries())
+chain = ROOT.RooStats.GetAsTTree("chainTreeData", "chainTreeData", chainData)
+assert chain


Suggested change

assert chain

No asserts in tutorials

guitargeek · 2022-06-29T17:46:02Z

tutorials/roostats/rs101_limitexample.py

+
+# the points used in the profile construction
+parScanData = fc.GetPointsToScan()
+assert parScanData


Suggested change

assert parScanData

No asserts in tutorials

guitargeek · 2022-06-29T17:46:31Z

tutorials/roostats/rs101_limitexample.py

+    y = evt.getRealValue("ratioSigEff")
+    z = evt.getRealValue("s")
+    gr.SetPoint(ievt, x, y, z)
+    # print("{}  {}  {}  {}".format(ievt, x, y, z))


Suggested change

# print("{} {} {} {}".format(ievt, x, y, z))

It's not good to have commented out code in tutorials.

guitargeek · 2022-06-29T17:47:24Z

tutorials/roostats/rs101_limitexample.py

+data = ROOT.RooDataSet("exampleData", "exampleData", ROOT.RooArgSet(obs))
+data.add(obs)
+
+all = ROOT.RooArgSet(s, ratioBkgEff, ratioSigEff)


Suggested change

all = ROOT.RooArgSet(s, ratioBkgEff, ratioSigEff)

all is never used. You should delete this line

guitargeek · 2022-06-29T17:49:08Z

tutorials/roostats/rs101_limitexample.py

+modelWithConstraints = wspace.pdf("modelWithConstraints") # get the model
+obs = wspace.var("obs") # get the observable


Suggested change

modelWithConstraints = wspace.pdf("modelWithConstraints") # get the model

obs = wspace.var("obs") # get the observable

modelWithConstraints = wspace["modelWithConstraints"] # get the model

obs = wspace["obs"] # get the observable

In PyRoot, you can universally get anything from a RooWorkspace with the []-operator. You should use this everytime you see var(), pdf(), or function() in C++ (there are also more occurences of var() to replace in the lines below).

guitargeek · 2022-06-29T17:49:44Z

tutorials/roostats/rs101_limitexample.py

+
+ratioSigEff = wspace.var("ratioSigEff") # get uncertain parameter to constrain
+ratioBkgEff = wspace.var("ratioBkgEff") # get uncertain parameter to constrain
+constrainedParams = ROOT.RooArgSet(ratioSigEff, ratioBkgEff) # need to constrain these in the fit (should change default behavior)


Suggested change

constrainedParams = ROOT.RooArgSet(ratioSigEff, ratioBkgEff) # need to constrain these in the fit (should change default behavior)

constrainedParams = {ratioSigEff, ratioBkgEff} # need to constrain these in the fit (should change default behavior)

In PyROOT, you can directly use Python sets instead of RooArgSets

guitargeek · 2022-06-29T17:50:10Z

tutorials/roostats/rs101_limitexample.py

+
+# Create an example dataset with 160 observed events
+obs.setVal(160.)
+data = ROOT.RooDataSet("exampleData", "exampleData", ROOT.RooArgSet(obs))


Suggested change

data = ROOT.RooDataSet("exampleData", "exampleData", ROOT.RooArgSet(obs))

data = ROOT.RooDataSet("exampleData", "exampleData", {obs})

guitargeek · 2022-06-29T17:50:47Z

tutorials/roostats/rs101_limitexample.py

+all = ROOT.RooArgSet(s, ratioBkgEff, ratioSigEff)
+
+# not necessary
+modelWithConstraints.fitTo(data, ROOT.RooFit.Constrain(ROOT.RooArgSet(ratioSigEff, ratioBkgEff)))


Suggested change

modelWithConstraints.fitTo(data, ROOT.RooFit.Constrain(ROOT.RooArgSet(ratioSigEff, ratioBkgEff)))

modelWithConstraints.fitTo(data, Constrain=constrainedParams)

There was already a set defined for them before. You should reuse that here.

guitargeek · 2022-06-29T17:51:20Z

tutorials/roostats/rs101_limitexample.py

+modelConfig.SetParametersOfInterest(paramOfInterest)
+modelConfig.SetNuisanceParameters(constrainedParams)
+modelConfig.SetObservables(obs)
+modelConfig.SetGlobalObservables(ROOT.RooArgSet(gSigEff, gSigBkg))


Suggested change

modelConfig.SetGlobalObservables(ROOT.RooArgSet(gSigEff, gSigBkg))

modelConfig.SetGlobalObservables({gSigEff, gSigBkg})

guitargeek · 2022-06-29T17:52:06Z

tutorials/roostats/rs101_limitexample.py

+
+modelConfig = ROOT.RooStats.ModelConfig(wspace)
+modelConfig.SetPdf(modelWithConstraints)
+modelConfig.SetParametersOfInterest(paramOfInterest)


Suggested change

modelConfig.SetParametersOfInterest(paramOfInterest)

modelConfig.SetParametersOfInterest({s})

The paramOfInterest variable above is used nowhere else. You should better delete it, and directly use {s} here to be more explicit.

guitargeek · 2022-06-29T17:53:04Z

tutorials/roostats/rs101_limitexample.py

+wspace.writeToFile("rs101_ws.root")
+
+# First, let's use a Calculator based on the Profile Likelihood Ratio
+# plc = ROOT.RooStats.ProfileLikelihoodCalculator(data, modelWithConstraints, paramOfInterest)


Suggested change

# plc = ROOT.RooStats.ProfileLikelihoodCalculator(data, modelWithConstraints, paramOfInterest)

You should remove this commented out code, as there is not good reason to have it here

guitargeek

Hi @busorgin thank you so much for this contribution, which is really appreciated 👍 I which context are you doing this? Are you a RooFit user who wants to contribute, or you were just browsing GitHub for issues you could address?

Please address the requests I have (mostly about Pythonic code, and one fix for the destruction order fiasco).

Also, after applying all these suggestions, could you please format the new tutorial with black?

black --line-length=120 rs101_limitexample.py

busorgin · 2022-06-29T21:52:32Z

Hi @guitargeek, thanks for the review. I addressed the suggestions and formatted the code.

I was looking for an open source project to contribute to that would be interesting to me, settled on root

guitargeek · 2022-06-29T23:00:23Z

@phsft-bot build

phsft-bot · 2022-06-29T23:00:29Z

Starting build on ROOT-debian10-i386/soversion, ROOT-performance-centos8-multicore/cxx17, ROOT-ubuntu18.04/nortcxxmod, ROOT-ubuntu2004/python3, mac1015/cxx17, mac11/cxx14, windows10/cxx14
How to customize builds

guitargeek · 2022-06-29T23:03:22Z

Alright, thank you very much! That's very nice you chose ROOT to contribute to!

If you need any guidance on what you could help with, just ask me. Sometimes it's not easy to find the issues that are actually easy to work on, as the code is quite complex.

Are you looking only Python-related issues, or would you also be interested in doing some C++ contributions to beef up your C++ skills and experience?

As for this PR, I'll merge it if the CI bot tests pass.

busorgin · 2022-06-30T10:50:00Z

Thanks, I'll keep that in mind.

I was looking for both C++ and Python. I found #9859, seemed like a good starter item on the C++ side. Open to suggestions if there are other good issues to look into.

guitargeek · 2022-06-30T11:34:10Z

Hi @busorgin, okay sounds great! The issue you linked is indeed a very good one to get started with RooFit C++ development. If you have any questions about it, we can follow up in the issue thread there.

Of course, it would also be of great help if you continue doing these RooStats tutorial translations! People use more Python than C++ now, and a better coverage of the Python tutorials is something requested by many users, in particular for RooStats. So among all the "easier" RooFit issues, the RooStats tutorial translations have the highest priority.

guitargeek

LGTM, thank you for this contribution!

new tutorial translation runs on my local machine and on the CI
the output of the tutorial is the same as for the C++ version

Convert rs101_limitexample.C to Python

29ffad4

busorgin requested a review from couet as a code owner June 29, 2022 14:02

guitargeek self-requested a review June 29, 2022 17:39

guitargeek self-assigned this Jun 29, 2022

guitargeek added the in:RooFit/RooStats label Jun 29, 2022

guitargeek reviewed Jun 29, 2022

View reviewed changes

guitargeek requested changes Jun 29, 2022

View reviewed changes

guitargeek changed the title ~~Convert rs101_limitexample.C to Python~~ [RF] Convert rs101_limitexample.C to Python Jun 29, 2022

Fix code issues and format

f770604

guitargeek changed the title ~~[RF] Convert rs101_limitexample.C to Python~~ [RF] Translate RooStats tutorial rs101_limitexample.C to Python Jun 30, 2022

guitargeek approved these changes Jun 30, 2022

View reviewed changes

guitargeek merged commit 134e2a1 into root-project:master Jun 30, 2022

guitargeek added the new contributor label Jul 8, 2022

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[RF] Translate RooStats tutorial rs101_limitexample.C to Python #10853

[RF] Translate RooStats tutorial rs101_limitexample.C to Python #10853

busorgin commented Jun 29, 2022 •

edited

phsft-bot commented Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek Jun 29, 2022

guitargeek left a comment •

edited

busorgin commented Jun 29, 2022

guitargeek commented Jun 29, 2022

phsft-bot commented Jun 29, 2022

guitargeek commented Jun 29, 2022

busorgin commented Jun 30, 2022

guitargeek commented Jun 30, 2022

guitargeek left a comment

	## \author Kyle Cranmer
	## \date June 2022
	## \authors Artem Busorgin, Kyle Cranmer (C++ version)

-dataCanvas.SaveAs("rs101_limitexample.png")
+dataCanvas.SaveAs("rs101_limitexample.png")
+# TODO: The MCMCCalculator has to be destructed first. Otherwise, we can get
+# segmentation faults depending on the destruction order, which is random in
+# Python. Probably the issue is that some object has a non-owning pointer to
+# another object, which it uses in its destructor. This should be fixed either
+# in the design of RooStats in C++, or with phythonizations.
+del mc

		modelWithConstraints = wspace.pdf("modelWithConstraints") # get the model
		obs = wspace.var("obs") # get the observable

	constrainedParams = ROOT.RooArgSet(ratioSigEff, ratioBkgEff) # need to constrain these in the fit (should change default behavior)
	constrainedParams = {ratioSigEff, ratioBkgEff} # need to constrain these in the fit (should change default behavior)

	data = ROOT.RooDataSet("exampleData", "exampleData", ROOT.RooArgSet(obs))
	data = ROOT.RooDataSet("exampleData", "exampleData", {obs})

	modelWithConstraints.fitTo(data, ROOT.RooFit.Constrain(ROOT.RooArgSet(ratioSigEff, ratioBkgEff)))
	modelWithConstraints.fitTo(data, Constrain=constrainedParams)

	modelConfig.SetGlobalObservables(ROOT.RooArgSet(gSigEff, gSigBkg))
	modelConfig.SetGlobalObservables({gSigEff, gSigBkg})

	modelConfig.SetParametersOfInterest(paramOfInterest)
	modelConfig.SetParametersOfInterest({s})

[RF] Translate RooStats tutorial rs101_limitexample.C to Python #10853

[RF] Translate RooStats tutorial rs101_limitexample.C to Python #10853

Conversation

busorgin commented Jun 29, 2022 • edited

This Pull request:

Changes or fixes:

Checklist:

phsft-bot commented Jun 29, 2022

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

guitargeek left a comment • edited

Choose a reason for hiding this comment

busorgin commented Jun 29, 2022

guitargeek commented Jun 29, 2022

phsft-bot commented Jun 29, 2022

guitargeek commented Jun 29, 2022

busorgin commented Jun 30, 2022

guitargeek commented Jun 30, 2022

guitargeek left a comment

Choose a reason for hiding this comment

busorgin commented Jun 29, 2022 •

edited

guitargeek left a comment •

edited