13.10.3. Macro “calibrationLinBayesFlowrate1D.C

13.10.3.1. Objective

The goal here is to calibrate the parameter \(H_l\) within the flowrate model, while varying only two inputs (\(r_{\omega}\) and \(L\)). The remaining variables are fixed to the following values: \(r=25050\), \(T_u=89335\), \(T_l=89.55\), \(H_u=1050\), \(K_{\omega}=10950\). The context of this example has already been presented in Use-case for this chapter, including the model (implemented here as a C++ function) and the initial lines defining the TDataServer objects. This macro demonstrates how to apply a linear Bayesian estimation technique using a Relauncher-based architecture.

13.10.3.2. Macro Uranie

    // 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","tdsPar");
    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->setLikelihood("gauss_lin",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->drawResiduals("Residuals","*","","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("Parameters","*","nonewcanvas,transformed");

Much of this code has already been covered in the previous section Macro “calibrationMinimisationFlowrate1D.C” (up to the sequential run). The main difference here is that the input parameter is now defined as a TStochasticDistribution, representing the a priori chosen distribution. In this example, it can be either a TNormalDistribution or a TUniformDistribution (see [Bla17]), with the latter being selected here:

    tdsPar->addAttribute( new TUniformDistribution("hl", 700.0, 760.0) );

Another difference from the previous example is that the method must be slightly adapted to obtain values for the regressor. As discussed in [Bla17], linear Bayesian estimation can be applied when the model is approximately linear. This requires linearising the flowrate function, as demonstrated here:

\[f_{\theta}(x) = \left(2 \pi T_u\right)\left(\ln(\frac{r}{r_{\omega}})\left[ 1 + \frac{2 L T_u}{\ln(\frac{r}{r_{\omega}}) r_{\omega}^2 K_{\omega}} + \frac{T_u}{T_l}\right]\right)^{-1} \theta = H \times \theta\]

Here, the regressor can be expressed as

\[H = \left(2 \pi T_u\right)\left(\ln(\frac{r}{r_{\omega}}) \left[ 1 + \frac{2 L T_u}{\ln(\frac{r}{r_{\omega}}) r_{\omega}^2 K_{\omega}} + \frac{T_u}{T_l}\right]\right)^{-1}\]

From this expression, it is clear that we will be calibrating a newly defined parameter \(\theta = \left( H_u - H_l\right)\), which must later be converted back into the parameter of interest. To obtain the regressor estimation, we simply use another C++ function, flowrateModelnoH, together 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();
    }

This method also requires the input covariance matrix. The provided dataset (Ex2DoE_n100_sd1.75.dat) contains an estimate of the uncertainty affecting the reference measurements. Since this uncertainty is constant across all samples, the input covariance matrix is diagonal, with each diagonal element equal to the square of the standard deviation, as illustrated 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 defined (from a TCIntEval instance with the three input variables discussed above, in the correct order) along with the computation distribution method (sequential).

The calibration object is then created using a Mahalanobis distance function, which is used here for illustrative purposes (see Analytical linear Bayesian estimation for more details). The three key steps are:

  • Providing the input covariance matrix via the setObservationCovarianceMatrix method;

  • Specifying the regressors’ names using setRegressorName;

  • Defining the parameter transformation function (optional) with setParameterTransformationFunction.

The final step is more involved: since we are calibrating \(\theta\), we need to recover the corresponding parameter value at the end. This is achieved by providing a C++ function that transforms the parameter estimated from the linearisation back into the parameter of interest. For illustration, this is implemented in UserFunctions.C via the transf function, shown below:

void transf(double *x, double *res)
{
    res[0] = 1050 - x[0];  // simply H_l = \theta - H_u
}

The complete block of code discussed in this section is as follows:

    // Set the calibration object
    TLinearBayesian *cal = new TLinearBayesian(tdsPar,runner,1,"");
    cal->setLikelihood("gauss_lin",tdsRef,"rw:l","Qexp");
    cal->setObservationCovarianceMatrix(mat);
    cal->setRegressorName("H");
    cal->setParameterTransformationFunction(transf);
    cal->estimateParameters();

The final part demonstrates how to display the results. Since this method produces normal a posteriori distributions, they are represented by a vector of means and a covariance structure, both easily accessible. The means are displayed on screen, as illustrated in Console. Two additional pieces of a posteriori information are presented as plots: the residuals (Figure 13.63), which show the expected normal distribution behavior, as discussed in [Bla17] and the posterior distribution (Figure 13.64).

    // 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->drawResiduals("Residuals","*","","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("Parameters","*","nonewcanvas,transformed");

13.10.3.3. Console

--- Uranie v4.11/0 --- Developed with ROOT (6.36.06)
                      Copyright (C) 2013-2026 CEA/DES 
                      Contact: support-uranie@cea.fr 
                      Date: Thu Feb 12, 2026

***  TLinearBayesian::computeParameters(); Transformed parameters estimated are

1x1 matrix is as follows

     |      0    |
------------------
   0 |      749.7 

13.10.3.4. Graphs

../../_images/calibrationLinBayesFlowrate1D_Res.png

Figure 13.63 Residuals graph of the macro “calibrationLinBayesFlowrate1D.C”

../../_images/calibrationLinBayesFlowrate1D_Par.png

Figure 13.64 Parameter graph of the macro “calibrationLinBayesFlowrate1D.C”