Skip to content

Commit 030f599

Browse files
akeemlhguitargeek
akeemlh
authored and
guitargeek
committed
Translated StandardProfileInspectorDemo and StandardProfileLikelihoodDemo to pyROOT
1 parent 84b5b3a commit 030f599

File tree

4 files changed

+248
-12
lines changed

4 files changed

+248
-12
lines changed

‎tutorials/roostats/StandardProfileInspectorDemo.C‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -72,7 +72,7 @@ void StandardProfileInspectorDemo(const char *infile = "", const char *workspace
7272
// Try to open the file
7373
TFile *file = TFile::Open(filename);
7474

75-
// if input file was specified byt not found, quit
75+
// if input file was specified but not found, quit
7676
if (!file) {
7777
cout << "StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
7878
return;

‎tutorials/roostats/StandardProfileInspectorDemo.py‎

Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,103 @@
1+
## \file
2+
## \ingroup tutorial_roostats
3+
## \notebook -js
4+
## Standard demo of the ProfileInspector class
5+
## StandardProfileInspectorDemo
6+
##
7+
## This is a standard demo that can be used with any ROOT file
8+
## prepared in the standard way. You specify:
9+
## - name for input ROOT file
10+
## - name of workspace inside ROOT file that holds model and data
11+
## - name of ModelConfig that specifies details for calculator tools
12+
## - name of dataset
13+
##
14+
## With the values provided below this script will attempt to run the
15+
## standard hist2workspace example and read the ROOT file
16+
## that it produces.
17+
##
18+
## The actual heart of the demo is only about 10 lines long.
19+
##
20+
## The ProfileInspector plots the conditional maximum likelihood estimate
21+
## of each nuisance parameter in the model vs. the parameter of interest.
22+
## (aka. profiled value of nuisance parameter vs. parameter of interest)
23+
## (aka. best fit nuisance parameter with p.o.i fixed vs. parameter of interest)
24+
##
25+
## \macro_image
26+
## \macro_output
27+
## \macro_code
28+
##
29+
## \authors Akeem Hart, Kyle Cranmer (C++ Version)
30+
31+
import ROOT
32+
import sys
33+
import math
34+
35+
# -------------------------------------------------------
36+
# First part is just to access a user-defined file
37+
# or create the standard example file if it doesn't exist
38+
39+
workspaceName = "combined"
40+
modelConfigName = "ModelConfig"
41+
dataName = "obsData"
42+
filename = "results/example_combined_GaussExample_model.root"
43+
# if file does not exists generate with histfactory
44+
if ROOT.gSystem.AccessPathName(filename):
45+
# Normally this would be run on the command line
46+
print("will run standard hist2workspace example")
47+
ROOT.gROOT.ProcessLine(".! prepareHistFactory .")
48+
ROOT.gROOT.ProcessLine(".! hist2workspace config/example.xml")
49+
print("\n\n---------------------")
50+
print("Done creating example input")
51+
print("---------------------\n\n")
52+
53+
# Try to open the file
54+
try:
55+
file = ROOT.TFile.Open(filename)
56+
except:
57+
# if input file was specified but not found, quit
58+
print("StandardRooStatsDemoMacro: Input file %s is not found" % filename)
59+
sys.exit()
60+
61+
# -------------------------------------------------------
62+
# Tutorial starts here
63+
# -------------------------------------------------------
64+
65+
# get the workspace out of the file
66+
67+
w = file.Get(workspaceName)
68+
if not w:
69+
print("Workspace not found")
70+
sys.exit()
71+
72+
# get the modelConfig out of the file
73+
mc = w.obj(modelConfigName)
74+
75+
# get the modelConfig out of the file
76+
data = w.data(dataName)
77+
78+
# make sure ingredients are found
79+
if not data or not mc:
80+
w.Print()
81+
print("data or ModelConfig was not found")
82+
sys.exit()
83+
84+
# -----------------------------
85+
# now use the profile inspector
86+
p = ROOT.RooStats.ProfileInspector()
87+
list = p.GetListOfProfilePlots(data, mc)
88+
89+
# now make plots
90+
c1 = ROOT.TCanvas("c1", "ProfileInspectorDemo", 800, 200)
91+
n = list.GetSize()
92+
if n > 4:
93+
nx = int(math.sqrt(n))
94+
ny = ROOT.TMath.CeilNint(n / nx)
95+
nx = ROOT.TMath.CeilNint(math.sqrt(n))
96+
c1.Divide(ny, nx)
97+
else:
98+
c1.Divide(n)
99+
for i in range(n):
100+
c1.cd(i + 1)
101+
list.At(i).Draw("al")
102+
103+
c1.Update()

‎tutorials/roostats/StandardProfileLikelihoodDemo.py‎

Lines changed: 132 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,132 @@
1+
## \file
2+
## \ingroup tutorial_roostats
3+
## \notebook
4+
## Standard demo of the Profile Likelihood calculator
5+
## StandardProfileLikelihoodDemo
6+
##
7+
## This is a standard demo that can be used with any ROOT file
8+
## prepared in the standard way. You specify:
9+
## - name for input ROOT file
10+
## - name of workspace inside ROOT file that holds model and data
11+
## - name of ModelConfig that specifies details for calculator tools
12+
## - name of dataset
13+
## With the values provided below the macro will attempt to run the
14+
## standard hist2workspace example and read the ROOT file
15+
## that it produces.
16+
##
17+
## The actual heart of the demo is only about 10 lines long.
18+
##
19+
## The ProfileLikelihoodCalculator is based on Wilks's theorem
20+
## and the asymptotic properties of the profile likelihood ratio
21+
## (eg. that it is chi-square distributed for the true value).
22+
##
23+
## \macro_image
24+
## \macro_output
25+
## \macro_code
26+
##
27+
## \authors Akeem Hart, Kyle Cranmer (C++ Version)
28+
29+
import ROOT
30+
import sys
31+
32+
workspaceName = "combined"
33+
modelConfigName = "ModelConfig"
34+
dataName = "obsData"
35+
confLevel = 0.95
36+
nScanPoints = 50
37+
plotAsTF1 = False
38+
poiXMin = 1
39+
poiXMax = 0
40+
doHypoTest = False
41+
nullParamValue = 0
42+
filename = "results/example_combined_GaussExample_model.root"
43+
# if file does not exists generate with histfactory
44+
if ROOT.gSystem.AccessPathName(filename):
45+
# Normally this would be run on the command line
46+
print("will run standard hist2workspace example")
47+
ROOT.gROOT.ProcessLine(".! prepareHistFactory .")
48+
ROOT.gROOT.ProcessLine(".! hist2workspace config/example.xml")
49+
print("\n\n---------------------")
50+
print("Done creating example input")
51+
print("---------------------\n\n")
52+
53+
# Try to open the file
54+
try:
55+
file = ROOT.TFile.Open(filename)
56+
except:
57+
# if input file was specified but not found, quit
58+
print("StandardRooStatsDemoMacro: Input file %s is not found" % filename)
59+
sys.exit()
60+
61+
# -------------------------------------------------------
62+
# Tutorial starts here
63+
# -------------------------------------------------------
64+
65+
# get the workspace out of the file
66+
67+
w = file.Get(workspaceName)
68+
if not w:
69+
print("Workspace not found")
70+
sys.exit()
71+
72+
# get the modelConfig out of the file
73+
mc = w.obj(modelConfigName)
74+
75+
# get the modelConfig out of the file
76+
data = w.data(dataName)
77+
78+
# make sure ingredients are found
79+
if not data or not mc:
80+
w.Print()
81+
print("data or ModelConfig was not found")
82+
sys.exit()
83+
84+
# ---------------------------------------------
85+
# create and use the ProfileLikelihoodCalculator
86+
# to find and plot the 95% confidence interval
87+
# on the parameter of interest as specified
88+
# in the model config
89+
90+
pl = ROOT.RooStats.ProfileLikelihoodCalculator(data, mc)
91+
pl.SetConfidenceLevel(confLevel)
92+
interval = pl.GetInterval()
93+
94+
# print out the interval on the first Parameter of Interest
95+
firstPOI = mc.GetParametersOfInterest().first()
96+
print(
97+
"\n>>>> RESULT : "
98+
+ str(confLevel * 100)
99+
+ "% interval on "
100+
+ firstPOI.GetName()
101+
+ " is : ["
102+
+ str(interval.LowerLimit(firstPOI))
103+
+ ", "
104+
+ str(interval.UpperLimit(firstPOI))
105+
+ "]\n\n"
106+
)
107+
108+
# make a plot
109+
110+
print(
111+
"making a plot of the profile likelihood function ....(if it is taking a lot of time use less points or the "
112+
"TF1 drawing option)\n"
113+
)
114+
plot = ROOT.RooStats.LikelihoodIntervalPlot(interval)
115+
plot.SetNPoints(nScanPoints) # do not use too many points, it could become very slow for some models
116+
if poiXMin < poiXMax:
117+
plot.SetRange(poiXMin, poiXMax)
118+
opt = ROOT.TString("")
119+
if plotAsTF1:
120+
opt += ROOT.TString("tf1")
121+
plot.Draw(opt.Data()) # use option TF1 if too slow (plot.Draw("tf1")
122+
123+
# if requested perform also an hypothesis test for the significance
124+
if doHypoTest:
125+
nullparams = ROOT.RooArgSet("nullparams")
126+
nullparams.addClone(firstPOI)
127+
nullparams.setRealValue(firstPOI.GetName(), nullParamValue)
128+
pl.SetNullParameters(nullparams)
129+
print("Perform Test of Hypothesis : null Hypothesis is " + firstPOI.GetName() + str(nullParamValue))
130+
result = pl.GetHypoTest()
131+
print("\n>>>> Hypotheis Test Result ")
132+
result.Print()

‎tutorials/tmva/TMVA_Higgs_Classification.py‎

Lines changed: 12 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,7 @@
4444
useMLP = False # Multi Layer Perceptron (old TMVA NN implementation)
4545
useBDT = True # Boosted Decision Tree
4646
useDL = True # TMVA Deep learning ( CPU or GPU)
47-
useKeras = True # Use Keras Deep Learning via PyMVA
47+
useKeras = True # Use Keras Deep Learning via PyMVA
4848

4949
if ROOT.gSystem.GetFromPipe("root-config --has-tmva-pymva") == "yes":
5050
TMVA.PyMethodBase.PyInitialize()
@@ -73,13 +73,13 @@
7373

7474

7575
if ROOT.gSystem.AccessPathName(inputFileName):
76-
# file exists
7776
ROOT.Info("TMVA_Higgs_Classification", "Download Higgs_data.root file")
7877
TFile.SetCacheFileDir(".")
7978
inputFile = TFile.Open(inputFileLink, "CACHEREAD")
8079
if inputFile is None:
8180
raise FileNotFoundError("Input file cannot be downloaded - exit")
8281
else:
82+
# file exists
8383
inputFile = TFile.Open(inputFileName)
8484

8585

@@ -319,7 +319,7 @@
319319
Architecture=arch,
320320
)
321321

322-
#Keras DL
322+
# Keras DL
323323
if useKeras:
324324
ROOT.Info("TMVA_Higgs_Classification", "Building Deep Learning keras model")
325325
# create Keras model with 4 layers of 64 units and relu activations
@@ -329,13 +329,13 @@
329329
from tensorflow.keras.layers import Input, Dense
330330

331331
model = Sequential()
332-
model.add(Dense(64, activation='relu',input_dim=7))
333-
model.add(Dense(64, activation='relu'))
334-
model.add(Dense(64, activation='relu'))
335-
model.add(Dense(64, activation='relu'))
336-
model.add(Dense(2, activation='sigmoid'))
337-
model.compile(loss = 'binary_crossentropy', optimizer = Adam(learning_rate = 0.001), weighted_metrics = ['accuracy'])
338-
model.save('model_higgs.h5')
332+
model.add(Dense(64, activation="relu", input_dim=7))
333+
model.add(Dense(64, activation="relu"))
334+
model.add(Dense(64, activation="relu"))
335+
model.add(Dense(64, activation="relu"))
336+
model.add(Dense(2, activation="sigmoid"))
337+
model.compile(loss="binary_crossentropy", optimizer=Adam(learning_rate=0.001), weighted_metrics=["accuracy"])
338+
model.save("model_higgs.h5")
339339
model.summary()
340340

341341
if not os.path.exists("model_higgs.h5"):
@@ -353,7 +353,8 @@
353353
FilenameModel="model_higgs.h5",
354354
FilenameTrainedModel="trained_model_higgs.h5",
355355
NumEpochs=20,
356-
BatchSize=100 )
356+
BatchSize=100,
357+
)
357358
# GpuOptions="allow_growth=True",
358359
# ) # needed for RTX NVidia card and to avoid TF allocates all GPU memory
359360

0 commit comments

Comments
 (0)