Documentation
/ User's manual in C++
: 
| XIV.11. Macros Calibration | ||
|---|---|---|
 | Chapter XIV. Use-cases in C++ |  |
This section introduces few examples dealing with calibration in order to illustrate the different techniques introduced in Chapter XI and few of the options available as well.
The purpose here is to calibrate the value of 
that entered the flowrate model, when only two inputs have been varied
(
and
) while the rest of the
variables are set to a frozen value: 
, 
, 
, 
, 
. The context has been already discussed in Section XI.2.4
(including discussing the model, here using the model as a C++-function and the first few lines defining the TDataServer
objects). This macro shows how to use simple minimisation technique, with a Relauncher-architecture.
{
// Load the function flowrateCalib1D
gROOT->LoadMacro("UserFunctions.C");
// Input reference file
TString ExpData="Ex2DoE_n100_sd1.75.dat";
// define the reference
TDataServer *tdsRef = new TDataServer("tdsRef","doe_exp_Re_Pr");
tdsRef->fileDataRead(ExpData.Data());
// define the parameters
TDataServer *tdsPar = new TDataServer("tdsPar","pouet");
tdsPar->addAttribute( new TAttribute("hl", 700.0, 760.0) );
tdsPar->getAttribute("hl")->setDefaultValue(728.0);
// Create the output attribute
TAttribute *out = new TAttribute("out");
// Create interface to assessors
TCIntEval *Model = new TCIntEval("flowrateCalib1D");
Model->addInput(tdsPar->getAttribute("hl"));
Model->addInput(tdsRef->getAttribute("rw"));
Model->addInput(tdsRef->getAttribute("l"));
Model->addOutput(out);
// Set the runner
TSequentialRun *runner = new TSequentialRun(Model);
// Set the calibration object
TMinimisation *cal = new TMinimisation(tdsPar,runner,1);
cal->setDistanceAndReference("LS",tdsRef,"rw:l","Qexp");
TNloptSubplexe solv;
cal->setOptimProperties(&solv);
cal->estimateParameters();
// Draw the residuals
TCanvas *canRes = new TCanvas("CanRes","CanRes",1200,800);
TPad *apad = new TPad("apad","apad",0, 0.03, 1, 1); apad->Draw(); apad->cd();
cal->drawResidues("Residual title","*","","nonewcanvas");
}
Apart from the first lines discussed in Section XI.2.4, the important line is the
one defining the starting point of the minimisation. This can be done by calling the
setStartingPoint method of the TNlopt class, or simply by defining
default value for the parameter attributes. This is done here:
tdsPar->getAttribute("hl")->setDefaultValue(728.0);
This macro continues by defining the model and the way to run it. The instance created here, is a
TCIntEval which simply request the three input variables discussed above in the correct order. Here the first one has to be , the parameter that we want to calibrate, because of the way the
C++-function has been defined and then the two varying ones, ( and ) whose values are coming from the reference input file. Once done, the output
attribute is added (as our model computes only one variable) and the chosen distribution strategy is chosen to be
sequential.
// Create interface to assessors
TCIntEval *Model = new TCIntEval("flowrateCalib1D");
Model->addInput(tdsPar->getAttribute("hl"));
Model->addInput(tdsRef->getAttribute("rw"));
Model->addInput(tdsRef->getAttribute("l"));
Model->addOutput(out);
// Set the runner
TSequentialRun *runner = new TSequentialRun(Model);
Once done the calibration object (TMinimisation) is created and, as discussed in Section XI.2.2.1, the first object to be created is the distance function (here the
least-square one) through the setDistanceAndReference, that also defines the TDataServer that
contains reference data, the name of the reference inputs and the reference variable to which the output of the
model should be compared with. Finally the optimisation algorithm is defined by creating an instance of
TNloptSubplexe and then the parameters are estimated.
// Set the calibration object
TMinimisation *cal = new TMinimisation(tdsPar,runner,1);
cal->setDistanceAndReference("LS",tdsRef,"rw:l","Qexp");
TNloptSubplexe solv;
cal->setOptimProperties(&solv);
cal->estimateParameters();
The final part is how to represents part of the results. As this method is a point-estimation there is only one value so it is always displayed on screen, as shown in Section XIV.11.1.3. The other interesting point is to look at the residual, as discussed in [metho] and this is done in Figure XIV.96 which shows normally-distributed residual for the a posteriori estimations.
Processing calibrationMinimisationFlowrate1D.C...
--- Uranie v0.0/0 --- Developed with ROOT (6.32.02)
Copyright (C) 2013-2024 CEA/DES
Contact: support-uranie@cea.fr
Date: Tue Jan 09, 2024
|....:....|....:....|....:....|....:
.************************************************
* Row * tdsPar__n * hl.hl * agreement *
************************************************
* 0 * 0 * 749.72363 * 18.149368 *
************************************************
The purpose here is to calibrate the value of that entered the flowrate model, when only two inputs have been varied
( and
) while the rest of the
variables are set to a frozen value: , , , , . The context has been already discussed in Section XI.2.4
(including discussing the model, here using the model as a C++-function and the first few lines defining the TDataServer
objects). This macro shows how to use linear Bayesian estimation technique, with a Relauncher-architecture.
{
// Load the function flowrateCalib1D
gROOT->LoadMacro("UserFunctions.C");
// Input reference file
TString ExpData="Ex2DoE_n100_sd1.75.dat";
// define the reference
TDataServer *tdsRef = new TDataServer("tdsRef","doe_exp_Re_Pr");
tdsRef->fileDataRead(ExpData.Data());
// define the parameters
TDataServer *tdsPar = new TDataServer("tdsPar","pouet");
tdsPar->addAttribute( new TUniformDistribution("hl", 700.0, 760.0) );
// Create the output attribute
TAttribute *out = new TAttribute("out");
// Create interface to assessors
TCIntEval *Reg = new TCIntEval("flowrateModelnoH");
Reg->addInput(tdsRef->getAttribute("rw"));
Reg->addInput(tdsRef->getAttribute("l"));
Reg->addOutput(new TAttribute("H") );
TSequentialRun *runnoH = new TSequentialRun(Reg);
runnoH->startSlave();
if(runnoH->onMaster())
{
TLauncher2 l(tdsRef, runnoH);
l.solverLoop();
runnoH->stopSlave();
}
// Create interface to assessors
TCIntEval *Model = new TCIntEval("flowrateCalib1D");
Model->addInput(tdsPar->getAttribute("hl"));
Model->addInput(tdsRef->getAttribute("rw"));
Model->addInput(tdsRef->getAttribute("l"));
Model->addOutput(out);
// Set the runner
TSequentialRun *runner = new TSequentialRun(Model);
// Set the covariance matrix of the input reference
double sd=tdsRef->getValue("sd_eps",0);
TMatrixD mat(100,100);
for(unsigned int ival=0; ival<tdsRef->getNPatterns(); ival++)
mat(ival,ival)=(sd*sd);
// Set the calibration object
TLinearBayesian *cal = new TLinearBayesian(tdsPar,runner,1,"");
cal->setDistanceAndReference("Mahalanobis",tdsRef,"rw:l","Qexp");
cal->setObservationCovarianceMatrix(mat);
cal->setRegressorName("H");
cal->setParameterTransformationFunction(transf);
cal->estimateParameters();
// Draw the residuals
TCanvas *canRes = new TCanvas("CanRes","CanRes",1200,800);
TPad *padRes = new TPad("padRes","padRes",0, 0.03, 1, 1); padRes->Draw(); padRes->cd();
cal->drawResidues("Residual title","*","","nonewcanvas");
// Draw the parameters
TCanvas *canPar = new TCanvas("CanPar","CanPar",1200,800);
TPad *padPar = new TPad("padPar","padPar",0, 0.03, 1, 1); padPar->Draw(); padPar->cd();
cal->drawParameters("Parameter title","*","","nonewcanvas,transformed");
}
A very large fraction of this code has been discussed in Section XIV.11.1.2 (from the very start to the sequential run used). The main
difference is the fact that the input parameter is now defined as a TStochasticDistribution
- inheriting object which is the a priori chosen distribution. It can, in this case, only be a
either a TNormalDistribution or a TUniformDistribution (see [metho])
and the latter is the chosen solution here:
tdsPar->addAttribute( new TUniformDistribution("hl", 700.0, 760.0) );
Another difference with previous example is the fact that we need to prepare a bit the method as we need to get
values for the regressor. As seen in [metho], the linear Bayesian estimation is available when the model can be
considered linear. This means that one should linearise the flowrate function as done here by writing:
where the
regressor can be expressed as 
. From there, it is clear that we will be calibrating a newly defined parameter 
, so we will have to transform that back into our parameter of interest at some point. To get the
regressor estimation we simply use another C++-function flowrateModelnoH along with the
standard Relauncher approach:
// Create interface to assessors
TCIntEval *Reg = new TCIntEval("flowrateModelnoH");
Reg->addInput(tdsRef->getAttribute("rw"));
Reg->addInput(tdsRef->getAttribute("l"));
Reg->addOutput(new TAttribute("H") );
TSequentialRun *runnoH = new TSequentialRun(Reg);
runnoH->startSlave();
if(runnoH->onMaster())
{
TLauncher2 l(tdsRef, runnoH);
l.solverLoop();
runnoH->stopSlave();
}
Moving on, this method also needs the input covariance matrix. The input datasets provided
(Ex2DoE_n100_sd1.75.dat) does contain an estimation of the uncertainty affecting the reference
measurement, this uncertainty being constant throughout the sample, the input covariance is stated to be diagonal
with a constant value set to the standard deviation squared, as, done below:
// Set the covariance matrix of the input reference
double sd=tdsRef->getValue("sd_eps",0);
TMatrixD mat(100,100);
for(unsigned int ival=0; ival<tdsRef->getNPatterns(); ival++)
mat(ival,ival)=(sd*sd);
The model is anyway defined along with the way to distribute the computation, and then the calibration object is constructed with a Mahalanobis distance function (which is not so relevant as discussed in Section XI.4, as it is only used for illustration purpose). The three important steps are then providing
the input covariance matrix through the
setObservationCovarianceMatrixmethod;the regressors name, by calling
setRegressorName;the parameter transformation function (not compulsory) with the
setParameterTransformationFunction.
The last step is tricky: as we will be performing calibration to get 
, it would be nice to get the proper parameter value in the
end. This is possible if one provides a C++-function that transform the parameter estimated from the linearisation
back into our parameter of interest. This is done in UserFunctions.C for illustration purpose
as it contains this transf function shown below
void transf(double *x, double *res) { res[0] = 1050 - x[0]; } // simply H_l = \theta - H_u The full block of code discussed here is this one
// Set the calibration object
TLinearBayesian *cal = new TLinearBayesian(tdsPar,runner,1,"");
cal->setDistanceAndReference("Mahalanobis",tdsRef,"rw:l","Qexp");
cal->setObservationCovarianceMatrix(mat);
cal->setRegressorName("H");
cal->setParameterTransformationFunction(transf);
cal->estimateParameters();
The final part is how to represents part of the results. As this method gives normal a posteriori laws, there defined by a vector of means and a covariance structure which could simply be accessed and the former is shown in on screen, as shown in Section XIV.11.2.3. Two other a posteriori information can be seen as plots: the parameter distribution (shown in Figure XIV.98) and the residual, as discussed in [metho], shown in Figure XIV.97 which shows normally-distributed behaviour.
Processing calibrationLinBayesFlowrate1D.C...
--- Uranie v0.0/0 --- Developed with ROOT (6.32.02)
Copyright (C) 2013-2024 CEA/DES
Contact: support-uranie@cea.fr
Date: Tue Jan 09, 2024
*** TLinearBayesian::computeParameters(); Transformed parameters estimated are
1x1 matrix is as follows
| 0 |
------------------
0 | 749.7
The purpose here is to calibrate the value of that entered the flowrate model, when only two inputs have been varied
( and
) while the rest of the
variables are set to a frozen value: , , , , . The context has been already discussed in Section XI.2.4
(including discussing the model, here using the model as a C++-function and the first few lines defining the TDataServer
objects). This macro shows how to use the rejection ABC method, with a Relauncher-architecture.
{
// Load the function flowrateCalib1D
gROOT->LoadMacro("UserFunctions.C");
// Input reference file
TString ExpData="Ex2DoE_n100_sd1.75.dat";
// define the reference
TDataServer *tdsRef = new TDataServer("tdsRef","doe_exp_Re_Pr");
tdsRef->fileDataRead(ExpData.Data());
// define the parameters
TDataServer *tdsPar = new TDataServer("tdsPar","pouet");
tdsPar->addAttribute( new TUniformDistribution("hl", 700.0, 760.0) );
// Create the output attribute
TAttribute *out = new TAttribute("out");
// Create interface to assessors
TCIntEval *Model = new TCIntEval("flowrateCalib1D");
Model->addInput(tdsPar->getAttribute("hl"));
Model->addInput(tdsRef->getAttribute("rw"));
Model->addInput(tdsRef->getAttribute("l"));
Model->addOutput(out);
// Set the runner
TSequentialRun *runner = new TSequentialRun(Model);
// Set the calibration object
Int_t nABC = 100; Double_t eps = 0.05;
TRejectionABC *cal = new TRejectionABC(tdsPar, runner, nABC, "");
cal->setDistanceAndReference("LS",tdsRef,"rw:l","Qexp");
cal->setGaussianNoise("sd_eps");
cal->setPercentile(eps);
cal->estimateParameters();
// Compute statistics
tdsPar->computeStatistic();
cout << "The mean of hl is " << tdsPar->getAttribute("hl")->getMean() << endl;
cout << "The std of hl is " << tdsPar->getAttribute("hl")->getStd() << endl;
// Draw the parameters
TCanvas *canPar = new TCanvas("CanPar","CanPar",1200,800);
TPad *padPar = new TPad("padPar","padPar",0, 0.03, 1, 1); padPar->Draw(); padPar->cd();
cal->drawParameters("Parameter title","*","","nonewcanvas");
// Draw the residuals
TCanvas *canRes = new TCanvas("CanRes","CanRes",1200,800);
TPad *padRes = new TPad("padRes","padRes",0, 0.03, 1, 1); padRes->Draw(); padRes->cd();
cal->drawResidues("Residual title","*","","nonewcanvas");
}
A very large fraction of this code has been discussed in Section XIV.11.1.2 (from the very start to the sequential run used). The main
difference is the fact that the input parameter is now defined as a TStochasticDistribution
- inheriting object, as a sample will be generated to test locations:
tdsPar->addAttribute( new TUniformDistribution("hl", 700.0, 760.0) );
Apart from this, the model is defined along with the way to distribute the computation, and then the calibration
object is constructed by defining the number of elements in the final sample (nABC set to 100) and,
here, the percentile of events kept (eps set to 5 percent, which means of total number of estimation
of 2000 locations). As the code is deterministic, the uncertainty model is inserted through a random gaussian noise
whose standard deviation is defined event-by-event thanks to a variable in the observation dataserver. The distance
is also define and the estimation is performed.
// Set the calibration object
Int_t nABC = 100; Double_t eps = 0.05;
TRejectionABC *cal = new TRejectionABC(tdsPar, runner, nABC, "");
cal->setDistanceAndReference("LS",tdsRef,"rw:l","Qexp");
cal->setGaussianNoise("sd_eps");
cal->setPercentile(eps);
cal->estimateParameters();
The final part is how to represents part of the results. As this method gives a sample, the first two lines give basic statistical information, directly on screen, as shown in Section XIV.11.3.3. Two other a posteriori information can be seen as plots: the parameter distribution (shown in Figure XIV.100) and the residual, as discussed in [metho], shown in Figure XIV.99 which shows normally-distributed behaviour.
Processing calibrationRejectionABCFlowrate1D.C...
--- Uranie v0.0/0 --- Developed with ROOT (6.32.02)
Copyright (C) 2013-2024 CEA/DES
Contact: support-uranie@cea.fr
Date: Tue Jan 09, 2024
<URANIE::INFO>
<URANIE::INFO> *** URANIE INFORMATION ***
<URANIE::INFO> *** File[${SOURCEDIR}/meTIER/calibration/souRCE/TDistanceFunction.cxx] Line[601]
<URANIE::INFO> TDistanceFunction::setGaussianRandomNoise: gaussian random noise(s) is added using information in [sd_eps] to modify the output variable(s) [out].
<URANIE::INFO> *** END of URANIE INFORMATION ***
<URANIE::INFO>
<URANIE::INFO>
<URANIE::INFO> *** URANIE INFORMATION ***
<URANIE::INFO> *** File[${SOURCEDIR}/meTIER/calibration/souRCE/TRejectionABC.cxx] Line[107]
<URANIE::INFO> TRejectionABC::computeParameters:: 2000 evaluations will be performed !
<URANIE::INFO> *** END of URANIE INFORMATION ***
<URANIE::INFO>
A posteriori mean coordinates are : (749.926)
The mean of hl is 749.926
The std of hl is 1.97042
The purpose here is to calibrate the value of that entered the flowrate model, when only two inputs have been varied
( and
) while the rest of the
variables are set to a frozen value: , , , , . The context has been already discussed in Section XI.2.4
(including discussing the model, here using the model as a C++-function and the first few lines defining the TDataServer
objects). This macro shows how to use Markov-chain approach with the metropolis-hasting algorithm, using a
Relauncher-architecture.
{
// Load the function flowrateCalib1D
gROOT->LoadMacro("UserFunctions.C");
// Input reference file
TString ExpData="Ex2DoE_n100_sd1.75.dat";
// define the reference
TDataServer *tdsRef = new TDataServer("tdsRef","doe_exp_Re_Pr");
tdsRef->fileDataRead(ExpData.Data());
tdsRef->addAttribute("wei_exp","1./(sd_eps*sd_eps)");
// define the parameters
TDataServer *tdsPar = new TDataServer("tdsPar","pouet");
tdsPar->addAttribute( new TUniformDistribution("hl", 700.0, 760.0) );
// Create the output attribute
TAttribute *out = new TAttribute("out");
// Create interface to assessors
TCIntEval *Model = new TCIntEval("flowrateCalib1D");
Model->addInput(tdsPar->getAttribute("hl"));
Model->addInput(tdsRef->getAttribute("rw"));
Model->addInput(tdsRef->getAttribute("l"));
Model->addOutput(out);
// Set the runner
TSequentialRun *runner = new TSequentialRun(Model);
// Set the calibration object
TMetropHasting *cal = new TMetropHasting(tdsPar,runner,2000,"");
cal->setDistanceAndReference("weightedLS",tdsRef,"rw:l","Qexp","wei_exp");
cal->setNbDump(400);
cal->setAcceptationRatioRange(0.4, 0.45);
cal->estimateParameters();
// Quality assessment : Draw the trace the MCMC
TCanvas *canTr = new TCanvas("CanTr","CanTr",1200,800);
TPad *padTr = new TPad("padTr","padTr",0, 0.03, 1, 1); padTr->Draw(); padTr->cd();
cal->drawTrace("Trace title","*","","nonewcanvas");
// Draw the parameters
TCanvas *canPar = new TCanvas("CanPar","CanPar",1200,800);
TPad *padPar = new TPad("padPar","padPar",0, 0.03, 1, 1); padPar->Draw(); padPar->cd();
cal->drawParameters("Parameter title","*","","nonewcanvas");
// Draw the residuals
TCanvas *canRes = new TCanvas("CanRes","CanRes",1200,800);
TPad *padRes = new TPad("padRes","padRes",0, 0.03, 1, 1); padRes->Draw(); padRes->cd();
cal->drawResidues("Residual title","*","","nonewcanvas");
// Compute the auto-correlation
int burn=20; // Remove first 20 elements
vector<int> lag={1,5,10,20};
vector<double> autoCorr;
cal->getAutoCorrelation(lag, &autoCorr, burn);
cout<<"Autocorrelation are "<<autoCorr.size()<<":"<<endl;
for(unsigned il=0; il<lag.size(); il++)
{
cout<<"*** for lag="<<lag.at(il)<<": ";
for(unsigned ip=0; ip<cal->getNPar(); ip++)
cout<<autoCorr.at(ip*lag.size()+il)<<"; ";
cout<<endl;
}
}
A very large fraction of this code has been discussed in Section XIV.11.1.2 (from the very start to the sequential run used). The main
difference is the fact that the input parameter is now defined as a TStochasticDistribution
- inheriting object, as a sample will be generated to test locations:
tdsPar->addAttribute( new TUniformDistribution("hl", 700.0, 760.0) );
Apart from this, the model is defined along with the way to distribute the computation, and then the calibration object is constructed by defining the number of elements in the final sample (set to 2000). The distance function is then defined and two properties are set along: the threshold to which a new line is dump on screen to provide information and the acceptation ratio to be kept by playing on the standard deviation of the research (see Section XI.6). Finally the estimation is performed.
// Set the calibration object
TMetropHasting *cal = new TMetropHasting(tdsPar,runner,2000,"");
cal->setDistanceAndReference("weightedLS",tdsRef,"rw:l","Qexp","wei_exp");
cal->setNbDump(400);
cal->setAcceptationRatioRange(0.4, 0.45);
cal->estimateParameters();
The final part is how to represents part of the results. At first one should look at the trace to check for any peculiar trend and choose a burn-in threshold if needed (see Section XI.6), which is shown in Figure XIV.101. As this method gives a sample, the first two lines give basic statistical information, directly on screen, as shown in Section XIV.11.4.3. Two other a posteriori information can be seen as plots: the parameter distribution (shown in Figure XIV.103) and the residual, as discussed in [metho], shown in Figure XIV.102 which shows normally-distributed behaviour.
Finally, the auto-correlation of the resulting sample can be computed with different lag values (see Section XI.6), as by definition, elements from a Markov-chain are not fully independent. This is done
here by calling the getAutoCorrelation method which provides as many estimations as one
request, for the lag values. The results are show on screen (see Section XIV.11.4.3) and are used for post-processing analysis, as the trace plot
discussed above.
// Compute the auto-correlation
int burn=20;
vector<int> lag={1,5,10,20};
vector<double> autoCorr;
cal->getAutoCorrelation(lag, &autoCorr, burn);
cout<<"Autocorrelation are "<<autoCorr.size()<<":"<<endl;
for(unsigned il=0; il<lag.size(); il++)
{
cout<<"*** for lag="<<lag.at(il)<<": ";
for(unsigned ip=0; ip<cal->getNPar(); ip++)
cout<<autoCorr.at(ip*lag.size()+il)<<"; ";
cout<<endl;
}
Processing calibrationMetropHastingFlowrate1D.C...
--- Uranie v0.0/0 --- Developed with ROOT (6.32.02)
Copyright (C) 2013-2024 CEA/DES
Contact: support-uranie@cea.fr
Date: Tue Jan 09, 2024
400 events done
800 events done
1200 events done
1600 events done
A posteriori mean coordinates are : (749.66)
Autocorrelation are 4:
*** for lag=1: 0.409977;
*** for lag=5: 0.0365494;
*** for lag=10: 0.0138947;
*** for lag=20: 0.0773913;
The purpose here is to calibrate the value of 
and 
in a
very simple linear model where both parameters are respectively the constant term and the multiplicative
coefficient of the only input variable. The input file is linReg_Database.dat and the toy
model is stored in UserFunctions.C. This macro illustrates the steps in two dimensions to get
a sample, estimate the burn-in and lag (if needed) and plot the resulting distribution.
{
// Load the function flowrateCalib1D
gROOT->LoadMacro("UserFunctions.C");
// Input reference file
TString ExpData="linReg_Database.dat";
// Input reference file loaded
TDataServer *tdsRef = new TDataServer("tdsRef","tdsRef");
tdsRef->fileDataRead(ExpData.Data());
// Define the uncertainty model wih a guess
double sd_exp=0.2;
tdsRef->addAttribute("wei_exp",Form("1./(%g*%g)",sd_exp,sd_exp));
// Define the parameters
TDataServer *tdsPar = new TDataServer("tdsPar","poute");
double binf_search=-2.0, bsup_search=2.0;
tdsPar->addAttribute(new TUniformDistribution("t0", binf_search, bsup_search) );
tdsPar->addAttribute(new TUniformDistribution("t1", binf_search, bsup_search)) ;
// Create the output attribute
TAttribute *out = new TAttribute("out");
// Create interface to assessors
TCIntEval eval("Toy");
eval.addInput(tdsRef->getAttribute("x"));
eval.addInput(tdsPar->getAttribute("t0"));
eval.addInput(tdsPar->getAttribute("t1"));
eval.addOutput(out);
// Set the runner
TSequentialRun run(&eval);
// Set the calibration object
// Providing wild guess for value and variation range
vector<double> inval={0.8, -0.6}, std={0.4, 0.5};
int ns=12000;
TMetropHasting *cal = new TMetropHasting(tdsPar, &run, ns,"");
cal->setDistanceAndReference("weightedLS",tdsRef,"x","yExp","wei_exp");
cal->setNbDump(4000);
cal->setInitialisation(inval, std);
cal->estimateParameters();
// Quality assessment : Draw the trace the MCMC
TCanvas *canTr = new TCanvas("CanTr","CanTr",1200,800);
TPad *padTr = new TPad("padTr","padTr",0, 0.03, 1, 1); padTr->Draw(); padTr->cd();
cal->drawTrace("Trace title","*","","nonewcanvas");
// Quality assessment : Draw the trace the MCMC
TCanvas *canAcc = new TCanvas("CanAcc","CanAcc",1200,800);
TPad *padAcc = new TPad("padAcc","padAcc",0, 0.03, 1, 1); padAcc->Draw(); padAcc->cd();
cal->drawAcceptationRatio("AcceptRatio title","*","","nonewcanvas");
int burn=100;
// Compute the auto-correlation
vector<int> lag={1,3,6,10,20};
vector<double> autoCorr;
cal->getAutoCorrelation(lag, &autoCorr, burn);
cout<<"Autocorrelation are:"<<endl;
for(unsigned il=0; il<lag.size(); il++)
{
cout<<"*** for lag="<<lag.at(il)<<": ";
for(unsigned ip=0; ip<cal->getNPar(); ip++)
cout<<autoCorr.at(ip*lag.size()+il)<<"; ";
cout<<endl;
}
int mylag=6;
//Define a selection based on burn-in and lag
TString mycut=Form("(%s > %d) && ((%s %% %d) == 0)", tdsPar->getIteratorName(), burn, tdsPar->getIteratorName(), mylag);
// Draw the parameters
TCanvas *canPar = new TCanvas("CanPar","CanPar",1200,800);
TPad *padPar = new TPad("padPar","padPar",0, 0.03, 1, 1); padPar->Draw(); padPar->cd();
cal->drawParameters("Parameter title","*",mycut.Data(),"nonewcanvas");
// Draw the residuals
TCanvas *canRes = new TCanvas("CanRes","CanRes",1200,800);
TPad *padRes = new TPad("padRes","padRes",0, 0.03, 1, 1); padRes->Draw(); padRes->cd();
cal->drawResidues("Residual title","*","","nonewcanvas");
}This macro starts, as usual by defining both reference and parameter dataservers. The only specific lines here are these lines used later-on in which we define the uncertainty hypothesis, meaning a guess of the uncertainty by creating the weight variable (constant throughout the 30 reference bservations)
// Define the uncertainty model wih a guess
double sd_exp=0.2;
tdsRef->addAttribute("wei_exp",Form("1./(%g*%g)",sd_exp,sd_exp));
This macro continues by defining the model and the way to run it. The instance created here, is a
TCIntEval which simply request the three input variables discussed above in the correct order. Here the first one has to be the input variable, whose values are coming
from the reference datasets, while the other ones are the parameters to be calibrated, because of the way the
C++-function has been defined. Once done, the output attribute is added (as our model computes only one variable)
and the chosen distribution strategy is chosen to be sequential.
// Create interface to assessors
TCIntEval eval("Toy");
eval.addInput(tdsRef->getAttribute("x"));
eval.addInput(tdsPar->getAttribute("t0"));
eval.addInput(tdsPar->getAttribute("t1"));
eval.addOutput(out);
// Set the runner
TSequentialRun run(&eval);
Apart from this, the model is defined along with the way to distribute the computation, and then the calibration object is constructed by defining the number of elements in the final sample (set to 12000). The distance function is then defined and two properties are set along: the threshold to which a new line is dump on screen to provide information and the initialisation properties (values and variation ranges, see Section XI.6). Finally the estimation is performed.
// Set the calibration object
// Providing wild guess for value and variation range
vector<double> inval={0.8, -0.6}, std={0.4, 0.5};
int ns=12000;
TMetropHasting *cal = new TMetropHasting(tdsPar, &run, ns,"");
cal->setDistanceAndReference("weightedLS",tdsRef,"x","yExp","wei_exp");
cal->setNbDump(4000);
cal->setInitialisation(inval, std);
cal->estimateParameters();
The final part is how to represents part of the results. At first one should look at the trace to check for any peculiar trend and choose a burn-in threshold if needed (see Section XI.6), which is shown in Figure XIV.104, but one can also look at the acceptation ratio plots show in Figure XIV.105. As this method gives a sample, the first two lines give basic statistical information, directly on screen, as shown in Section XIV.11.5.3. One can also look at the autocorrelation and this might lead to the choice of a lag value to get low autocorrelation values (as shown in the console below).
Given both burn-in and lag values set, two other a posteriori information can be seen as plots: the parameter distribution (shown in Figure XIV.107) and the residual, as discussed in [metho], shown in Figure XIV.106 which shows normally-distributed behaviour.
Processing calibrationMetropHastingLinReg.C...
--- Uranie v0.0/0 --- Developed with ROOT (6.32.02)
Copyright (C) 2013-2024 CEA/DES
Contact: support-uranie@cea.fr
Date: Tue Jan 09, 2024
4000 events done
8000 events done
12000 events done
A posteriori mean coordinates are : (-0.441188,0.373825)
Autocorrelation are:
*** for lag=1: 0.598494; 0.629943;
*** for lag=3: 0.189151; 0.205056;
*** for lag=6: 0.0542879; 0.0667108;
*** for lag=10: 0.0160068; 0.0195018;
*** for lag=20: 0.0363133; 0.00163761;
The purpose here is to calibrate the value both of 
and that entered the flowrate model, when only two inputs have been varied
( and
) while the rest of the
variables are set to a frozen value: , , , . The context is the same as the one discussed in Section XIV.11.1 but it describes two things:
using Vizir instead of a more simple
TNloptSolver-inheriting instancediscuss the identifiability of a problem, introduced in [metho]
The model is the function flowrateCalib2D which is the same as the
flowrateClib1D just requesting the variable as first input.
{
// Load the function flowrateCalib2DVizir
gROOT->LoadMacro("UserFunctions.C");
// Input reference file
TString ExpData="Ex2DoE_n100_sd1.75.dat";
// define the reference
TDataServer *tdsRef = new TDataServer("tdsRef","doe_exp_Re_Pr");
tdsRef->fileDataRead(ExpData.Data());
// define the parameters
TDataServer *tdsPar = new TDataServer("tdsPar","pouet");
tdsPar->addAttribute( new TAttribute("hu", 1020.0, 1080.0) );
tdsPar->addAttribute( new TAttribute("hl", 720.0, 780.0) );
// Create the output attribute
TAttribute *out = new TAttribute("out");
// Create interface to assessors
TCIntEval *Model = new TCIntEval("flowrateCalib2D");
Model->addInput(tdsPar->getAttribute("hu"));
Model->addInput(tdsPar->getAttribute("hl"));
Model->addInput(tdsRef->getAttribute("rw"));
Model->addInput(tdsRef->getAttribute("l"));
Model->addOutput(out);
// Set the runner
TSequentialRun *runner = new TSequentialRun(Model);
// Set the calibration object
TMinimisation *cal = new TMinimisation(tdsPar,runner,1);
cal->setDistanceAndReference("relativeLS",tdsRef,"rw:l","Qexp");
// Set optimisaiton properties
URANIE::Reoptimizer::TVizirGenetic solv;
solv.setSize(24,15000,100);
cal->setOptimProperties(&solv);
// ((URANIE::Reoptimizer::TVizir2*)cal->getOptimMaster()->setTolerance(1e-6);
cal->estimateParameters();
// Draw the Residual
TCanvas *canRes = new TCanvas("CanRes","CanRes",1200,800);
TPad *padRes = new TPad("padRes","padRes",0, 0.03, 1, 1); padRes->Draw(); padRes->cd();
cal->drawResidues("tutu","*","","nonewcanvas");
// Draw the box plot of parameters
TCanvas *canPar = new TCanvas("CanPar","CanPar",1200,800);
tdsPar->getTuple()->SetMarkerStyle(20); tdsPar->getTuple()->SetMarkerSize(0.8);
tdsPar->Draw("hu:hl");
// Look at the correlation and statistic
tdsPar->computeStatistic("hu:hl");
TMatrixD corr=tdsPar->computeCorrelationMatrix("hu:hl");
corr.Print();
cout<<"hl is "<<tdsPar->getAttribute("hl")->getMean()<<" +- "<<tdsPar->getAttribute("hl")->getStd()<<endl;
cout<<"hu is "<<tdsPar->getAttribute("hu")->getMean()<<" +- "<<tdsPar->getAttribute("hu")->getStd()<<endl;
}
Apart from the first lines discussed in Section XI.2.4, the important line is the
one defining the variable, here as TAttribute with boundaries to define the phase space in
which the algorithm will look for:
tdsPar->addAttribute( new TAttribute("hu", 1020.0, 1080.0) );
tdsPar->addAttribute( new TAttribute("hl", 720.0, 780.0) );
This macro continues by defining the model and the way to run it. The instance created here, is a
TCIntEval which simply request the three input variables discussed above in the correct order. Here the first ones have to be and , the parameter that we want to calibrate, because of the way the C++-function has been defined
and then the two varying ones, ( and ) whose
values are coming from the reference input file. Once done, the output attribute is added (as our model computes
only one variable) and the chosen distribution strategy is chosen to be sequential.
// Create interface to assessors
TCIntEval *Model = new TCIntEval("flowrateCalib2D");
Model->addInput(tdsPar->getAttribute("hu"));
Model->addInput(tdsPar->getAttribute("hl"));
Model->addInput(tdsRef->getAttribute("rw"));
Model->addInput(tdsRef->getAttribute("l"));
Model->addOutput(out);
// Set the runner
TSequentialRun *runner = new TSequentialRun(Model);
Once done the calibration object (TMinimisation) is created and, as discussed in Section XI.2.2.1, the first object to be created is the distance function (here the
least-square one) through the setDistanceAndReference, that also defines the TDataServer that
contains reference data, the name of the reference inputs and the reference variable to which the output of the
model should be compared with. Finally the optimisation algorithm is defined by creating an instance of
TVizirGenetic and then the parameters are estimated.
// Set the calibration object
TMinimisation *cal = new TMinimisation(tdsPar,runner,1);
cal->setDistanceAndReference("relativeLS",tdsRef,"rw:l","Qexp");
// Set optimisaiton properties
URANIE::Reoptimizer::TVizirGenetic solv;
solv.setSize(24,15000,100);
cal->setOptimProperties(&solv);
// ((URANIE::Reoptimizer::TVizir2*)cal->getOptimMaster()->setTolerance(1e-6);
cal->estimateParameters();
The final part is how to represents part of the results. There are many interesting point in this discussion: the residual, which are estimated using the mean of both parameters, are shown in Figure XIV.108. The fact that it is normally-distributed residual for the a posteriori estimations shows that the model is correct even though looking at the second plots, the parameters distribution, shows that there is large variety of solutions possible, see Figure XIV.109. This is a problem of identifiability as there are an infinity of solutions that could give the same results, and this can be seen by looking at the correlation matrix shown in Section XIV.11.6.3.
Processing calibrationMinimisationFlowrate2DVizir.C...
--- Uranie v0.0/0 --- Developed with ROOT (6.32.02)
Copyright (C) 2013-2024 CEA/DES
Contact: support-uranie@cea.fr
Date: Tue Jan 09, 2024
first 100
Genetic 1
Generation : 1, rang max 23
Nb d'evaluation : 100, taille de la Z.P. : 0
Generation : 2, rang max 23
Nb d'evaluation : 465, taille de la Z.P. : 1
Generation : 3, rang max 8
Nb d'evaluation : 963, taille de la Z.P. : 6
Generation : 4, rang max 0
Nb d'evaluation : 1617, taille de la Z.P. : 24
Genetic converge 1617
************************************************************************************
* Row * tdsPar__n * hu.hu * hl.hl * agreement * rgpareto. * generatio *
************************************************************************************
* 0 * 0 * 1038.1302 * 738.13383 * 0.3639076 * 0 * 3 *
* 1 * 1 * 1080 * 780 * 0.3639083 * 0 * 3 *
* 2 * 2 * 1040.7058 * 740.77116 * 0.3639018 * 0 * 3 *
* 3 * 3 * 1079.9545 * 780 * 0.3639024 * 0 * 3 *
* 4 * 4 * 1040.7058 * 740.77116 * 0.3639018 * 0 * 1 *
* 5 * 5 * 1038.6638 * 738.70819 * 0.3639025 * 0 * 3 *
* 6 * 6 * 1040.7058 * 740.77116 * 0.3639018 * 0 * 3 *
* 7 * 7 * 1036.9776 * 736.98122 * 0.3639076 * 0 * 3 *
* 8 * 8 * 1036.9776 * 736.98122 * 0.3639076 * 0 * 3 *
* 9 * 9 * 1038.6638 * 738.70819 * 0.3639025 * 0 * 3 *
* 10 * 10 * 1036.9776 * 736.98122 * 0.3639076 * 0 * 2 *
* 11 * 11 * 1080 * 780 * 0.3639083 * 0 * 3 *
* 12 * 12 * 1036.9776 * 736.98122 * 0.3639076 * 0 * 3 *
* 13 * 13 * 1040.7058 * 740.77116 * 0.3639018 * 0 * 2 *
* 14 * 14 * 1080 * 780 * 0.3639083 * 0 * 3 *
* 15 * 15 * 1080 * 780 * 0.3639083 * 0 * 3 *
* 16 * 16 * 1038.6638 * 738.70819 * 0.3639025 * 0 * 3 *
* 17 * 17 * 1040.7058 * 740.77116 * 0.3639018 * 0 * 3 *
* 18 * 18 * 1080 * 780 * 0.3639083 * 0 * 3 *
* 19 * 19 * 1036.9776 * 736.98122 * 0.3639076 * 0 * 3 *
* 20 * 20 * 1036.9042 * 736.96108 * 0.3639019 * 0 * 3 *
* 21 * 21 * 1038.6638 * 738.70819 * 0.3639025 * 0 * 3 *
* 22 * 22 * 1080 * 780 * 0.3639083 * 0 * 3 *
* 23 * 23 * 1080 * 780 * 0.3639083 * 0 * 3 *
************************************************************************************
2x2 matrix is as follows
| 0 | 1 |
-------------------------------
0 | 1 1
1 | 1 1
hl is 752.445 +- 19.9463
hu is 1052.42 +- 19.9598
| |  | |
| XIV.10. Macros MetaModelOptim |  | XIV.12. Macros UncertModeler |