TKriging.h

Go to the documentation of this file.

// Copyright (C) 2013-2024 CEA/DES
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published
// by the Free Software Foundation, either version 3 of the License, or any
// later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
// $Id$
// $Author$
// $Date$
// $Revision$
// $State$

#ifndef TKRIGING_H
#define TKRIGING_H

#include <cmath>

#include "TObjArray.h"
#include "TFormula.h"
#include "TMath.h"

#include "TCorrelationFunction.h"
#include "TStandardEval.h"
#include "TDataServer.h"
#include "TAttribute.h"
#include "TDataServer.h"

#include <memory>

extern "C"
{
    void gpPred0(double *Kp, double variance, double *gamma, double *iL, 
                 int nnew, int n, double *yp, double *vp);
    void gpPred0Cov(double *Kp, double *Rp, double variance, double *gamma, double *iL,
                    int nnew, int n, double *yp, double *vp);
    void gpPredh(double *Hp, double *Kp, double variance, double *hbeta, double *gamma, double *iL, double *iR, double *M, 
                 int nnew, int n, int p, double *yp, double *vp);
    void gpPredhCov(double *Hp, double *Kp, double *Rp, double variance, double *hbeta, double *gamma, double *iL, double *iR, double *M, 
                    int nnew, int n, int p, double *yp, double *vp);
}
;

namespace URANIE
{
namespace Modeler
{
class TKriging: public URANIE::Relauncher::TDoubleEval
{

    // Attributes
private:
    Bool_t _blog; 
public:
    TString    _outputVarName;  
    TString    _inputVarNames;  
    TObjArray* _inputNamesList; 
    
    Int_t     _nX;          
    Int_t     _nS;          
    Int_t     _nYDeclared;  
    Int_t     _nXDeclared;  
    Double_t* _xObs;        
    Double_t* _xNormParams; 
    Double_t* _yObs;        
    Double_t* _iL;          
    Double_t* _gamma;       
    
    Double_t              _sigma2;              
    Double_t              _vErrMes;             
    TCorrelationFunction* _correlationFunction; 
    Double_t*             _corrLengths;         
    Double_t*             _normCorrLengths;     
    
    Bool_t     bUse_normalisation; 
    Bool_t     bHas_trend;     
    TString    _trendString;   
    vector<TFormula*> _trendCoefList; 
    vector<vector<int>> _trendIndex;
    std::unique_ptr<TObjArray> _inpname;
    Int_t      _nP;            
    Double_t*  _hbeta;         
    Double_t*  _vbeta;         
    Double_t*  _M;             
    Double_t*  _iR;            
    
    Double_t* _errorLoo; 
    Double_t* _varLoo;   
    Double_t  _mseLoo;   

    // Operations
public:
    //---------------------------------------------

     TKriging(Int_t nX, Int_t nS, Double_t* xObs, Double_t* xNormParams, Double_t* yObs, Double_t* gamma, Double_t* iL,
              Double_t gpVar, Double_t vErrMes, URANIE::Modeler::TCorrelationFunction* corrFunc, TString varNames,
              Int_t nP = 0, Double_t* hbeta = NULL, Double_t* vbeta= NULL, Double_t* M = NULL, Double_t* iR = NULL, TString trend = "");

#ifdef WITH_LIBXML2

     TKriging(TString filename, TString modelName);
#endif
    void initMatrices();

    virtual ~TKriging();

    //---------------------------------------------
    const char* getInputNames()
    {
        return _inputVarNames.Data();
    }

    const char* getOutputName()
    {
        return _outputVarName.Data();
    }
    
    Bool_t hasTrend()
    {
        return bHas_trend;
    }

    TCorrelationFunction* getCorrFunction()
    {
        return _correlationFunction;
    }

    Int_t getNInputs()
    {
        return _nX;
    }

    Int_t getNObs()
    {
        return _nS;
    }

    Int_t getNCorrFunctionParams()
    {
        return _correlationFunction->getNParameters();
    }

    Int_t getNCorrLengths()
    {
        return _correlationFunction->getNCorrLengths();
    }
    
    Int_t getNTrendParams()
    {
        return _nP;
    }

    Double_t* getObsInputs()
    {
        return _xObs;
    }

    Double_t* getObsOutputs()
    {
        return _yObs;
    }

    Double_t* getObsNormParams()
    {
        return _xNormParams;
    }
    
    Double_t* getGPParams()
    {
        return _gamma;
    }

    Double_t* getTrendParams()
    {
        return _hbeta;
    }

    Double_t* getTrendParamsCovMatrix()
    {
        return _vbeta;
    }
    
    vector<TFormula*> getTrendFormula()
    {
        return _trendCoefList;
    }

    TString getTrendString()
    {
        return _trendString;
    }
    
    Double_t getVariance()
    {
        return _sigma2;
    }


    Double_t getMeasurementErrorVariance()
    {
        return _vErrMes;
    }

    Double_t* getiL()
    {
        return _iL;
    }

    Double_t* getM()
    {
        return _M;
    }

    Double_t* getiR()
    {
        return _iR;
    }
    
    
    Double_t* getLooErrors()
    {
        return _errorLoo;
    }
    
    Double_t* getLooVariances()
    {
        return _varLoo;
    }
    
    void getLooErrors(double *arr, int size){
        try{
            if(size!=_nS)
                throw URANIE::Exceptions::UErrorExceptions(__FILE__, __LINE__,
                                                           Form("TKriging::getLooErrors method\n The size of provided array [%d] is not the size of the output to be dumped [%d]",
                                                                size, _nS));
        } catch (URANIE::Exceptions::UErrorExceptions& ue)
        {
            ue.printMessage();
            throw ue;
        }

        for(int ival=0; ival<size; ival++)
            arr[ival]=_errorLoo[ival];
        
    }
    
    void getLooVariances(double *arr, int size){
         try{
            if(size!=_nS)
                throw URANIE::Exceptions::UErrorExceptions(__FILE__, __LINE__,
                                                           Form("TKriging::getLooErrors method\n The size of provided array [%d] is not the size of the output to be dumped [%d]",
                                                                size, _nS));
        } catch (URANIE::Exceptions::UErrorExceptions& ue)
        {
            ue.printMessage();
            throw ue;
        }

        for(int ival=0; ival<size; ival++)
            arr[ival]=_varLoo[ival];
    }
    Double_t getLooRMSE()
    {
        return sqrt(_mseLoo);
    }

    Double_t getLooQ2()
    {
        Double_t mean = TMath::Mean(_nS, _yObs);
        Double_t sumVar = 0.0;
        for(Int_t i = 0; i < _nS; i++) sumVar += (_yObs[i]-mean)*(_yObs[i]-mean);
        if(_nS > 1)
        {
            Double_t varYObs = sumVar/(_nS-1);
            return 1.0 - _mseLoo/varYObs;
        }
        else
            return 1.0 - _mseLoo/sumVar;                
    }
    

    URANIE::DataServer::TDataServer* getLooData();
    

    Double_t* getCorrLengthsNormalised()
    {
        return _normCorrLengths;
    }
    

    Double_t* getCorrLengths()
    {
        return _corrLengths;
    }


    Double_t* getCorrFunctionParams()
    {
        return _correlationFunction->getParameters();
    }
    
    //---------------------------------------------

     void setCorrFunction(TCorrelationFunction* corrFunc)
     {
         _correlationFunction = corrFunc;
     }


    void setVariance(Double_t var)
    {
        _sigma2 = var;
    }


    void setLooErrors(Double_t *err, Double_t *var, Double_t mse);
    
    void createTrend();

    static void StatMat(std::string title, int n1, int n2, Double_t *v);
    static void StatMat(std::string title, int n1, int n2, Double_t *v,
                        std::vector<double> &tvmin, std::vector<double>  &tvmax,
                        std::vector<double> &tmean, std::vector<double>  &tvar,
                        std::vector<double> &tnorm);
    static void StatVect(std::string title, int n, Double_t *v);
    static void StatVect(std::string title, int n, Double_t *v, double &vmin, double &vmax, double &mean, double &var, double &norm);
    void computeCovarianceMatrix(int ns, double *xObs,
                                 bool bVariance_in_factor, 
                                 bool bHas_measurement_error, double alpha, double reg,
                                 double sigma2, double *Vmes,
                                 double *C, double *K);


    //---------------------------------------------
    

    Int_t eval(Double_t* x0, Double_t* y0, int=0);

#ifdef WITH_GPUANN
      void estimate(URANIE::DataServer::TDataServer *tdsQuery,
                    TString listOfInputs = "",
                    TString listOfOutputs = "",
                    Option_t * option = "",
                    Bool_t useGPU = true);
      void estimate_GPU(URANIE::DataServer::TDataServer *tdsQuery,
                        TString listOfInputs = "",
                        TString listOfOutputs = "",
                        Option_t * option = "");
#else
      void estimate(URANIE::DataServer::TDataServer *tdsQuery,
                    TString listOfInputs = "",
                    TString listOfOutputs = "",
                    Option_t * option = "",
                    Bool_t useGPU = false);
#endif
      void estimate_CPU(URANIE::DataServer::TDataServer *tdsQuery,
                        TString listOfInputs = "",
                        TString listOfOutputs = "",
                        Option_t * option = "");

    void estimateWithCov(URANIE::DataServer::TDataServer *tdsQuery,
                         TString listOfInputs = "",
                         TString listOfOutputs = "",
                         TString realValue = "",
                         Option_t * option = "");

    //---------------------------------------------

    void exportFunction(const char* lang, const char* file, const char* name, Option_t *option="");
    
    virtual void exportModelCplusplus(std::ofstream *sourcefile, const char* name = "");

    virtual void exportModelFortran(std::ofstream *sourcefile)
    {
        UNUSED(sourcefile);
        std::cout << "TKriging::exportModelFortran not yet implemented" << std::endl;
    }
    virtual void exportModelPMML(const char* file, const char* name, Option_t *option) const;

    void freeze();

    //---------------------------------------------
    virtual void printLog(Option_t *option = "");

    ClassDef(URANIE::Modeler::TKriging, ID_MODELER)
    //Classe de
};

} // Fin du namespace Modeler
} // Fin du namespace URANIE
#endif