13.6.12. Macro “sensitivitySobolLeveLE.py”

Warning

The levele command will be installed on your machine only if a Fortran compiler is found

13.6.12.1. Objective

The objective of this macro is to perform a full Sobol analysis on the temporal use-case levele. This use-case is an example of code that takes a dozen of entries in order to compute the evolution of dose as a function of time. The result of every computation consists in 3 vectors: the time (always the same value disregarding all entries), the dose called “y” and a third useless information.

13.6.12.2. Macro Uranie

"""
Example of Sobol estimation on the levelE code
"""
import sys
import numpy as np
from URANIE import Launcher, Sensitivity
from URANIE import DataServer as DS
import ROOT

# Define the DataServer
tds = DS.TDataServer("tdsLevelE", "levele")
tds.addAttribute(DS.TUniformDistribution("t", 100, 1000))
tds.addAttribute(DS.TLogUniformDistribution("kl", 0.001, 0.01))
tds.addAttribute(DS.TLogUniformDistribution("kc", 0.000001, 0.00001))
tds.addAttribute(DS.TLogUniformDistribution("v1", 0.001, 0.1))
tds.addAttribute(DS.TUniformDistribution("l1", 100, 500))
tds.addAttribute(DS.TUniformDistribution("r1", 1, 5))
tds.addAttribute(DS.TUniformDistribution("rc1", 3, 30))
tds.addAttribute(DS.TLogUniformDistribution("v2", 0.01, 0.1))
tds.addAttribute(DS.TUniformDistribution("l2", 50, 200))
tds.addAttribute(DS.TUniformDistribution("r2", 1, 5))
tds.addAttribute(DS.TUniformDistribution("rc2", 3, 30))
tds.addAttribute(DS.TLogUniformDistribution("w", 100000, 10000000))

# Tell the code where to find attribute value in input file
sIn = "levelE_input_with_values_rows.in"
tds.getAttribute("t").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)
tds.getAttribute("kl").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)
tds.getAttribute("kc").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)
tds.getAttribute("v1").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)
tds.getAttribute("l1").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)
tds.getAttribute("r1").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)
tds.getAttribute("rc1").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)
tds.getAttribute("v2").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)
tds.getAttribute("l2").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)
tds.getAttribute("r2").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)
tds.getAttribute("rc2").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)
tds.getAttribute("w").setFileKey(sIn, "", "%e", DS.TAttributeFileKey.kNewRow)

# How to read ouput files
out = Launcher.TOutputFileRow("_output_levelE_withRow_.dat")
# Tell the output file that attribute IS a vector and is SECOND column
out.addAttribute(DS.TAttribute("y", DS.TAttribute.kVector), 2)

# Creation of Launcher.TCode
myc = Launcher.TCode(tds, "levele 2> /dev/null")
myc.addOutputFile(out)

# Run Sobol analysis
tsobol = Sensitivity.TSobol(tds, myc, 10000)
tsobol.computeIndexes(args.computeIdxOpt) # args.computeIdxOpt = "" or "nointermed"

ntresu = tsobol.getResultTuple()

# Plotting mess
colors = [1, 2, 3, 4, 6, 7, 8, 15, 30, 38, 41, 46]
tps = np.array([20000, 30000, 40000, 50000, 60000, 70000, 80000, 90000, 100000,
                200000, 300000, 400000, 500000, 600000, 700000, 800000, 900000,
                1e+06, 2e+06, 3e+06, 4e+06, 5e+06, 6e+06, 7e+06, 8e+06, 9e+06],
               dtype=float)

c2 = ROOT.TCanvas("c2", "c2", 5, 64, 1200, 900)
pad = ROOT.TPad("pad", "pad", 0, 0.03, 1, 1)
pad.Draw()
pad.Divide(1, 2)
pad.cd(1)
ROOT.gPad.SetLogx()
ROOT.gPad.SetGrid()

# LegendandMArker
ROOT.gStyle.SetMarkerStyle(3)
ROOT.gStyle.SetLegendBorderSize(0)
ROOT.gStyle.SetFillStyle(0)

mg2 = ROOT.TMultiGraph()
names = ["t", "kl", "kc", "v1", "l1", "r1", "rc1",
         "v2", "l2", "r2", "rc2", "w"]
fdeg = [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]
leg = [0., 0., 0., 0., 0., 0.]
Zeros = np.zeros([26])

for igr in range(12):

    sel = "Inp==\""+names[igr]+"\""
    sel += " && Order==\"First\" &&  Algo==\"martinez11\""
    ntresu.Draw("CILower:CIUpper:Value", sel, "goff")
    data = ntresu.GetV3()
    themin = ntresu.GetV1()
    themax = ntresu.GetV2()
    for i in range(26):
        themin[i] = data[i] - themin[i]
        themax[i] = - data[i] + themax[i]

    if (igr % 2) == 0:
        leg[int(igr/2)] = ROOT.TLegend(0.1+0.15*(igr/2), 0.91,
                                       0.25+0.15*(igr/2), 0.98)
        leg[int(igr/2)].SetTextSize(0.045)

    fdeg[igr] = ROOT.TGraphAsymmErrors(26, tps, data,
                                       Zeros, Zeros, themin, themax)
    fdeg[igr].SetMarkerColor(colors[igr])
    fdeg[igr].SetLineColor(colors[igr])
    fdeg[igr].SetFillColor(colors[igr])
    leg[int(igr/2)].AddEntry(fdeg[igr], names[igr], "pl")

for igr2 in range(11, -1, -1):
    mg2.Add(fdeg[igr2])

mg2.Draw("APC")
mg2.GetXaxis().SetTitle("Time")
mg2.GetYaxis().SetTitle("S_{1}[martinez11]")
for igr3 in range(6):
    leg[igr3].Draw()

mg = ROOT.TMultiGraph()
tdeg = [0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]
pad.cd(2)
ROOT.gPad.SetLogx()
ROOT.gPad.SetGrid()
for igr in range(12):

    sel = "Inp==\""+names[igr]+"\""
    sel += " && Order==\"Total\" &&  Algo==\"martinez11\""
    ntresu.Draw("CILower:CIUpper:Value", sel, "goff")
    data = ntresu.GetV3()
    themin = ntresu.GetV1()
    themax = ntresu.GetV2()
    for i in range(26):
        themin[i] = data[i] - themin[i]
        themax[i] = - data[i] + themax[i]

    for ip in range(26):
        if ip == 0:
            tdeg[igr] = ROOT.TGraphAsymmErrors()
            tdeg[igr].SetMarkerColor(colors[igr])
            tdeg[igr].SetLineColor(colors[igr])
            tdeg[igr].SetFillColor(colors[igr])

        tdeg[igr].SetPoint(ip, tps[ip], data[ip])
        tdeg[igr].SetPointError(ip, 0, 0, themin[ip], themax[ip])

for igr2 in range(11, -1, -1):
    mg.Add(tdeg[igr2])

mg.Draw("APC")
mg.GetXaxis().SetTitle("Time")
mg.GetYaxis().SetTitle("S_{T}[martinez11]")
for igr3 in range(6):
    leg[igr3].Draw()

The levele external code is located in the bin directory of the Uranie installation.

When looking at the code and comparing it to an usual Sobol estimation, the organisation is completely transparent. The only noticeable (and compulsory) thing to do is to change the default type of the attribute read at the end of the job. This is done in this line:

out.addAttribute(DS.TAttribute("y", DS.TAttribute.kVector), 2)

where the output attribute is provided, changing its nature to a vector, thanks to the second argument of the TAttribute constructor from the default (kReal) to the desired nature (kVector). Once this is done, this information is broadcast internally to the code that knows how to deal with this type of attribute.

The rest of the code is the graphical part, leading to the figure below (it is provided to illustrate how to represent results).

13.6.12.3. Graph

The results of the previous macro is shown in Figure 13.28, where the evolution of the sobol coefficient is shown for all inputs with the uncertainty band, for the first order coefficient and the total one, respectively on the top and bottom panel.

Figure 13.28 Graph of the macro “sensitivitySobolLeveLE.py”