XSM_MLP.cxx

#include "XSModel.hxx"
#include "XSM_MLP.hxx"
#include "CLASSLogger.hxx"
#include "CLASSMethod.hxx"
#include "CLASSReader.hxx"
#include "StringLine.hxx"

#include "TMVA/Reader.h"
#include "TMVA/Tools.h"
#include "TMVA/MethodCuts.h"

#include <TGraph.h>
#include <TString.h>
#include "TFile.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TStopwatch.h"

#include <dirent.h>
#include <errno.h>
#include <sstream>
#include <string>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>

//________________________________________________________________________
//
//        XSM_MLP
//
//
//
//
//________________________________________________________________________
XSM_MLP::XSM_MLP(string TMVA_Weight_Directory,string InformationFile):XSModel(new CLASSLogger("XSM_MLP.log"))
{
    DBGL
    fTMVAWeightFolder = TMVA_Weight_Directory;
    
    fInformationFile = fTMVAWeightFolder+InformationFile;
    
    DBGL
    GetMLPWeightFiles();
    DBGL
    
    INFO << "__A cross section interpolator using" << endl;
    INFO << "Multi Layer Perceptron has been define__" << endl;
    INFO << " \t His TMVA folder is : \" " << fTMVAWeightFolder << "\"" << endl;
    
    LoadKeyword();
    DBGL
    ReadNFO();
    DBGL
    BookTMVAReader(); 
    DBGL
    
}
//________________________________________________________________________
XSM_MLP::XSM_MLP(CLASSLogger* Log,string TMVA_Weight_Directory,string InformationFile):XSModel(Log)
{
    DBGL
    fTMVAWeightFolder = TMVA_Weight_Directory;

    fInformationFile = fTMVAWeightFolder+InformationFile;
    
    GetMLPWeightFiles();
    DBGL
    
    INFO << "__A cross section interpolator using" << endl;
    INFO << "Multi Layer Perceptron has been define__" << endl;
    INFO << " \t His TMVA folder is : \" " << fTMVAWeightFolder << "\"" << endl;
    
    DBGL
    LoadKeyword();
    DBGL
    ReadNFO();
    DBGL
    BookTMVAReader(); 
    DBGL
}
//________________________________________________________________________
XSM_MLP::~XSM_MLP()
{
    DBGL
    fMapOfTMVAVariableNames.clear();
    fDKeyword.clear();
    //delete fReader;
    DBGL
}


//________________________________________________________________________
void XSM_MLP::LoadKeyword()
{
    DBGL
    XSModel::LoadKeyword();
    fDKeyword.insert( pair<string, XS_MLP_DMthPtr>( "k_timestep",    & XSM_MLP::ReadTimeSteps));
    fDKeyword.insert( pair<string, XS_MLP_DMthPtr>( "k_zainame",    & XSM_MLP::ReadZAIName)     );
    DBGL
}


//________________________________________________________________________
void XSM_MLP::ReadZAIName(const string &line)
{
    DBGL
    int pos = 0;
    string keyword = tlc(StringLine::NextWord(line, pos, ' '));
    if( keyword != "k_zainame" )    // Check the keyword
    {
        ERROR << " Bad keyword : \"k_zainame\" not found !" << endl;
        exit(1);
    }
    
    int Z = atoi(StringLine::NextWord(line, pos, ' ').c_str());
    int A = atoi(StringLine::NextWord(line, pos, ' ').c_str());
    int I = atoi(StringLine::NextWord(line, pos, ' ').c_str());
    
    string name = StringLine::NextWord(line, pos, ' ');
    
    fMapOfTMVAVariableNames.insert( pair<ZAI,string>( ZAI(Z, A, I), name ) );
    
    DBGL
}

//________________________________________________________________________
void XSM_MLP::ReadTimeSteps(const string &line)
{
    DBGL
    int pos = 0;
    string keyword = tlc(StringLine::NextWord(line, pos, ' '));
    if( keyword != "k_timestep" )    // Check the keyword
    {
        ERROR << " Bad keyword : \"k_timestep\" not found !" << endl;
        exit(1);
    }
    
    while( pos < (int)line.size() )
        fMLP_Time.push_back( atof( (StringLine::NextWord(line,pos,' ')).c_str() ));
    DBGL
}

//________________________________________________________________________
void XSM_MLP::ReadLine(string line)
{
    DBGL
    
    int pos = 0;
    string keyword = tlc(StringLine::NextWord(line, pos, ' '));
    
    map<string, XS_MLP_DMthPtr>::iterator it = fDKeyword.find(keyword);
    
    if(it != fDKeyword.end())
        (this->*(it->second))( line );
    
    DBGL
}


 
//________________________________________________________________________
void XSM_MLP::GetMLPWeightFiles()
{
    DBGL
    /**********Get All TMVA weight files*******************/
    //check if the folder containing weights exists
    DIR* rep = NULL;
    struct dirent* fichierLu = NULL;
    rep = opendir(fTMVAWeightFolder.c_str());
    
    if (rep ==  NULL)
    {
        ERROR << " Reading error for TMVA weight folder  " << fTMVAWeightFolder.c_str() << " : " << strerror(errno) << endl;
        exit(1);
    }
    
    /**Save file names of TMVA  weights*/
    fWeightFiles.resize(0);
    while ((fichierLu = readdir(rep)) != NULL)
    {
        string FileName = fichierLu->d_name ;
        if(FileName != "." && FileName != ".." )
        {
            if(FileName[FileName.size()-3] == 'x'  &&  FileName[FileName.size()-2] == 'm' && FileName[FileName.size()-1] == 'l' && FileName[0] != '.' )
                fWeightFiles.push_back(FileName);
            
        }
    }
    DBGL
}
//________________________________________________________________________
//________________________________________________________________________
//
//    Time  (MLP take time as parameter)
//________________________________________________________________________
//________________________________________________________________________
void XSM_MLP::ReadWeightFile(string Filename, int &Z, int &A, int &I, int &Reaction)
{
    Z = -1;
    A = -1;
    I = -1;
    Reaction = -1;
    
    size_t found = Filename.find("XS");
    
    string NameJOB;
    NameJOB = Filename.substr(found);
    int pos = 0;
    
    StringLine::NextWord(NameJOB, pos, '_');
    
    Z = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
    
    A = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
    
    I = atof( (StringLine::NextWord(NameJOB,pos,'_') ).c_str() );
    
    string  sReaction = (StringLine::NextWord(NameJOB,pos,'_') ).c_str() ;
    size_t foundext = sReaction.find(".weights.xml");
    sReaction = sReaction.substr(0,foundext);
    
    if(sReaction == "fis")
        Reaction = 0;
    if(sReaction == "cap")
        Reaction = 1;
    if(sReaction == "n2n")
        Reaction = 2;
    
    if(Z<= 0 || A<= 0 || I<0 || Reaction == -1)
    {
        ERROR << " wrong TMVA weight format " << endl;
        exit(0);
    }
    
}

void XSM_MLP::BookTMVAReader() {

    DBGL
    string dir = fTMVAWeightFolder;
    if (dir[dir.size() - 1] != '/') {
        dir += "/";
    }

    for (int i = 0; i < int(fWeightFiles.size()); i++) {
        int Z = -2;
        int A = -2;
        int I = -2;
        int Reaction = -2;
        ReadWeightFile(fWeightFiles[i], Z, A, I, Reaction);
        if (Z >= GetZAIThreshold()) {
            fReader.push_back(new CLASSReader(fMapOfTMVAVariableNames));
            fReader[i]->AddVariable("Time");
            fReader[i]->BookMVA("MLP method", dir + fWeightFiles[i]);
        }
    }
    DBGL
}

//________________________________________________________________________
vector<float> XSM_MLP::CreateTMVAInput(IsotopicVector isotopicvector,int TimeStep)
{
    DBGL
    vector<float> tmva_input;

    float Time = 0;

    IsotopicVector IVInputTMVA;
    map<ZAI , string >::iterator it_ZAI_s;

    for ( it_ZAI_s = fMapOfTMVAVariableNames.begin()  ; it_ZAI_s != fMapOfTMVAVariableNames.end() ; it_ZAI_s++)
    {
        IVInputTMVA +=  ((*it_ZAI_s).first) * 1;
    }
    
    // build IV containing only the ZAI known by the Model
    IsotopicVector IVAccordingToUserInfoFile    = isotopicvector.GetThisComposition(IVInputTMVA);
    
    // Normalize the vector
    double Ntot                     = IVAccordingToUserInfoFile.GetSumOfAll();
    IVAccordingToUserInfoFile           = IVAccordingToUserInfoFile / Ntot;

    // Add value in the input vector
    for ( it_ZAI_s = fMapOfTMVAVariableNames.begin() ; it_ZAI_s != fMapOfTMVAVariableNames.end() ; it_ZAI_s++)
    {
        tmva_input.push_back( IVAccordingToUserInfoFile.GetZAIIsotopicQuantity( (*it_ZAI_s).first ) );
    }

    tmva_input.push_back(fMLP_Time[TimeStep]);

    return tmva_input;
}
//________________________________________________________________________
EvolutionData XSM_MLP::GetCrossSections(IsotopicVector IV, double t) {
    DBGL

    string dir = fTMVAWeightFolder;
    if (dir[dir.size() - 1] != '/') {
        dir += "/";
    }

    EvolutionData EvolutionDataFromMLP = EvolutionData();

    map<ZAI, TGraph*> ExtrapolatedXS[3];
    /*************DATA BASE INFO****************/
    EvolutionDataFromMLP.SetReactorType(fDBRType);
    EvolutionDataFromMLP.SetFuelType(fDBFType);
    EvolutionDataFromMLP.SetPower(fDBPower);
    EvolutionDataFromMLP.SetHeavyMetalMass(fDBHMMass);
    /************* The Cross sections***********/
    int Z[int(fWeightFiles.size())];
    int A[int(fWeightFiles.size())];
    int I[int(fWeightFiles.size())];
    int Reaction[int(fWeightFiles.size())];

    for (int i = 0; i < int(fWeightFiles.size()); i++) {
        Z[i] = -2;
        A[i] = -2;
        I[i] = -2;
        Reaction[i] = -2;
        ReadWeightFile(fWeightFiles[i], Z[i], A[i], I[i], Reaction[i]);
    }
    
    vector<float> tmva_input = CreateTMVAInput(IV, 0);
    for (int TimeStep = 0; TimeStep < int(fMLP_Time.size()); TimeStep++) {
      tmva_input.back() = fMLP_Time[TimeStep];  
      for (int i = 0; i < int(fWeightFiles.size()); i++) {
            if (Z[i] >= GetZAIThreshold()) {
                fReader[i]->SetInputData(tmva_input);
                double XSValue = fReader[i]->EvaluateRegression("MLP method")[0];
                // std::cout << "XSValue " << XSValue << std::endl;
                pair<map<ZAI, TGraph*>::iterator, bool> IResult =
                        ExtrapolatedXS[Reaction[i]].insert(
                                pair<ZAI, TGraph*>(ZAI(Z[i], A[i], I[i]), new TGraph()));

                if (IResult.second) {
                    (IResult.first)
                            ->second->SetPoint(0, (double)fMLP_Time[TimeStep], XSValue);

                } else {
                    (IResult.first)
                            ->second->SetPoint((IResult.first)->second->GetN(),
                                                                 (double)fMLP_Time[TimeStep], XSValue);
                }
            }
        }
    }

    /**********Sorting TGraph*********/
    for (int x = 0; x < 3; x++) {
        map<ZAI, TGraph*>::iterator it;
        for (it = ExtrapolatedXS[x].begin(); it != ExtrapolatedXS[x].end(); it++)
            it->second->Sort();
    }
    /**********Filling Matrices*/
    EvolutionDataFromMLP.SetFissionXS(ExtrapolatedXS[0]);
    EvolutionDataFromMLP.SetCaptureXS(ExtrapolatedXS[1]);
    EvolutionDataFromMLP.Setn2nXS(ExtrapolatedXS[2]);

    DBGL return EvolutionDataFromMLP;
}

//________________________________________________________________________