Documentation
/ User's manual in C++
: 
| VI.6. Fourier-based methods | ||
|---|---|---|
 | Chapter VI. The Sensitivity module |  |
This section is a short excerpt of [metho]
The Fourier Amplitude Sensitivity Test (FAST) [MCRAE198215, SALTELLI1998445] is a
procedure that provides a way to estimate the expected value and variance of the output variable of a model, along
with the contribution of the input factors to this variance. An advantage of it, is that the evaluation of
sensitivity can be carried out independently for each factor using just a set of runs because all the terms in a
Fourier expansion are mutually orthogonal. The main idea behind this procedure is to transform the 
-dimensional integration into a
single-dimension one, by using the transformation
where ideally, 
is a set of angular frequencies said to be incommensurate (meaning that no frequency can be
obtained by linear combination of the other ones when using integer coefficients) and 
is a transformation function chosen in order to ensure
that the variable is sampled accordingly with the probability density function of 
(meaning that they are all uniformly distributed in their
respective volume definition).
The first order coefficient is then obtained by estimating the variance for a fundamental 
and its harmonics. The important point to
notice, for a real computation, is the limitation of the sum that, in the previous equation, runs up to infinity. A
truncation is done by imposing a cut-off with a factor 
called the interference factor (whose default
value in Uranie is set to 6).
The Random Balance Design (RBD) [Tarantola2006717] method selects
design points over a curve
in the input space. The input space is explored here using the same frequency 
. However the curve is not space-filling, therefore, we
take random permutations of the coordinates of such points, to generate a set of scrambled points that cover the
input space. The model is then evaluated at each design point. Subsequently, the model outputs are re-ordered such
that the design points are in increasing order with respect to factor .
Thanks
to the use of permutations, the total cost is of the order of assessments instead of the order of 
for the FAST one.
Both TFAST and TRBD can be constructed either
from an external code or from a function. In the implementation done within Uranie there are several
modifiable parameters that can be considered before starting an analysis using the FAST method:
The transformation function
chosen among the following list:
Cukier:
SaltelliA:
SaltelliB:
In this list,
is the nominal value of the factor , 
denotes the endpoints that define the estimated range of uncertainty of , 
is a random phase shift taken value in 
and
evolves in 
.
The interference factor:
can be changed as well.
The frequencies: by providing a vector, it is possible to set a default at the frequencies' value used instead of having them determined by a specific algorithm to avoid, as best as possible, the interference.
The only common parameter changeable for both methods (and directly in the construction) is the number of
samples. Once the object is constructed, running the computeIndexes
will compute the sensitivity indices and the usual drawing methods will allow to retrieve and represent the results,
as already explained for other methods.
Warning
If the Fourier-based object is constructed to run a function, the output used for the estimation will be the first output provided by the function, unless the sixth argument of the constructor is filled.
In Uranie, computing sensitivity with the FAST method is dealt with the eponymous class TFast
which inherits from the TSensitivity class.
This class handles all the steps to compute the Sobol indices:
- generating the deterministic sample;
- running the code or the analytic function to get the response on the sample;
- computing the first order index for each input variable.
The example script below uses the TFast class to compute and display Sobol sensitivity indices:
{
gROOT->LoadMacro("UserFunctions.C");
// Define the DataServer
TDataServer *tds = new TDataServer("tdsflowreate", "DataBase flowreate");
tds->addAttribute( new TUniformDistribution("rw", 0.05, 0.15));
tds->addAttribute( new TUniformDistribution("r", 100.0, 50000.0));
tds->addAttribute( new TUniformDistribution("tu", 63070.0, 115600.0));
tds->addAttribute( new TUniformDistribution("tl", 63.1, 116.0));
tds->addAttribute( new TUniformDistribution("hu", 990.0, 1110.0));
tds->addAttribute( new TUniformDistribution("hl", 700.0, 820.0));
tds->addAttribute( new TUniformDistribution("l", 1120.0, 1680.0));
tds->addAttribute( new TUniformDistribution("kw", 9855.0, 12045.0));
// \param Size of a sampling.
Int_t nS = 500;
// Graph
TFast * tfast = new TFast(tds, "flowrateModel", nS);
tfast->computeIndexes("graph");
TCanvas *cc = new TCanvas("canhist", "histgramme",1);
tfast->drawIndexes("Flowrate", "", "nonewcanv,hist,first");
TCanvas *ccc = new TCanvas("canpie", "TPie",1);
TPad *apad = new TPad("apad","apad",0, 0.03, 1, 1); apad->Draw(); apad->cd();
tfast->drawIndexes("Flowrate", "", "nonewcanv,pie,first");
}The figures resulting from this script are shown below, the first shows the frequencies:
The two figures below display the sensitivity indices in a histogram and in a pie chart:
First, define the uncertain parameters and add them to a TDataServer object:
// Define the DataServer
TDataServer *tds = new TDataServer("tdsflowreate", "DataBase flowreate");
tds->addAttribute( new TUniformDistribution("rw", 0.05, 0.15));
tds->addAttribute( new TUniformDistribution("r", 100.0, 50000.0));
tds->addAttribute( new TUniformDistribution("tu", 63070.0, 115600.0));
tds->addAttribute( new TUniformDistribution("tl", 63.1, 116.0));
tds->addAttribute( new TUniformDistribution("hu", 990.0, 1110.0));
tds->addAttribute( new TUniformDistribution("hl", 700.0, 820.0));
tds->addAttribute( new TUniformDistribution("l", 1120.0, 1680.0));
tds->addAttribute( new TUniformDistribution("kw", 9855.0, 12045.0));
There are three different constructors to build a TFast
object, each corresponding to a different problem:
- the model is an analytic function run by Uranie,
- the model is a code run by Uranie.
- the model is either a function or a code and the problem is specified through a Relauncher architecture.
The constructor prototype used with an analytic function is:
// Create a TFast object with an analytic function
TFast(TDataServer *tds, void *fcn(double*,double*), const char *inp, const char *out, Int_t ns)
TFast(TDataServer *tds, const char *fcn, Int_t ns, const char *inp="", const char *out="", Option_t *option="")This constructor takes five arguments:
- a pointer to a
TDataServerobject, - a pointer to an analytic function (either a
voidor aconst charthat represents the function's name when it has been loaded in ROOT's memory), - an integer to specify the number of simulations, its default value is 100,
- a string to specify the names of the input factors separated by ':' (ex. "rw:r:tu:tl:hu:hl:l:kw"), its default value might be the empty string "",
- a string to specify the names of the output variables of the model, its default value might be the empty string "".
Here is an example of how to use the constructor with an analytic function:
Int_t ns = 50;
TFast * tfast = new TFast(tds, "flowrateModel", ns, "rw:r:tu:tl:hu:hl:l:kw", "flowrateModel"); The constructor prototype used with a code is:
// Create a TFast object with a code
TFast(TDataServer *tds, TCode *fcode, Int_t ns, Option_t *option="")This constructor takes three arguments:
- a pointer to a
TDataServerobject, - a pointer to a
TCode, - an integer to specify the number of simulations.
Here is an example of the use of this constructor on the flowrate case:
// The reference input file
TString sJDDReference = TString("flowrate_input_with_keys.in");
// Set the reference input file and the key for each input attributes
tds->getAttribute("rw")->setFileKey(sJDDReference, "Rw");
tds->getAttribute("r")->setFileKey(sJDDReference, "R");
tds->getAttribute("tu")->setFileKey(sJDDReference, "Tu");
tds->getAttribute("tl")->setFileKey(sJDDReference, "Tl");
tds->getAttribute("hu")->setFileKey(sJDDReference, "Hu");
tds->getAttribute("hl")->setFileKey(sJDDReference, "Hl");
tds->getAttribute("l")->setFileKey(sJDDReference, "L");
tds->getAttribute("kw")->setFileKey(sJDDReference, "Kw");
// The output file of the code
TOutputFileRow *fout = new TOutputFileRow("_output_flowrate_withRow_.dat");
// The attribute in the output file
fout->addAttribute(new TAttribute("yhat"));
// Instanciation de mon code
TCode *mycode = new TCode(tds, "flowrate -s -k");
// Adding the outputfile to the tcode object
mycode->addOutputFile( fout );
// Size of a sampling.
Int_t nS = 50;
TFast * tfastC = new TFast(tds, mycode, nS); The constructor prototype used with a runner is:
// Create a TFast object with a runner
TFast(TDataServer *tds, TRun *run, Int_t ns, Option_t *option="")This constructor takes three arguments:
- a pointer to a
TDataServerobject, - a pointer to a
TRun, - an integer to specify the number of simulations.
Here is an example of the use of this constructor on the flowrate code in a sequential mode:
// The input file
TKeyScript infile("flowrate_input_with_keys.in");
// provide the input and their key
infile.setInputs(8, tds->getAttribute("rw"), "Rw", tds->getAttribute("r"), "R",
tds->getAttribute("tu"), "Tu", tds->getAttribute("tl"), "Tl", tds->getAttribute("hu"), "Hu",
tds->getAttribute("hl"), "Hl", tds->getAttribute("l"), "L", tds->getAttribute("kw"), "Kw");
TAttribute yhat("yhat");
// The output file of the code
TKeyResult outfile("_output_flowrate_withKey_.dat");
// The attribute in the output file
outfile.addOutput(&yhat, "yhat");
// Instanciation de mon code
TCodeEval code("flowrate -s -k");
// Adding the intput/output file to the code
code.addInputFile(&infile);
code.addOutputFile(&outfile);
TSequentialRun run(&code);
run.startSlave();
if(run.onMaster())
{
Int_t nS = 50;
TFast * tfastR = new TFast(tds, &run, nS);
//...
}
To generate the FAST sample, use the generateSample method:
tfast->generateSample(); Then, the sample generated can be exported in a file and used outside of Uranie to compute the simulations associated.
To compute the indices, run the method computeIndexes:
tfast->computeIndexes(); Note that this method is all inclusive: it constructs the sample, launches the simulations and computes the indices. If the "graph" option is provided, a graph of frequency will be drawn, as the one shown in Figure VI.7
To display graphically the coefficients, use the method drawIndexes:
void TSensitivity::drawIndexes(TString sTitre, const char *select, Option_t * option)The method needs:
- a
TStringcontaining the title of the figure, - a string containing a selection (empty if no selection),
- a string containing the options of the graphics separated by commas.
Some of the options available are:
- "nonewcanv": to not create a new canvas,
- "pie": to display a pie chart,
- "hist": to display a histogram,
- "first": to display the indices of the first order (with "hist" only),
In our example the use of this method is:
TCanvas *cc = new TCanvas("canhist", "histgramme");
tfast->drawIndexes("Flowrate", "", "nonewcanv,hist,first");
cc->Print("appliUranieFlowrateFAST100Histogram.png");
TCanvas *ccc = new TCanvas("canpie", "TPie");
tfast->drawIndexes("Flowrate", "", "nonewcanv,pie");
ccc->Print("appliUranieFlowrateFAST100Pie.png");
The coefficients, once computed, are stored in a TTree. To get this
TTree, use the method TSensitivity::getResultTuple():
TTree * results = tfast->getResultTuple();
Several methods exist in ROOT to extract data from a TTree, it is advised to look for them
into the ROOT documentation. We propose two ways of extracting the value of each coefficient from the
TTree.
The first method uses 3 steps to extract an index:
-
scan the
TTreefor the chosen input variable (with a selection) in order to obtain its row number. In our example, if we chose the variable "hl", we'll use the command:
and in the resulting figure below we see that the first order index of "hl" is in the row 10:results->Scan("*","Inp==\"hl\"");************************************************************************* * Row * Out * Inp * Order * Method * Algo * Value * ************************************************************************* * 10 * flowrateM * hl * First * FAST * --first-- * 0.0413399 * * 11 * flowrateM * hl * Total * FAST * --total-- * 0.0413399 * *************************************************************************
-
set the entry of the
TTreeon this row with the methodGetEntry; -
get the value of the index with
GetValuemethod on the "Value" leaf of theTTree.
Below is an example of extraction of the index for "hl" in our flowrate case:
results->Scan("*","Inp==\"hl\"");
results->GetEntry(10);
Double_t S_Rw_Indexe = results->GetLeaf("Value")->GetValue(); The second method uses 2 steps to extract an index:
-
use the
Drawmethod with a selection to select the index, for example the selection for the first order index of "rw" is "Inp==\"rw\" && Algo==\"--first--\""; -
get the pointer on the value of the index with
GetV1method on theTTree.
Below is another example of extraction of the first order index for "rw" in our flowrate case:
results->Draw("Value","Inp==\"rw\" && Algo==\"--first--\"","goff");
Double_t S_Rw_IndexeS = results->GetV1()[0];
In Uranie, computing sensitivity with the method RBD is dealt with the eponymous class TRBD
which inherits from the TSensitivity class and the TFast.
This class handles all the steps to compute the Sobol indices:
- generating the deterministic sample;
- running the code or the analytic function to get the response on the sample;
- computing the first order index for each input variable.
The example script below uses the TRBD class to compute and display Sobol sensitivity
indices:
{
gROOT->LoadMacro("UserFunctions.C");
// Define the DataServer
TDataServer *tds = new TDataServer("tdsflowreate", "DataBase flowreate");
tds->addAttribute( new TUniformDistribution("rw", 0.05, 0.15));
tds->addAttribute( new TUniformDistribution("r", 100.0, 50000.0));
tds->addAttribute( new TUniformDistribution("tu", 63070.0, 115600.0));
tds->addAttribute( new TUniformDistribution("tl", 63.1, 116.0));
tds->addAttribute( new TUniformDistribution("hu", 990.0, 1110.0));
tds->addAttribute( new TUniformDistribution("hl", 700.0, 820.0));
tds->addAttribute( new TUniformDistribution("l", 1120.0, 1680.0));
tds->addAttribute( new TUniformDistribution("kw", 9855.0, 12045.0));
// Size of a sampling.
Int_t nS = 500;
// Graph
TRBD * trbd = new TRBD(tds, "flowrateModel", nS);
trbd->computeIndexes("graph");
TCanvas *cc = new TCanvas("canhist", "histogramme",1);
trbd->drawIndexes("Flowrate", "", "nonewcanv,hist,first");
TCanvas *ccc = new TCanvas("canpie", "TPie",1);
TPad *apad = new TPad("apad","apad",0, 0.03, 1, 1); apad->Draw(); apad->cd();
trbd->drawIndexes("Flowrate", "", "nonewcanv,pie");
}The figures resulting from this script are shown below, the first shows the frequencies:
The two figures below display the sensitivity indices in a histogram and in a pie chart:
First, define the uncertain parameters and add them to a TDataServer object:
// Define the DataServer
TDataServer *tds = new TDataServer("tdsflowreate", "DataBase flowreate");
tds->addAttribute( new TUniformDistribution("rw", 0.05, 0.15));
tds->addAttribute( new TUniformDistribution("r", 100.0, 50000.0));
tds->addAttribute( new TUniformDistribution("tu", 63070.0, 115600.0));
tds->addAttribute( new TUniformDistribution("tl", 63.1, 116.0));
tds->addAttribute( new TUniformDistribution("hu", 990.0, 1110.0));
tds->addAttribute( new TUniformDistribution("hl", 700.0, 820.0));
tds->addAttribute( new TUniformDistribution("l", 1120.0, 1680.0));
tds->addAttribute( new TUniformDistribution("kw", 9855.0, 12045.0));
There are three different constructors to build a TRBD
object, each corresponding to a different problem:
- the model is an analytic function run by Uranie,
- the model is a code run by Uranie.
- the model is either a function or a code and the problem is specified through a Relauncher architecture.
The constructor prototype used with an analytic function is:
// Create a TRBD object with an analytic function
TRBD(TDataServer *tds, void *fcn, const char *inp, const char *out, Int_t ns)
TRBD(TDataServer *tds, const char *fcn, Int_t ns, const char *inp="", const char *out="")This constructor takes five arguments:
- a pointer to a
TDataServerobject, - a pointer to an analytic function (either a
voidor aconst charthat represents the function's name when it has been loaded in ROOT's memory), - an integer to specify the number of simulations, its default value is 100,
- a string to specify the names of the input factors separated by ':' (ex. "rw:r:tu:tl:hu:hl:l:kw"), its default value might be the empty string "",
- a string to specify the names of the output variables of the model, its default value might be the empty string "".
Here is an example of how to use the constructor with an analytic function:
Int_t nS = 50;
TRBD * trbd = new TRBD(tds, "flowrateModel", nS, "rw:r:tu:tl:hu:hl:l:kw", "flowrateModel"); The constructor prototype used with a code is:
// Create a TRBD object with a code
TRBD(TDataServer *tds, TCode *fcode, Int_t ns, Option_t *option="")This constructor takes three arguments:
- a pointer to a
TDataServerobject, - a pointer to a
TCode, - an integer to specify the number of simulations.
Here is an example of the use of this constructor on the flowrate case:
// The reference input file
TString sJDDReference = TString("flowrate_input_with_keys.in");
// Set the reference input file and the key for each input attributes
tds->getAttribute("rw")->setFileKey(sJDDReference, "Rw");
tds->getAttribute("r")->setFileKey(sJDDReference, "R");
tds->getAttribute("tu")->setFileKey(sJDDReference, "Tu");
tds->getAttribute("tl")->setFileKey(sJDDReference, "Tl");
tds->getAttribute("hu")->setFileKey(sJDDReference, "Hu");
tds->getAttribute("hl")->setFileKey(sJDDReference, "Hl");
tds->getAttribute("l")->setFileKey(sJDDReference, "L");
tds->getAttribute("kw")->setFileKey(sJDDReference, "Kw");
// The output file of the code
TOutputFileRow *fout = new TOutputFileRow("_output_flowrate_withRow_.dat");
// The attribute in the output file
fout->addAttribute(new TAttribute("yhat"));
// Instanciation de mon code
TCode *mycode = new TCode(tds, "flowrate -s -k");
// Adding the outputfile to the tcode object
mycode->addOutputFile( fout );
// Size of a sampling.
Int_t nS2 = 500;
TRBD * trbd2 = new TRBD(tds, mycode, nS2); The constructor prototype used with a runner is:
// Create a TRBD object with a runner
TRBD(TDataServer *tds, TRun *frun, Int_t ns, Option_t *option="")This constructor takes three arguments:
- a pointer to a
TDataServerobject, - a pointer to a
TRun, - an integer to specify the number of simulations.
Here is an example of the use of this constructor on the flowrate code in sequential mode:
// The input file
TKeyScript infile("flowrate_input_with_keys.in");
// provide the input and their key
infile.setInputs(8, tds->getAttribute("rw"), "Rw", tds->getAttribute("r"), "R",
tds->getAttribute("tu"), "Tu", tds->getAttribute("tl"), "Tl", tds->getAttribute("hu"), "Hu",
tds->getAttribute("hl"), "Hl", tds->getAttribute("l"), "L", tds->getAttribute("kw"), "Kw");
TAttribute yhat("yhat");
// The output file of the code
TKeyResult outfile("_output_flowrate_withKey_.dat");
// The attribute in the output file
outfile.addOutput(&yhat, "yhat");
// Instanciation de mon code
TCodeEval code("flowrate -s -k");
// Adding the intput/output file to the code
code.addInputFile(&infile);
code.addOutputFile(&outfile);
TSequentialRun run(&code);
run.startSlave();
if(run.onMaster())
{
// Size of a sampling.
Int_t nS2 = 500;
TRBD * trbdR = new TRBD(tds, &run, nS2);
//...
}
To generate the RBD sample, use the generateSample method:
trbd->generateSample(); Then, the sample generated can be exported in a file and used outside of Uranie to compute the simulations associated.
To compute the indices, run the method computeIndexes:
trbd->computeIndexes();Note that this method is all inclusive: it constructs the sample, launches the simulations and computes the indices. If the "graph" option is provided, a graph of frequency will be drawn, as the one shown in Figure VI.10
To display graphically the coefficients, use the method drawIndexes:
void TSensitivity::drawIndexes(TString sTitre, const char *select, Option_t * option)The method needs:
- a
TStringcontaining the title of the figure, - a string containing a selection (empty if no selection),
- a string containing the options of the graphics separated by commas.
Some of the options available are:
- "nonewcanv": to not create a new canvas,
- "pie": to display a pie chart,
- "hist": to display a histogram,
- "first": to display the indices of the first order (with "hist" only),
- ...
In our example the use of this method is:
TCanvas *cc = new TCanvas("canhist", "histogramme");
trbd->drawIndexes("Flowrate", "", "nonewcanv,hist,first");
cc->Print("appliUranieFlowrate_RBD_50Histogram.png");
TCanvas *ccc = new TCanvas("canpie", "TPie",5,64,1270,560);
trbd->drawIndexes("Flowrate", "", "nonewcanv,pie");
ccc->Print("appliUranieFlowrate_RBD_50Pie.png");
The coefficients, once computed, are stored in a TTree. To get this
TTree, use the method TSensitivity::getResultTuple():
TTree * results = trbd->getResultTuple();
Several methods exist in ROOT to extract data from a TTree, it is advised to look for them
into the ROOT documentation. We propose two ways of extracting the value of each coefficient from the
TTree.
The first method uses 3 steps to extract an index:
-
scan the
TTreefor the chosen input variable (with a selection) in order to obtain its row number. In our example, if we chose the variable "hl", we'll use the command:
and in the resulting figure below we see that the first order index of "hl" is in the row 10:results->Scan("*","Inp==\"hl\"");************************************************************************* * Row * Out * Inp * Order * Method * Algo * Value * ************************************************************************* * 10 * flowrateM * hl * First * RBD * --first-- * 0.0392761 * * 11 * flowrateM * hl * Total * RBD * --total-- * 0.0392761 * *************************************************************************
-
set the entry of the
TTreeon this row with the methodGetEntry; -
get the value of the index with
GetValuemethod on the "Value" leaf of theTTree.
Below is an example of extraction of the index for "hl" in our flowrate case:
results->Scan("*","Inp==\"hl\"");
results->GetEntry(100);
Double_t S_Rw_Indexe = results->GetLeaf("Value")->GetValue(); The second method uses 2 steps to extract an index:
-
use the
Drawmethod with a selection to select the index, for example the selection for the first order index of "rw" is "Inp==\"rw\" && Algo==\"--first--\""; -
get the pointer on the value of the index with
GetV1method on theTTree.
Below is another example of extraction of the first order index for "rw" in our flowrate case:
results->Draw("Value","Inp==\"rw\" && Algo==\"--first--\"","goff");
Double_t S_Rw_IndexeS = results->GetV1()[0]; | |  | |
| VI.5. The Sobol method |  | VI.7. The Johnson relative weight |