Skip to content

Commit

Permalink
style: Update validation scripts to be more Pythonic (#1565)
Browse files Browse the repository at this point in the history
* Ensure `if __name__ == "__main__"` is used
* Apply isort and black to validation scripts
* Use f-strings
* Use list comprehensions over loops
* Use snake_case
  • Loading branch information
matthewfeickert authored Aug 27, 2021
1 parent 40bd916 commit ce2ffab
Show file tree
Hide file tree
Showing 6 changed files with 181 additions and 179 deletions.
52 changes: 26 additions & 26 deletions validation/makedata.py
Original file line number Diff line number Diff line change
@@ -1,37 +1,37 @@
import ROOT

import json
import sys

source_data = json.load(open(sys.argv[1]))
root_file = sys.argv[2]
import ROOT

binning = source_data['binning']
bindata = source_data['bindata']
if __name__ == "__main__":
source_data = json.load(open(sys.argv[1]))
root_file = sys.argv[2]

binning = source_data["binning"]
bin_data = source_data["bin_data"]

f = ROOT.TFile(root_file, 'RECREATE')
data = ROOT.TH1F('data', 'data', *binning)
for i, v in enumerate(bindata['data']):
data.SetBinContent(i + 1, v)
data.Sumw2()
out_file = ROOT.TFile(root_file, "RECREATE")
data = ROOT.TH1F("data", "data", *binning)
for idx, value in enumerate(bin_data["data"]):
data.SetBinContent(idx + 1, value)
data.Sumw2()

bkg = ROOT.TH1F('bkg', 'bkg', *binning)
for i, v in enumerate(bindata['bkg']):
bkg.SetBinContent(i + 1, v)
bkg.Sumw2()
bkg = ROOT.TH1F("bkg", "bkg", *binning)
for idx, value in enumerate(bin_data["bkg"]):
bkg.SetBinContent(idx + 1, value)
bkg.Sumw2()

if "bkgerr" in bin_data:
bkgerr = ROOT.TH1F("bkgerr", "bkgerr", *binning)

if 'bkgerr' in bindata:
bkgerr = ROOT.TH1F('bkgerr', 'bkgerr', *binning)
# shapesys must be as multiplicative factor
for idx, value in enumerate(bin_data["bkgerr"]):
bkgerr.SetBinContent(idx + 1, value / bkg.GetBinContent(idx + 1))
bkgerr.Sumw2()

# shapesys must be as multiplicative factor
for i, v in enumerate(bindata['bkgerr']):
bkgerr.SetBinContent(i + 1, v / bkg.GetBinContent(i + 1))
bkgerr.Sumw2()
sig = ROOT.TH1F("sig", "sig", *binning)
for idx, value in enumerate(bin_data["sig"]):
sig.SetBinContent(idx + 1, value)
sig.Sumw2()

sig = ROOT.TH1F('sig', 'sig', *binning)
for i, v in enumerate(bindata['sig']):
sig.SetBinContent(i + 1, v)
sig.Sumw2()
f.Write()
out_file.Write()
63 changes: 31 additions & 32 deletions validation/run_cls.py
Original file line number Diff line number Diff line change
@@ -1,45 +1,44 @@
import ROOT
import sys

infile = sys.argv[1]

infile = ROOT.TFile.Open(infile)
workspace = infile.Get("combined")
data = workspace.data("obsData")

import ROOT

sbModel = workspace.obj("ModelConfig")
poi = sbModel.GetParametersOfInterest().first()
if __name__ == "__main__":
infile = sys.argv[1]

sbModel.SetSnapshot(ROOT.RooArgSet(poi))
infile = ROOT.TFile.Open(infile)
workspace = infile.Get("combined")
data = workspace.data("obsData")

bModel = sbModel.Clone()
bModel.SetName("bonly")
poi.setVal(0)
bModel.SetSnapshot(ROOT.RooArgSet(poi))
sb_model = workspace.obj("ModelConfig")
poi = sb_model.GetParametersOfInterest().first()

ac = ROOT.RooStats.AsymptoticCalculator(data, bModel, sbModel)
ac.SetPrintLevel(10)
ac.SetOneSided(True)
ac.SetQTilde(True)
sb_model.SetSnapshot(ROOT.RooArgSet(poi))

calc = ROOT.RooStats.HypoTestInverter(ac)
calc.RunFixedScan(51, 0, 5)
calc.SetConfidenceLevel(0.95)
calc.UseCLs(True)
bkg_model = sb_model.Clone()
bkg_model.SetName("bonly")
poi.setVal(0)
bkg_model.SetSnapshot(ROOT.RooArgSet(poi))

calc = ROOT.RooStats.AsymptoticCalculator(data, bkg_model, sb_model)
calc.SetPrintLevel(10)
calc.SetOneSided(True)
calc.SetQTilde(True)

result = calc.GetInterval()
test_inverter = ROOT.RooStats.HypoTestInverter(calc)
test_inverter.RunFixedScan(51, 0, 5)
test_inverter.SetConfidenceLevel(0.95)
test_inverter.UseCLs(True)

plot = ROOT.RooStats.HypoTestInverterPlot("plot", "plot", result)
c = ROOT.TCanvas()
c.SetLogy(False)
plot.Draw("OBS EXP CLb 2CL")
c.Draw()
c.SaveAs('scan.pdf')
result = test_inverter.GetInterval()

plot = ROOT.RooStats.HypoTestInverterPlot("plot", "plot", result)
canvas = ROOT.TCanvas()
canvas.SetLogy(False)
plot.Draw("OBS EXP CLb 2CL")
canvas.Draw()
canvas.SaveAs("scan.pdf")

print(f'observed: {result.UpperLimit()}')
print(f"observed: {result.UpperLimit()}")

for i in [-2, -1, 0, 1, 2]:
print(f'expected {i}: {result.GetExpectedUpperLimit(i)}')
for n_sigma in [-2, -1, 0, 1, 2]:
print(f"expected {n_sigma}: {result.GetExpectedUpperLimit(n_sigma)}")
61 changes: 30 additions & 31 deletions validation/run_single.py
Original file line number Diff line number Diff line change
@@ -1,43 +1,42 @@
import ROOT
import json
import sys

infile = sys.argv[1]

infile = ROOT.TFile.Open(infile)
workspace = infile.Get("combined")
data = workspace.data("obsData")

sbModel = workspace.obj("ModelConfig")
poi = sbModel.GetParametersOfInterest().first()

sbModel.SetSnapshot(ROOT.RooArgSet(poi))
import ROOT

bModel = sbModel.Clone()
bModel.SetName("bonly")
poi.setVal(0)
bModel.SetSnapshot(ROOT.RooArgSet(poi))
if __name__ == "__main__":
infile = sys.argv[1]

ac = ROOT.RooStats.AsymptoticCalculator(data, bModel, sbModel)
ac.SetPrintLevel(10)
ac.SetOneSided(True)
ac.SetQTilde(True)
infile = ROOT.TFile.Open(infile)
workspace = infile.Get("combined")
data = workspace.data("obsData")

sb_model = workspace.obj("ModelConfig")
poi = sb_model.GetParametersOfInterest().first()

calc = ROOT.RooStats.HypoTestInverter(ac)
calc.SetConfidenceLevel(0.95)
calc.UseCLs(True)
calc.RunFixedScan(1, 1, 1)
sb_model.SetSnapshot(ROOT.RooArgSet(poi))

result = calc.GetInterval()
bkg_model = sb_model.Clone()
bkg_model.SetName("bonly")
poi.setVal(0)
bkg_model.SetSnapshot(ROOT.RooArgSet(poi))

index = 0
calc = ROOT.RooStats.AsymptoticCalculator(data, bkg_model, sb_model)
calc.SetPrintLevel(10)
calc.SetOneSided(True)
calc.SetQTilde(True)

w = result.GetExpectedPValueDist(index)
v = w.GetSamplingDistribution()
test_inverter = ROOT.RooStats.HypoTestInverter(calc)
test_inverter.SetConfidenceLevel(0.95)
test_inverter.UseCLs(True)
test_inverter.RunFixedScan(1, 1, 1)

CLs_obs = result.CLs(index)
CLs_exp = list(v)[3:-3]
result = test_inverter.GetInterval()

import json
index = 0
dist = result.GetExpectedPValueDist(index)
CLs_obs = result.CLs(index)
CLs_exp = list(dist.GetSamplingDistribution())[3:-3]

print(json.dumps({'CLs_obs': CLs_obs, 'CLs_exp': CLs_exp}, sort_keys=True, indent=4))
print(
json.dumps({"CLs_obs": CLs_obs, "CLs_exp": CLs_exp}, sort_keys=True, indent=4)
)
40 changes: 21 additions & 19 deletions validation/run_single_q0.py
Original file line number Diff line number Diff line change
@@ -1,36 +1,38 @@
if __name__ == "__main__":
import ROOT
import sys
import json
import json
import sys

import ROOT

if __name__ == "__main__":
infile = sys.argv[1]

infile = ROOT.TFile.Open(infile)
workspace = infile.Get("combined")
data = workspace.data("obsData")

sbModel = workspace.obj("ModelConfig")
poi = sbModel.GetParametersOfInterest().first()
sb_model = workspace.obj("ModelConfig")
poi = sb_model.GetParametersOfInterest().first()
poi.setVal(1)
sbModel.SetSnapshot(ROOT.RooArgSet(poi))
sb_model.SetSnapshot(ROOT.RooArgSet(poi))

bModel = sbModel.Clone()
bModel.SetName("bonly")
bkg_model = sb_model.Clone()
bkg_model.SetName("bonly")
poi.setVal(0)
bModel.SetSnapshot(ROOT.RooArgSet(poi))
bkg_model.SetSnapshot(ROOT.RooArgSet(poi))

ac = ROOT.RooStats.AsymptoticCalculator(data, sbModel, bModel)
ac.SetPrintLevel(10)
ac.SetOneSidedDiscovery(True)
calc = ROOT.RooStats.AsymptoticCalculator(data, sb_model, bkg_model)
calc.SetPrintLevel(10)
calc.SetOneSidedDiscovery(True)

result = ac.GetHypoTest()
result = calc.GetHypoTest()
pnull_obs = result.NullPValue()
palt_obs = result.AlternatePValue()
pnull_exp = []
for sigma in [-2, -1, 0, 1, 2]:
usecls = 0
pnull_exp.append(ac.GetExpectedPValues(pnull_obs, palt_obs, sigma, usecls))
usecls = 0
pnull_exp = [
calc.GetExpectedPValues(pnull_obs, palt_obs, sigma, usecls)
for sigma in [-2, -1, 0, 1, 2]
]

print(
json.dumps({'p0_obs': pnull_obs, 'p0_exp': pnull_exp}, sort_keys=True, indent=4)
json.dumps({"p0_obs": pnull_obs, "p0_exp": pnull_exp}, sort_keys=True, indent=4)
)
61 changes: 31 additions & 30 deletions validation/run_toys.py
Original file line number Diff line number Diff line change
@@ -1,55 +1,57 @@
import json
import numpy as np
import sys

import matplotlib.pyplot as plt
import numpy as np
import ROOT

import pyhf
import pyhf.contrib.viz.brazil as brazil
import ROOT
import sys


def run_toys_ROOT(infile, ntoys):
infile = ROOT.TFile.Open(infile)
workspace = infile.Get("combined")
data = workspace.data("obsData")

sbModel = workspace.obj("ModelConfig")
poi = sbModel.GetParametersOfInterest().first()
sb_model = workspace.obj("ModelConfig")
poi = sb_model.GetParametersOfInterest().first()

sbModel.SetSnapshot(ROOT.RooArgSet(poi))
sb_model.SetSnapshot(ROOT.RooArgSet(poi))

bModel = sbModel.Clone()
bModel.SetName("bonly")
bkg_model = sb_model.Clone()
bkg_model.SetName("bonly")
poi.setVal(0)
bModel.SetSnapshot(ROOT.RooArgSet(poi))
bkg_model.SetSnapshot(ROOT.RooArgSet(poi))

ac = ROOT.RooStats.FrequentistCalculator(data, bModel, sbModel)
ac.SetToys(ntoys, ntoys)
calc = ROOT.RooStats.FrequentistCalculator(data, bkg_model, sb_model)
calc.SetToys(ntoys, ntoys)

profll = ROOT.RooStats.ProfileLikelihoodTestStat(bModel.GetPdf())
profll.SetOneSidedDiscovery(False)
profll.SetOneSided(True)
ac.GetTestStatSampler().SetTestStatistic(profll)
profile_ll = ROOT.RooStats.ProfileLikelihoodTestStat(bkg_model.GetPdf())
profile_ll.SetOneSidedDiscovery(False)
profile_ll.SetOneSided(True)
calc.GetTestStatSampler().SetTestStatistic(profile_ll)

calc = ROOT.RooStats.HypoTestInverter(ac)
calc.SetConfidenceLevel(0.95)
calc.UseCLs(True)
test_inverter = ROOT.RooStats.HypoTestInverter(calc)
test_inverter.SetConfidenceLevel(0.95)
test_inverter.UseCLs(True)

npoints = 2
n_points = 2

calc.RunFixedScan(npoints, 1.0, 1.2)
test_inverter.RunFixedScan(n_points, 1.0, 1.2)

result = calc.GetInterval()
result = test_inverter.GetInterval()

plot = ROOT.RooStats.HypoTestInverterPlot("plot", "plot", result)

data = []
for i in range(npoints):
d = {
"test_b": list(result.GetAltTestStatDist(i).GetSamplingDistribution()),
"test_s": list(result.GetNullTestStatDist(i).GetSamplingDistribution()),
"pvals": list(result.GetExpectedPValueDist(i).GetSamplingDistribution()),
data = [
{
"test_b": list(result.GetAltTestStatDist(idx).GetSamplingDistribution()),
"test_s": list(result.GetNullTestStatDist(idx).GetSamplingDistribution()),
"pvals": list(result.GetExpectedPValueDist(idx).GetSamplingDistribution()),
}
data.append(d)
for idx in range(n_points)
]

json.dump(data, open("scan.json", "w"))

Expand All @@ -58,6 +60,7 @@ def run_toys_ROOT(infile, ntoys):
plot.Draw("OBS EXP CLb 2CL")
canvas.GetListOfPrimitives().At(0).GetYaxis().SetRangeUser(0, 0.2)
canvas.Draw()

extensions = ["pdf", "png"]
for ext in extensions:
canvas.SaveAs(f"poi_scan_ROOT.{ext}")
Expand All @@ -83,8 +86,6 @@ def run_toys_pyhf(ntoys=2_000, seed=0):

fig, ax = plt.subplots()
brazil.plot_results(test_mus, fit_results, ax=ax)
ax.set_xlabel(r"$\mu$")
ax.set_ylabel(r"$\mathrm{CL}_{s}$")
_buffer = 0.02
ax.set_xlim(1.0 - _buffer, 1.2 + _buffer)
ax.set_ylim(0.0, 0.2)
Expand Down
Loading

0 comments on commit ce2ffab

Please sign in to comment.