TSumOfSquares.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 TSUMOFSQUARES_H
#define TSUMOFSQUARES_H
 
// ROOT
#include "TH1.h"
#include "TF1.h"
#include "TPad.h"
#include "TRandom3.h"
#include "TVirtualFitter.h"
#include "TStyle.h"
#include "Minuit2/FCNBase.h"
//fga 2015/01/12 #include "TFitterMinuit.h"
#include "TSystem.h"
#include "TMatrixT.h"
 
#include <vector>
#include <iostream>
 
// URANIE
#include "Optimizer.h"
#include "TCode.h"
#include "TDataServer.h"
#include "TObjective.h"
 
namespace URANIE
{
namespace Optimizer
{
class TSumOfSquares: public ROOT::Minuit2::FCNBase
{
private:
    vector<TObjective *> _objectives;              
 
public:
 
    TSumOfSquares() :
             _objectives(0), _ninput(0), _noutput(0), _tds(NULL), _code(NULL), _sCriteria("")
    {
    }
    TSumOfSquares(URANIE::DataServer::TDataServer *tds,
            URANIE::Launcher::TCode *code) :
            _objectives(0), _tds(tds), _code(code), _sCriteria("") 
    {
        init();
    }
 
    ~TSumOfSquares()
    {
        _dinputValue = NULL;
        delete[] _dinputValue;
        _doutputValue = NULL;
        delete[] _doutputValue;
 
    }
 
    void init();
 
    double operator()(const vector<double> & par) const
    {
        Int_t niter = _tds->getTuple()->GetEntries();
        Double_t fval = -0.123456789;
        _tds->getAttribute(_sCriteria)->getDefaultValue(fval);
        //      cout << " ** fval [ " << fval << "]" << endl;
 
        Double_t dObjectives = 0.;
 
        if (_code != NULL)
        {
            if (_code->getLog())
            {
                cout << endl << " ****** Step[" << niter << "]" << endl;
            }
            _dinputValue[0] = 1.0 + niter;
            // Remark: If we want to save all computations, we need a constant Bool_t _bSave,
            // and call the the method _code->init() at each new computation.
            _code->init();
            // Pretraitment of the data
            Double_t valcrt = 1.23456789;
            Int_t ncpt = 0;
 
            // fga::modif le 2008-10-01 07:49:37
            // attributes kConstant must be forgotten
            for (int i = 0; i < _ninput; i++)
            {
                // end of fga::modif le 2008-10-01 07:49:37
                URANIE::DataServer::TAttribute *att = _tds->getAttribute(i);
                if (_code->getLog())
                    cout << " ** " << i << "/" << _ninput << " name ["
                            << att->GetName() << "]";
                URANIE::DataServer::TAttribute::EOrigin norigin =
                        att->getOrigin();
                switch (norigin)
                {
                case URANIE::DataServer::TAttribute::kConstant:
                    if (_code->getLog())
                        cout << " Origin[kConstant]";
                    att->getDefaultValue(valcrt);
                    break;
                case URANIE::DataServer::TAttribute::kAttribute:
                    if (_code->getLog())
                        cout << " Origin[kAttribute]";
                    valcrt = par[ncpt++];
                    break;
                default:
                    break;
                }
                _dinputValue[1 + i] = valcrt;
                if (_code->getLog())
                    cout << " val[" << _dinputValue[1 + i] << "]" << endl;
            }
            int ind = 1;
            _code->preTraitment(ind, _dinputValue);
            // Launch the code
            _code->run();
            // Post-traitment of the data
            ind = 0;
            _code->postTraitment(ind, _doutputValue, "", _listOfDataOutput);
            if (_code->getLog())
            {
                cout << endl << "*************************************" << endl
                        << " *** TSumOfSquares operator () _noutput["
                        << _noutput << "] ind[" << ind << "]" << endl;
                for (Int_t ii = 0; ii <= _noutput; ii++)
                    cout << " ** ii(" << ii << "/" << _noutput << ") val("
                            << _doutputValue[ii] << ")" << endl;
                cout << " *** End Of TFCNCode operator ()" << endl
                        << "*************************************" << endl;
            }
 
            TMatrixT<double> *lastmat =
                    (TMatrixT<double> *) _listOfDataOutput->Last();
            if (_code->getLog())
                lastmat->Print(Form("Matrice de sortie iteration[%d]", niter));
            Int_t nrows = lastmat->GetNrows();
            Int_t ncols = lastmat->GetNcols();
            if (_code->getLog())
                cout << " ** nrows[" << nrows << "] ncols[" << ncols << "]"
                        << endl;
            Double_t *allvalues = new Double_t[nrows * ncols];
            Double_t *dvalues = new Double_t[nrows];
            allvalues = lastmat->GetMatrixArray();
            //     TMatrixT<double> *matOut = new TMatrixT<double>(nrows, 1);
            Int_t ncrt = 0;
            cout
                    << " **********************************************************"
                    << endl;
            Double_t dweight;
            Int_t nObjectives = (Int_t)(_objectives.size());
            for (Int_t indx = 0; indx < nObjectives; indx++)
            {
                dweight = _objectives[indx]->getWeight();
                if (_objectives[indx]->isActive())
                {
                    if (_code->getLog())
                        _objectives[indx]->printLog();
                    ncrt = indx;
                    for (Int_t indy = 0; indy < nrows; indy++)
                    {
                        dvalues[indy] = allvalues[ncrt];
                        ncrt++;
                        //               ncrt += nrows;
                        //               cout << " ** indy[" << indy << "] val[" << dvalues[indy] << "]" << endl;
                    }
                    //       matOut->Use(nrows, 1, dvalues);
                    Double_t dvalcrt = _objectives[indx]->getObjective(dvalues,
                            nrows);
                    dObjectives += dweight * dvalcrt;
                    cout << " ** Objective [" << indx + 1 << "/" << nObjectives
                            << "] name[" << _objectives[indx]->GetName()
                            << "] dweight[" << dweight << "] val[" << dvalcrt
                            << "] Sum [" << dObjectives << "]" << endl;
                }
            }
            cout
                    << " **********************************************************"
                    << endl;
            dObjectives /= (Double_t)(nObjectives);
            cout << " ** Objective [" << dObjectives << "] nObjectives["
                    << nObjectives << "]" << endl;
            cout
                    << " **********************************************************"
                    << endl;
 
            fval = dObjectives;
 
            //Fill the data
            _dinputValue[_ninput + 1] = fval;
            _tds->getTuple()->Fill(_dinputValue);
            if (niter % 10)
            {
                _tds->Draw(
                        Form("%s:%s", _sCriteria.Data(),
                                _tds->getIteratorName()), "", "lp");
                gPad->Modified();
                gPad->Update();
            }
        }
        return fval;
    }
 
    double Up() const
    {
        return 1.;
    }
 
    void clean();
 
    void addObjective(TObjective *tobj)
    {
        _objectives.push_back(tobj);
    }
    Int_t getNObjectives()
    {
        return (Int_t)(_objectives.size());
    }
 
private:
    int _ninput;
    double * _dinputValue;
    int _noutput;
    double * _doutputValue;
    URANIE::DataServer::TDataServer *_tds;               
    URANIE::Launcher::TCode *_code;                
    TList *_listOfDataOutput;    
    TString _sCriteria;      
 
    ClassDef(URANIE::Optimizer::TSumOfSquares, ID_OPTIMIZER)
    //Classe de
};
} // Fin du namespace Optimizer
} // Fin du namespace URANIE
#endif