diff --git a/.gitignore b/.gitignore index 310549d494eaf9e218072961ed141750332dd7e2..0a91ea85f3ab64b83b39056adee7240e4a9c8367 100644 --- a/.gitignore +++ b/.gitignore @@ -43,6 +43,7 @@ BUILD *.app *.pcm CLASS_R2D +CLASS_R2D_S # ROOT files *.root @@ -64,7 +65,7 @@ source/include/Irradiation config # folder containing the GUI binary : -gui/bin +bin # Built doxygen : documentation/doxygen/html/doxygen diff --git a/GTest/Makefile b/GTest/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..9ac74493ba55246597e90b898fae847ae358eab5 --- /dev/null +++ b/GTest/Makefile @@ -0,0 +1,82 @@ +# A sample Makefile for building Google Test and using it in user +# tests. Please tweak it to suit your environment and project. You +# may want to move it to your project's root directory. +# +# SYNOPSIS: +# +# make [all] - makes everything. +# make TARGET - makes the given target. +# make clean - removes all files generated by make. + +# Please tweak the following variable definitions as needed by your +# project, except GTEST_HEADERS, which you can use in your own targets +# but shouldn't modify. + +# Points to the root of Google Test, relative to where this file is. +# Remember to tweak this if you move this file. +GTEST_DIR = .. + +# Where to find user code. +USER_DIR = ../samples + +# Flags passed to the preprocessor. +# Set Google Test's header directory as a system directory, such that +# the compiler doesn't generate warnings in Google Test headers. +CPPFLAGS += -isystem $(GTEST_DIR)/include + +# Flags passed to the C++ compiler. +CXXFLAGS += -g -Wall -Wextra -pthread + +# All tests produced by this Makefile. Remember to add new tests you +# created to the list. +TESTS = sample1_unittest + +# All Google Test headers. Usually you shouldn't change this +# definition. +GTEST_HEADERS = $(GTEST_DIR)/include/gtest/*.h \ + $(GTEST_DIR)/include/gtest/internal/*.h + +# House-keeping build targets. + +all : $(TESTS) + +clean : + rm -f $(TESTS) gtest.a gtest_main.a *.o + +# Builds gtest.a and gtest_main.a. + +# Usually you shouldn't tweak such internal variables, indicated by a +# trailing _. +GTEST_SRCS_ = $(GTEST_DIR)/src/*.cc $(GTEST_DIR)/src/*.h $(GTEST_HEADERS) + +# For simplicity and to avoid depending on Google Test's +# implementation details, the dependencies specified below are +# conservative and not optimized. This is fine as Google Test +# compiles fast and for ordinary users its source rarely changes. +gtest-all.o : $(GTEST_SRCS_) + $(CXX) $(CPPFLAGS) -I$(GTEST_DIR) $(CXXFLAGS) -c \ + $(GTEST_DIR)/src/gtest-all.cc + +gtest_main.o : $(GTEST_SRCS_) + $(CXX) $(CPPFLAGS) -I$(GTEST_DIR) $(CXXFLAGS) -c \ + $(GTEST_DIR)/src/gtest_main.cc + +gtest.a : gtest-all.o + $(AR) $(ARFLAGS) $@ $^ + +gtest_main.a : gtest-all.o gtest_main.o + $(AR) $(ARFLAGS) $@ $^ + +# Builds a sample test. A test should link with either gtest.a or +# gtest_main.a, depending on whether it defines its own main() +# function. + +sample1.o : $(USER_DIR)/sample1.cc $(USER_DIR)/sample1.h $(GTEST_HEADERS) + $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $(USER_DIR)/sample1.cc + +sample1_unittest.o : $(USER_DIR)/sample1_unittest.cc \ + $(USER_DIR)/sample1.h $(GTEST_HEADERS) + $(CXX) $(CPPFLAGS) $(CXXFLAGS) -c $(USER_DIR)/sample1_unittest.cc + +sample1_unittest : sample1.o sample1_unittest.o gtest_main.a + $(CXX) $(CPPFLAGS) $(CXXFLAGS) -lpthread $^ -o $@ diff --git a/GTest/main_test.cxx b/GTest/main_test.cxx new file mode 100644 index 0000000000000000000000000000000000000000..080077fce7289f5e0a3ee4dad1ae1531aa406e67 --- /dev/null +++ b/GTest/main_test.cxx @@ -0,0 +1,17 @@ +#include <gtest/gtest.h> +#include "test_iv.inl" + +int main(int argc,char * argv[]) { + ::testing::InitGoogleTest(&argc,argv); + return RUN_ALL_TESTS(); +} + +/* +Google test library compilation + +g++ -isystem ${GTEST_DIR}/include -I${GTEST_DIR} -pthread -c ${GTEST_DIR}/src/gtest-all.cc +ar -rv libgtest.a gtest-all.o + +Main Test compilation +g++ -isystem ${GTEST_DIR}/include -pthread main_test.cxx libgtest.a -o MyTest -I${GTEST_DIR}/include -I$CLASS_include -I$CLASS_external -I$CLASS_Equivalence -I$CLASS_Irradiation -I$CLASS_XS -L$CLASS_lib -lCLASSpkg `root-config --cflags` `root-config --libs` +*/ diff --git a/test/test_iv.inl b/GTest/test_iv.inl similarity index 54% rename from test/test_iv.inl rename to GTest/test_iv.inl index b2fbe1e40d0ae678fcd642fd87bbf9c178b119a0..62c50ef798ffbdcd3590415bef56e02b99db5a60 100644 --- a/test/test_iv.inl +++ b/GTest/test_iv.inl @@ -3,9 +3,9 @@ #include <iomanip> #include <math.h> #include <string> -#include "XS/XSM_MLP.hxx" //Load the include for Neural network cross section predictor -#include "Irradiation/IM_RK4.hxx" //Load the include for Runge Kutta 4 resolution -#include "Equivalence/EQM_MLP_Kinf.hxx" //Load the include for Neural Network Equivalence Model (PWRMOX) +#include "XSM_MLP.hxx" //Load the include for Neural network cross section predictor +#include "IM_RK4.hxx" //Load the include for Runge Kutta 4 resolution +#include "EQ_OneParameter.hxx" //Load the include for Neural Network Equivalence Model (PWRMOX) TEST ( TestIV, getSize ) { const int n=10; diff --git a/Utils/ROOT2DAT/CLASS_ROOT2DAT.cxx b/Utils/ROOT2DAT/CLASS_ROOT2DAT.cxx index 3d1048d1f3d5b339a6173d3776bc3cdb403a7718..9e34b4729f61a21788c779d09f58d5a938563023 100644 --- a/Utils/ROOT2DAT/CLASS_ROOT2DAT.cxx +++ b/Utils/ROOT2DAT/CLASS_ROOT2DAT.cxx @@ -276,5 +276,5 @@ DataFileName<<"C"<<endl; } /* - g++ -o CLASS_R2D CLASS_ROOT2DAT.cxx -I $CLASS_include -L $CLASS_lib -lCLASSpkg `root-config --cflags` `root-config --libs` -fopenmp -lgomp -Wunused-result +g++ -o CLASS_R2D CLASS_ROOT2DAT.cxx -I$CLASS_PATH/source/include -I$CLASS_PATH/source/external -I$CLASS_PATH/source/Model/Equivalence -I$CLASS_PATH/source/Model/Irradiation -I$CLASS_PATH/source/Model/XS -L$CLASS_lib -lCLASSpkg `root-config --cflags` `root-config --libs` */ diff --git a/Utils/ROOT2DAT/CLASS_ROOT2DAT_Simple.cxx b/Utils/ROOT2DAT/CLASS_ROOT2DAT_Simple.cxx index 847b971e7f9163fc88d7811544889687cd71d855..95973b3e593001c840a903936c23b745561a9bb4 100644 --- a/Utils/ROOT2DAT/CLASS_ROOT2DAT_Simple.cxx +++ b/Utils/ROOT2DAT/CLASS_ROOT2DAT_Simple.cxx @@ -250,5 +250,5 @@ DataFileName<<"C"<<endl; } /* - g++ -o CLASS_R2D_S CLASS_ROOT2DAT_Simple.cxx -I $CLASS_include -L $CLASS_lib -lCLASSpkg `root-config --cflags` `root-config --libs` -fopenmp -lgomp -Wunused-result +g++ -o CLASS_R2D CLASS_ROOT2DAT.cxx -I$CLASS_PATH/source/include -I$CLASS_PATH/source/external -I$CLASS_PATH/source/Model/Equivalence -I$CLASS_PATH/source/Model/Irradiation -I$CLASS_PATH/source/Model/XS -L$CLASS_lib -lCLASSpkg `root-config --cflags` `root-config --libs` */ diff --git a/source/Model/Equivalence/CMakeLists.txt b/source/Model/Equivalence/CMakeLists.txt index 0d37dde2384a9c9cdd2c134231b7d39d8dbcd5a4..09f9043af93c1ea2b62ac332f5d38a73cf8214d9 100644 --- a/source/Model/Equivalence/CMakeLists.txt +++ b/source/Model/Equivalence/CMakeLists.txt @@ -4,3 +4,9 @@ # include directories SET( CLASS_HEADERS_DIR_EQUIVALENCE ${CMAKE_CURRENT_SOURCE_DIR} CACHE INTERNAL "") + +LIST( APPEND CLASS_SRC_FILE_EQUIVALENCE "${CMAKE_CURRENT_SOURCE_DIR}/EQ_OneParameter.cxx") + +SET( CLASS_SRC_FILE_IRRADIATION ${CLASS_SRC_FILE_EQUIVALENCE} PARENT_SCOPE) + +MESSAGE( STATUS "CLASS_SRC_FILE_EQUIVALENCE=" ${CLASS_SRC_FILE_EQUIVALENCE}) diff --git a/source/Model/Equivalence/EQ_OneParameter.cxx b/source/Model/Equivalence/EQ_OneParameter.cxx new file mode 100644 index 0000000000000000000000000000000000000000..a65a83177aa0842cebc697940c6027c6256a9900 --- /dev/null +++ b/source/Model/Equivalence/EQ_OneParameter.cxx @@ -0,0 +1,1181 @@ +#include "EquivalenceModel.hxx" +#include "EQ_OneParameter.hxx" +#include "StringLine.hxx" +#include "CLASSMethod.hxx" + +#include <string> +#include <iostream> +#include <fstream> +#include <algorithm> +#include <cmath> +#include <cassert> + +#include "TSystem.h" +#include "TMVA/Reader.h" +#include "TMVA/Tools.h" +#include "TMVA/MethodCuts.h" + +#include "CLASSReader.hxx" + +//________________________________________________________________________ +//________________________________________________________________________ +EQ_OneParameter::EQ_OneParameter(string TMVAXMLFilePath, string TMVANFOFilePath):EquivalenceModel(new CLASSLogger("EQ_OneParameter.log")) +{ + fUseTMVAPredictor = true; + + fMaxIterration = 500; + fPCMprecision = 10; + + fTMVAXMLFilePath = TMVAXMLFilePath; + fTMVANFOFilePath = TMVANFOFilePath; + + fDBRType = ""; + fDBFType = ""; + fSpecificPower = 0; + fMaximalBU = 0; + fTargetParameter = ""; + fTargetParameterStDev = 0; + fBuffer = ""; + fPredictorType = ""; + fOutput = ""; + + LoadKeyword(); // Load Key words defineds in NFO file + ReadNFO(); //Getting information from file NFO + + //Check if any information is missing in NFO file + if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); } + if(fMapOfTMVAVariableNames.empty()) { ERROR<<"Missing information for k_zainame in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fTargetParameter.empty()) { ERROR<<"Missing information for k_targetparameter in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(!fTargetParameterStDev) { ERROR<<"Missing information for fTargetParameterStDev in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fModelParameter.empty()) { ERROR<<"Missing information for k_modelparameter in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fPredictorType.empty()) { ERROR<<"Missing information for k_predictortype in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fOutput.empty()) { ERROR<<"Missing information for k_output in : "<<fTMVANFOFilePath<<endl; exit(1);} + + INFO << "__An equivalence model has been define__" << endl; + INFO << "\tThe TMVA weights file is :" << fTMVAXMLFilePath << endl; + INFO << "\tThe TMVA NFO file is :" << fTMVANFOFilePath << endl; + PrintInfo(); + +} +//________________________________________________________________________ +EQ_OneParameter::EQ_OneParameter(CLASSLogger* log, string TMVAXMLFilePath, string TMVANFOFilePath):EquivalenceModel(log) +{ + fUseTMVAPredictor = true; + + fMaxIterration = 500; + freaded = false; + fPCMprecision = 10; + + fTMVAXMLFilePath = TMVAXMLFilePath; + fTMVANFOFilePath = TMVANFOFilePath; + + fDBRType = ""; + fDBFType = ""; + fSpecificPower = 0; + fMaximalBU = 0; + fTargetParameter = ""; + fTargetParameterStDev = 0; + fBuffer = ""; + fPredictorType = ""; + fOutput = ""; + + LoadKeyword(); // Load Key words defineds in NFO file + ReadNFO(); //Getting information from file NFO + + //Check if any information is missing in NFO file + if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); } + if(fMapOfTMVAVariableNames.empty()) { ERROR<<"Missing information for k_zainame in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fTargetParameter.empty()) { ERROR<<"Missing information for k_targetparameter in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(!fTargetParameterStDev) { ERROR<<"Missing information for fTargetParameterStDev in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fModelParameter.empty()) { ERROR<<"Missing information for k_modelparameter in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fPredictorType.empty()) { ERROR<<"Missing information for k_predictortype in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fOutput.empty()) { ERROR<<"Missing information for k_output in : "<<fTMVANFOFilePath<<endl; exit(1);} + + INFO << "__An equivalence model has been define__" << endl; + INFO << "\tThe TMVA weights file is :" << fTMVAXMLFilePath << endl; + INFO << "\tThe TMVA NFO file is :" << fTMVANFOFilePath << endl; + PrintInfo();} +//________________________________________________________________________ +EQ_OneParameter::EQ_OneParameter(string TMVANFOFilePath):EquivalenceModel(new CLASSLogger("EQ_OneParameter.log")) +{ + fUseTMVAPredictor = false; + + fTMVANFOFilePath = TMVANFOFilePath; + + fDBRType = ""; + fDBFType = ""; + fSpecificPower = 0; + fMaximalBU = 0; + fTargetParameter = ""; + fTargetParameterStDev = 0; + fBuffer = ""; + + LoadKeyword(); // Load Key words defineds in NFO file + ReadNFO(); //Getting information from file NFO + + //Check if any information is missing in NFO file + if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); } + if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);} + + INFO << "__An equivalence model without TMVA data has been define__" << endl; + INFO << "\tThe NFO file is :" << fTMVANFOFilePath << endl; + PrintInfo(); +} +//________________________________________________________________________ +EQ_OneParameter::EQ_OneParameter(CLASSLogger* log, string TMVANFOFilePath):EquivalenceModel(log) +{ + fUseTMVAPredictor = false; + + fTMVANFOFilePath = TMVANFOFilePath; + + fDBRType = ""; + fDBFType = ""; + fSpecificPower = 0; + fMaximalBU = 0; + fBuffer = ""; + + LoadKeyword(); // Load Key words defineds in NFO file + ReadNFO(); //Getting information from file NFO + + //Check if any information is missing in NFO file + if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);} + if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); } + if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);} + + INFO << "__An equivalence model has been define__" << endl; + INFO << "\tThe TMVA weights file is :" << fTMVAXMLFilePath << endl; + INFO << "\tThe TMVA NFO file is :" << fTMVANFOFilePath << endl; + PrintInfo();} +//________________________________________________________________________ +EQ_OneParameter::~EQ_OneParameter() +{ + +} +//________________________________________________________________________ +IsotopicVector EQ_OneParameter::BuildFuelToTest(map < string, vector<double> >& lambda, map < string , vector <IsotopicVector> > const& StreamArray, double HMMass, map <string, bool> StreamListIsBuffer) +{ + DBGL + //Iterators declaration + map < string , vector <IsotopicVector> >::const_iterator it_s_vIV; + map < string , bool >::iterator it_s_B; + + //Find the buffer and set its lambda to 0 + string BufferDenomination =""; + for( it_s_B = StreamListIsBuffer.begin(); it_s_B != StreamListIsBuffer.end(); it_s_B++) + { + if((*it_s_B ).second==true){ BufferDenomination = (*it_s_B).first; } + } + + for(int i = 0; i< lambda[BufferDenomination].size(); i++) + { + lambda[BufferDenomination][i]=0; + } + + //Build an IV with all materials besides buffer to get the total mass of others materials + IsotopicVector IV; + for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++) + { + for(int i=0; i < (int)StreamArray.at( it_s_vIV->first ).size(); i++) + { + IV += lambda[(*it_s_vIV).first][i] * StreamArray.at( it_s_vIV->first )[i]; + } + } + + //Calculate MassBuffer + double MassBuffer = HMMass - IV.GetTotalMass()*1e06; + + //Set buffer lambda according to MassBuffer + ConvertMassToLambdaVector(BufferDenomination, lambda[BufferDenomination], MassBuffer, StreamArray.at(BufferDenomination)); + + IV.Clear(); + + //Build fuel with all materials + for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++) + { + for(int i=0; i < (int)StreamArray.at( it_s_vIV->first ).size(); i++) + { + IV += lambda[(*it_s_vIV).first][i] * StreamArray.at( it_s_vIV->first )[i]; + } + } + DBGL + return IV; + +} + +//________________________________________________________________________ +map <string , vector<double> > EQ_OneParameter::BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray, map < string , double> StreamListFPMassFractionMin, map < string , double> StreamListFPMassFractionMax, map < int , string > StreamListPriority, map < string , bool> StreamListIsBuffer) +{ + DBGL + + HMMass *= 1e6; //Unit conversion : tons to gram + + map <string , vector<double> > lambda ; // map containing name of the list and associated vector of proportions taken from stocks + //Iterators declaration + map < string , vector <IsotopicVector> >::iterator it_s_vIV; + map < string , vector <double> >::iterator it_s_vD; + map < string , IsotopicVector >::iterator it_s_IV; + map < string , double >::iterator it_s_D; + map < int , string >::iterator it_i_s; + + // Initialize lambda to 0 // + for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++) + { + for(int i=0; i < (int)StreamArray[(*it_s_vIV).first].size(); i++) + { + lambda[(*it_s_vIV).first].push_back(0); + } + } + + // Test if there is at least one stock available in each list, otherwise fuel is not built // + bool BreakReturnLambda = false; + for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++) + { + if(StreamArray[(*it_s_vIV).first].size() == 0) + { + WARNING << " No stock available for stream : "<< (*it_s_vIV).first <<". Fuel not built." << endl; + SetLambdaToErrorCode(lambda[(*it_s_vIV).first]); + BreakReturnLambda = true; + } + } + if(BreakReturnLambda) { return lambda;} + + // Check if the targeted burn-up is lower than maximum burn-up of model // + if(BurnUp > fMaximalBU) + { + ERROR << " Targeted burn-up is higher than maximum burn-up defined in NFO file..."<< endl; + ERROR << " Targeted burn-up : "<<BurnUp<<" GWd/t"<<endl; + ERROR << " Maximum burn-up : "<<fMaximalBU<<" GWd/t"<<endl; + exit(1); + } + +// Fissile fraction calculation is needed. + if (fUseTMVAPredictor) + { + // Check if EQ_OneParameter->SetTMVAXMLFilePath() and/or EQ_OneParameter->SetTMVANFOFilePath() have been defined + if (fTMVAXMLFilePath.empty() || fTMVANFOFilePath.empty()) + { + ERROR << " TMVA XML and/or NFO File path are not defined..."<< endl; + ERROR << " You have to use EQ_OneParameter->SetTMVAXMLFilePath() and/or EQ_OneParameter->SetTMVANFOFilePath() methods."<<endl; + exit(1); + } + + double TargetParameterValue = 0; + + for( it_s_D = fModelParameter.begin(); it_s_D != fModelParameter.end(); it_s_D++) + { + if(fModelParameter[(*it_s_D).first] == -1) + { + ERROR<< "Model parameter ( "<<fModelParameter[(*it_s_D).first] << " ) value is not defined in the input." <<endl; + ERROR<< "Use EqM->SetModelParameter( \" "<<(*it_s_D).first<<" \", value) to define it." <<endl; + exit(1); + } + } + + if(fTargetParameter=="BurnUpMax") {TargetParameterValue = BurnUp;} + else if (fTargetParameter=="keffBOC") {TargetParameterValue = fModelParameter["keffBOC"];} + else + { + ERROR<< "Target parameter defined in InformationFile ( "<<fTargetParameter<<" ) doesn't exist." <<endl; + ERROR<< "Possible target parameters for the moment are : "<< endl; + ERROR<< " - BurnUpMax - Used for PWR" <<endl; + ERROR<< " - keffBOC - Used for SFR" <<endl; + exit(1); + } + + /// Search for the minimum and maximum fraction of each material in fuel /// + map < string, double > StreamListMassFractionMin ; + map < string, double > StreamListMassFractionMax ; + for( it_s_D = StreamListFPMassFractionMin.begin(); it_s_D != StreamListFPMassFractionMin.end(); it_s_D++) + { + if(StreamListFPMassFractionMin[(*it_s_D).first] < fStreamListEqMMassFractionMin[(*it_s_D).first]) // if limits FP are lower than limits EqM + { + ERROR << " User mass fraction min requirement is lower than the model mass fraction min for list : "<<(*it_s_D).first << endl; + ERROR << " User mass fraction min requirement : "<<StreamListFPMassFractionMin[(*it_s_D).first]<<endl; + ERROR << " Model mass fraction min requirement : "<<fStreamListEqMMassFractionMin[(*it_s_D).first]<<endl; + exit(1); + } + else + { + StreamListMassFractionMin[(*it_s_D).first] = StreamListFPMassFractionMin[(*it_s_D).first]; + } + } + + for( it_s_D = StreamListFPMassFractionMax.begin(); it_s_D != StreamListFPMassFractionMax.end(); it_s_D++) + { + if(StreamListFPMassFractionMax[(*it_s_D).first] > fStreamListEqMMassFractionMax[(*it_s_D).first]) // if limits FP are higher than limits EqM + { + ERROR << " User mass fraction max requirement is higher than the model mass fraction max for list : "<<(*it_s_D).first << endl; + ERROR << " User mass fraction max requirement : "<<StreamListFPMassFractionMax[(*it_s_D).first]<<endl; + ERROR << " Model mass fraction max requirement : "<<fStreamListEqMMassFractionMax[(*it_s_D).first]<<endl; + exit(1); + } + else + { + StreamListMassFractionMax[(*it_s_D).first] = StreamListFPMassFractionMax[(*it_s_D).first]; + } + + } + + //Calculate Total mass in stock for each stream and fill fTotalMassInStocks + StocksTotalMassCalculation(StreamArray); + + // Check if there is enough material in stock to satisfy mass fraction min // + BreakReturnLambda = false; + for( it_s_D = StreamListMassFractionMin.begin(); it_s_D != StreamListMassFractionMin.end(); it_s_D++) + { + if(fTotalMassInStocks[(*it_s_D).first]< HMMass*StreamListMassFractionMin[(*it_s_D).first]) + { + WARNING << " Not enough material : "<< (*it_s_D).first << " in stocks to reach the build fuel lower limit of "<<StreamListMassFractionMin[(*it_s_D).first]<<" reactor mass. Fuel not built." << endl; + SetLambdaToErrorCode(lambda[(*it_s_D).first]); + BreakReturnLambda = true; + } + } + if(BreakReturnLambda) { return lambda;} + + // Check if there is enough material in stock to satisfy mass fraction max, if not mass fraction max is set to MassINStock/MassReactor// + for( it_s_D = StreamListMassFractionMax.begin(); it_s_D != StreamListMassFractionMax.end(); it_s_D++) + { + if(fTotalMassInStocks[(*it_s_D).first]< HMMass*StreamListMassFractionMax[(*it_s_D).first]) + { + StreamListMassFractionMax[(*it_s_D).first] = fTotalMassInStocks[(*it_s_D).first]/HMMass; + WARNING << " Not enough material : "<< (*it_s_D).first << " in stocks to reach the build fuel higher limit of "<<StreamListMassFractionMax[(*it_s_D).first]<<" reactor mass. " << endl; + WARNING << " Mass fraction max of material : "<< (*it_s_D).first << " is set to MassInStock/HMMassReactor : "<< StreamListMassFractionMax[(*it_s_D).first]<< endl; + } + } + + //Check if TargetParameter is inside [TargetParameterMin, TargetParameterMax] associated to fraction Min et Max// + + map < string , double > MassMin; + map < string , double > MassMax; + + map < string , double > TargetParameterMin; + map < string , double > TargetParameterMax; + + IsotopicVector FuelToTest; + + bool TargetParameterIncluded = false; + for( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++) + { + //Calculate TargetParameterMin for each possibility : min1 ; max1 + min2 ; max1 + max2 + min3 .... + MassMin[(*it_i_s ).second] = HMMass * StreamListMassFractionMin[(*it_i_s).second]; + ConvertMassToLambdaVector((*it_i_s ).second, lambda[(*it_i_s ).second], MassMin[(*it_i_s ).second], StreamArray[(*it_i_s ).second]); + FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer); + FuelToTest = FuelToTest/FuelToTest.GetSumOfAll(); + TargetParameterMin[(*it_i_s ).second] = CalculateTargetParameter(FuelToTest, fTargetParameter); + + //Check is TargetParameterMin < TargetParameter + if(TargetParameterMin[(*it_i_s ).second]>TargetParameterValue) + { + if((*it_i_s).first ==1) //Minimum of first material is too high + { + WARNING << "CRITICAL ! Minimum parameter value associated to the first priority material ( "<<(*it_i_s ).second <<" ) is higher than targeted parameter."<< endl; + WARNING << "Targeted parameter : "<<fTargetParameter<<" = "<<TargetParameterValue<<endl; + WARNING << "Minimum parameter value : " <<TargetParameterMin[(*it_i_s ).second]<<endl; + WARNING << "Try to increase targeted parameter." <<endl; + SetLambdaToErrorCode(lambda[(*it_i_s).second]); + return lambda; + DBGL + } + else if ((*it_i_s).first >1) //TargetParameter is located between max n-1 and min n + { + WARNING << "CRITICAL ! Targeted parameter value ( "<<fTargetParameter<<" ) is located between 2 materials. "<<endl; + it_i_s --; + WARNING << fTargetParameter <<" of max fraction of material : "<< (*it_i_s).second<<" ---> "<<TargetParameterMax[(*it_i_s ).second]<<endl; + it_i_s ++; + WARNING << fTargetParameter<< " of min fraction of material : "<< (*it_i_s ).second<<" ---> "<<TargetParameterMin[(*it_i_s ).second]<<endl; + WARNING << "Targeted "<<fTargetParameter<<" : " <<TargetParameterValue<<endl; + WARNING << "Try to decrease mimimum fraction of : "<< (*it_i_s ).second<<endl; + SetLambdaToErrorCode(lambda[(*it_i_s).second]); + return lambda; + } + } + FuelToTest.Clear(); + + //Calculate TargetParameter max for each possibility : max1 ; max1 + max2 ; max1 + max2 + max3 .... + MassMax[(*it_i_s ).second] = HMMass * StreamListMassFractionMax[(*it_i_s).second]; + ConvertMassToLambdaVector((*it_i_s ).second, lambda[(*it_i_s ).second], MassMax[(*it_i_s ).second], StreamArray[(*it_i_s ).second]); + FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer); + FuelToTest = FuelToTest/FuelToTest.GetSumOfAll(); + TargetParameterMax[(*it_i_s ).second] = CalculateTargetParameter(FuelToTest, fTargetParameter); + + if(TargetParameterMax[(*it_i_s ).second]>=TargetParameterValue) + { + TargetParameterIncluded = true ; + break; + } + } + + //Check if target parameter increases monotously with the material mass + CheckTargetParameterConsistency(StreamListPriority, TargetParameterMin, TargetParameterMax); + + if(!TargetParameterIncluded) + { + WARNING << "CRITICAL ! Maximum reachable "<<fTargetParameter<<" is lower than targeted "<< fTargetParameter<<". "<< endl; + WARNING << "Targeted "<<fTargetParameter<<" = "<<TargetParameterValue<<endl; + WARNING << "Maximum reachable "<<fTargetParameter<<" : "<<TargetParameterMax[(*--StreamListPriority.end()).second]<<endl; + WARNING << "Try to increase maximum fraction of materials, or decrease "<< fTargetParameter<<" ." <<endl; + SetLambdaToErrorCode(lambda[(*--StreamListPriority.end()).second]); + return lambda; + } + + //Search the TargetParameterValue location in the mass damain // + string MaterialToSearch = (*it_i_s ).second; + double CalculatedTargetParameter = TargetParameterMax[MaterialToSearch] ; //Algo start with maximum point + double MassToAdd = MassMax[MaterialToSearch]; //Algo start with maximum point + + double LastMassMinus = MassMin[MaterialToSearch]; //Used in bissection method + double LastMassPlus = MassMax[MaterialToSearch]; //Used in bissection method + + int count = 0; + + FuelToTest.Clear(); + + /* + if (fDBFType == "MOX") + { + cout<<"------------------------------------------------------"<<endl; + cout<<"START ALGO -> BU, Mass "<<BurnUp<<" "<<HMMass<<endl; + cout<<"------------------------------------------------------"<<endl; + double MassTest = MassMin[MaterialToSearch]; + cout<<MaterialToSearch<<" "<<MassMax[MaterialToSearch]<<" "<<MassMin[MaterialToSearch]<<" "<<endl; + do + { + ConvertMassToLambdaVector(MaterialToSearch, lambda[MaterialToSearch], MassTest, StreamArray[MaterialToSearch]); + FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer); + FuelToTest = FuelToTest/FuelToTest.GetSumOfAll(); + CalculatedTargetParameter = CalculateTargetParameter(FuelToTest, fTargetParameter); + + cout<<"Lambda vector : "<<MaterialToSearch<<" - "; for(int i=0; i < (int)lambda[MaterialToSearch].size(); i++) cout<<lambda[MaterialToSearch][i]<<" "; + cout<<endl; + + + MassTest += (MassMax[MaterialToSearch] - MassMin[MaterialToSearch])/100.; + + cout<<MassTest<<" "<<CalculatedTargetParameter<<endl; + + } while (MassTest <= MassMax[MaterialToSearch]); + cout<<"------------------------------------------------------"<<endl; + cout<<"STOP ALGO EXIT(1)..."<<endl; exit(1); + cout<<"------------------------------------------------------"<<endl; + } + */ + + do + { + if(count > fMaxIterration) + { + ERROR << "CRITICAL ! Can't manage to predict fissile content\nHint : Try to decrease the precision on the target parameter using :\nYourEQ_OneParameter->SetTargetParameterStDev(Precision); " << endl; + ERROR << "Targeted "<<fTargetParameter<<" : "<<TargetParameterValue<<endl; + ERROR << "Last calculated "<<fTargetParameter<<" : "<<CalculatedTargetParameter<<endl; + ERROR << "Last Fresh fuel normalized composition : " <<endl; + ERROR << FuelToTest.sPrint()<<endl; + exit(1); + } + + if( (CalculatedTargetParameter - TargetParameterValue) < 0 ) //Need to add more fissile material in fuel + { + LastMassMinus = MassToAdd; + MassToAdd = MassToAdd + fabs(LastMassPlus - MassToAdd)/2.; + } + else if( (CalculatedTargetParameter - TargetParameterValue) > 0) //Need to add less fissile material in fuel + { + LastMassPlus = MassToAdd; + MassToAdd = MassToAdd - fabs(LastMassMinus - MassToAdd)/2.; + } + ConvertMassToLambdaVector(MaterialToSearch, lambda[MaterialToSearch], MassToAdd, StreamArray[MaterialToSearch]); + FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer); + FuelToTest = FuelToTest/FuelToTest.GetSumOfAll(); + CalculatedTargetParameter = CalculateTargetParameter(FuelToTest, fTargetParameter); + + count ++; + + }while(fabs(TargetParameterValue - CalculatedTargetParameter) > GetTargetParameterStDev()*TargetParameterValue); + } +// Fissile fraction is imposed by the FP +// No need to use algo +// Simplified fuel building + else + { + // Check if EQ_OneParameter->SetTMVANFOFilePath() have been defined + if (fTMVANFOFilePath.empty()) + { + ERROR << " TMVA NFO File path is not defined..."<< endl; + ERROR << " You have to use EQ_OneParameter->SetTMVANFOFilePath() methods."<<endl; + exit(1); + } + + /// Search for the fraction of each material in fuel /// + map < string, double > StreamListMassFraction; + for( it_s_D = StreamListFPMassFractionMin.begin(); it_s_D != StreamListFPMassFractionMin.end(); it_s_D++) + { + if(StreamListFPMassFractionMin[(*it_s_D).first] < fStreamListEqMMassFractionMin[(*it_s_D).first]) // if limits FP are lower than limits EqM + { + ERROR << " User mass fraction requirement is lower than the model mass fraction min for list : "<<(*it_s_D).first << endl; + ERROR << " User mass fraction requirement : "<<StreamListFPMassFractionMin[(*it_s_D).first]<<endl; + ERROR << " Model mass fraction min requirement : "<<fStreamListEqMMassFractionMin[(*it_s_D).first]<<endl; + exit(1); + } + else if(StreamListFPMassFractionMax[(*it_s_D).first] > fStreamListEqMMassFractionMax[(*it_s_D).first]) // if limits FP are higher than limits EqM + { + ERROR << " User mass fraction requirement is higher than the model mass fraction max for list : "<<(*it_s_D).first << endl; + ERROR << " User mass fraction requirement : "<<StreamListFPMassFractionMax[(*it_s_D).first]<<endl; + ERROR << " Model mass fraction max requirement : "<<fStreamListEqMMassFractionMax[(*it_s_D).first]<<endl; + exit(1); + } + else + { + StreamListMassFraction[(*it_s_D).first] = StreamListFPMassFractionMin[(*it_s_D).first]; // Because here, min = max + } + } + + //Calculate Total mass in stock for each stream and fill fTotalMassInStocks + StocksTotalMassCalculation(StreamArray); + + // Check if there is enough material in stock to satisfy requested mass fraction // + BreakReturnLambda = false; + for( it_s_D = StreamListMassFraction.begin(); it_s_D != StreamListMassFraction.end(); it_s_D++) + { + if(fTotalMassInStocks[(*it_s_D).first]< HMMass*StreamListMassFraction[(*it_s_D).first]) + { + WARNING << " Not enough material : "<< (*it_s_D).first << " in stocks to reach the build fuel limit of "<<StreamListMassFraction[(*it_s_D).first]<<" reactor mass. Fuel not built." << endl; + SetLambdaToErrorCode(lambda[(*it_s_D).first]); + BreakReturnLambda = true; + } + } + if(BreakReturnLambda) { return lambda;} + + IsotopicVector FuelToTest; + // Build Fuel + for( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++) + { + FuelToTest.Clear(); + ConvertMassToLambdaVector((*it_i_s).second, lambda[(*it_i_s).second], HMMass*StreamListMassFraction[(*it_i_s).second], StreamArray[(*it_i_s).second]); + FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer); + } + } + + //Final builded fuel + IsotopicVector IVStream; + for( it_s_vD = lambda.begin(); it_s_vD != lambda.end(); it_s_vD++) + { + for(int i=0; i<(int)lambda[(*it_s_vD).first].size(); i++) + { + IVStream +=lambda[(*it_s_vD).first][i] * StreamArray[(*it_s_vD).first][i]; + } + } + + //Check if BuildedFuel is in Model isotopic bounds + (*this).isIVInDomain(IVStream); + + for( it_s_vD = lambda.begin(); it_s_vD != lambda.end(); it_s_vD++) + { + DBGV( "Lambda vector : "<<(*it_s_vD).first ); + + for(int i=0; i < (int)lambda[(*it_s_vD).first].size(); i++) + { + DBGV(lambda[(*it_s_vD).first][i]); + } + } + + DBGL + return lambda; +} + +//________________________________________________________________________ +TTree* EQ_OneParameter::CreateTMVAInputTree(IsotopicVector TheFreshfuel, double ThisTime) +{ + /******Create Input data tree to be interpreted by TMVA::Reader***/ + TTree* InputTree = new TTree(fOutput.c_str(), fOutput.c_str()); + + vector<float> InputTMVA; + for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++) + InputTMVA.push_back(0); + + float Time = 0; + + IsotopicVector IVInputTMVA; + map<ZAI ,string >::iterator it_ZAI_s; + int j = 0; + + for( it_ZAI_s = fMapOfTMVAVariableNames.begin() ; it_ZAI_s != fMapOfTMVAVariableNames.end() ; it_ZAI_s++) + { + InputTree->Branch( ((*it_ZAI_s).second).c_str(), &InputTMVA[j], ((*it_ZAI_s).second + "/F").c_str()); + IVInputTMVA+= ((*it_ZAI_s).first)*1; + j++; + } + + if(ThisTime != -1) + InputTree->Branch("Time" ,&Time ,"Time/F"); + + IsotopicVector IVAccordingToUserInfoFile = TheFreshfuel.GetThisComposition(IVInputTMVA); + double Ntot = IVAccordingToUserInfoFile.GetSumOfAll(); + IVAccordingToUserInfoFile = IVAccordingToUserInfoFile/Ntot; + + j = 0; + + for( it_ZAI_s = fMapOfTMVAVariableNames.begin() ; it_ZAI_s != fMapOfTMVAVariableNames.end() ; it_ZAI_s++) + { + InputTMVA[j] = IVAccordingToUserInfoFile.GetZAIIsotopicQuantity( (*it_ZAI_s).first ) ; + j++; + } + + Time = ThisTime; + InputTree->Fill(); + + return InputTree; + +} +//________________________________________________________________________ +void EQ_OneParameter::CheckTargetParameterConsistency(map < int , string > StreamListPriority, map < string , double > TargetParameterMin, map < string , double > TargetParameterMax) +{ + map < int , string >::iterator it_i_s; + + //Loop on priority order to check if target parameter increases monotously with the material mass + for( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++) + { + double TargetParameterUp = -1.0; //to be sure BUMin is > to BUmax even if BUmin is zero + double TargetParameterDown = 0.0; + + if(TargetParameterMin.find((*it_i_s).second) == TargetParameterMin.end()) + { + break; //if material is not in map, break the loop + } + TargetParameterDown = TargetParameterMin[(*it_i_s).second]; + + if (TargetParameterDown < 0.0 ) + { + ERROR<< "Target parameter evolution should always be positive." <<endl; + ERROR<< "TargetParameterDown = "<< TargetParameterDown<<" is negative "<<endl; + ERROR<< "Check the evolution..." <<endl; + exit(1); + } + TargetParameterUp = TargetParameterMax[(*it_i_s).second]; + + if (TargetParameterDown > TargetParameterUp ) + { + ERROR<< "Target parameter evolution as a function of material mass is not monotonous." <<endl; + ERROR<< "TargetParameterDown = "<< TargetParameterDown<<" is greater than TargetParameterUp = "<< TargetParameterUp<<endl; + ERROR<< "Check the evolution..." <<endl; + exit(1); + } + } +} +//________________________________________________________________________ +double EQ_OneParameter::CalculateTargetParameter(IsotopicVector TheFuel, string TargetParameterName) +{ + double ParameterToCalculate = 0; + if(TargetParameterName=="BurnUpMax") ParameterToCalculate = CalculateBurnUpMax(TheFuel, fModelParameter); + else if(TargetParameterName=="keffBOC") ParameterToCalculate = CalculateKeffAtBOC(TheFuel); + else + { + ERROR<< "Target parameter defined in InformationFile ( "<<TargetParameterName<<" ) doesn't exist" <<endl; + ERROR<< "Possible target parameters for the moment are : BurnUpMax and keffBOC." <<endl; + exit(1); + } + + return ParameterToCalculate ; +} +//________________________________________________________________________ +double EQ_OneParameter::CalculateBurnUpMax(IsotopicVector TheFuel, map<string, double> ModelParameter) +{ + /**************************************************************************/ + //With a dichotomy, the maximal irradiation time (TheFinalTime) is calculated + //When average Kinf is very close (according "Precision") to the threshold + //then the corresponding irradiation time is convert in burnup and returned + /**************************************************************************/ + //Algorithm initialization + double KThreshold = fModelParameter["kThreshold"]; + int NumberOfBatch = (int)fModelParameter["NumberOfBatch"]; + double OldFinalTimeMinus = 0; + double MaximumBU = fMaximalBU; + double MinimumBU = 0 ; + double TheFinalTime = BurnupToSecond((MaximumBU-MinimumBU)/2.); + double OldFinalTimePlus = BurnupToSecond(MaximumBU); + double k_av = 0; //average kinf + double OldPredictedk_av = 0; + + CLASSReader * reader = new CLASSReader( fMapOfTMVAVariableNames ); + reader->AddVariable( "Time" ); + reader->BookMVA( "MLP method" , fTMVAXMLFilePath ); + + for(int b = 0;b<NumberOfBatch;b++) + { + float TheTime = (b+1)*TheFinalTime/NumberOfBatch; + + TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime); + reader->SetInputData( InputTree ); + + OldPredictedk_av += reader->EvaluateRegression( "MLP method" )[0]; + + delete InputTree; + } + OldPredictedk_av /= NumberOfBatch; + + //Algorithm control + int count = 0; + int MaximumLoopCount = 500; + do + { + if(count > MaximumLoopCount ) + { + ERROR << "CRITICAL ! Can't manage to predict burnup\nHint : Try to increase the precision on k effective using :\n YourEQM_MLP_Kinf->SetPCMPrecision(pcm); with pcm the precision in pcm (default 10) REDUCE IT\n If this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr " << endl; + exit(1); + } + + if( (OldPredictedk_av-KThreshold) > 0) //The burnup can be increased + { + OldFinalTimeMinus = TheFinalTime; + TheFinalTime = TheFinalTime + fabs(OldFinalTimePlus - TheFinalTime)/2.; + + if(SecondToBurnup(TheFinalTime) >= (MaximumBU-MaximumBU*GetTargetParameterStDev() ) ) + { delete reader; return MaximumBU; } + } + + else if( (OldPredictedk_av-KThreshold) < 0)//The burnup is too high + { + OldFinalTimePlus = TheFinalTime; + TheFinalTime = TheFinalTime - fabs(OldFinalTimeMinus-TheFinalTime)/2.; + if( SecondToBurnup(TheFinalTime) < (MaximumBU-MinimumBU)/2.*GetTargetParameterStDev() ) + { delete reader; return 0; } + } + + k_av = 0; + for(int b = 0;b<NumberOfBatch;b++) + { + float TheTime = (b+1)*TheFinalTime/NumberOfBatch; + TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime); + reader->SetInputData( InputTree ); + + k_av += reader->EvaluateRegression("MLP method")[0]; + + delete InputTree; + } + k_av/= NumberOfBatch; + //cout<<SecondToBurnup(TheFinalTime)<<" "; + OldPredictedk_av = k_av; + count++; +//std::clog << "-> " << k_av << "\t\t(" << count << ") \t [" << TheFinalTime << "]" << "\t" << OldPredictedk_av-KThreshold << "\t" << GetPCMPrecision() << std::endl; + } while( fabs(OldPredictedk_av-KThreshold) > GetPCMPrecision() ) ; + + delete reader; + //cout<<endl; + return SecondToBurnup(TheFinalTime); +} + +//________________________________________________________________________ +double EQ_OneParameter::CalculateKeffAtBOC(IsotopicVector FreshFuel) +{ + CLASSReader * reader = new CLASSReader( fMapOfTMVAVariableNames ); + reader->BookMVA( "MLP method" , fTMVAXMLFilePath ); + + TTree* InputTree = CreateTMVAInputTree(FreshFuel,-1) ; + reader->SetInputData( InputTree ); + + double keff = reader->EvaluateRegression( "MLP method" )[0]; + + delete InputTree; + + return keff; +} +//________________________________________________________________________ +void EQ_OneParameter::ReadNFO() +{ + DBGL + ifstream NFO(fTMVANFOFilePath.c_str()); + + if(!NFO) + { + ERROR << "Can't find/open file " << fTMVANFOFilePath << endl; + exit(0); + } + + do + { + string line; + getline(NFO,line); + + EQ_OneParameter::ReadLine(line); + + } while(!NFO.eof()); + + DBGL +} +//________________________________________________________________________ +void EQ_OneParameter::ReadLine(string line) +{ + DBGL + + if (!freaded) + { + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + + map<string, EQOP_MthPtr>::iterator it = fKeyword.find(keyword); + + if(it != fKeyword.end()) + (this->*(it->second))( line ); + + freaded = true; + ReadLine(line); + + } + + freaded = false; + + DBGL +} +//________________________________________________________________________ +void EQ_OneParameter::LoadKeyword() +{ + DBGL + fKeyword.insert( pair<string, EQOP_MthPtr>( "k_zail", & EQ_OneParameter::ReadZAIlimits) ); + fKeyword.insert( pair<string, EQOP_MthPtr>( "k_reactor", & EQ_OneParameter::ReadType) ); + fKeyword.insert( pair<string, EQOP_MthPtr>( "k_fuel", & EQ_OneParameter::ReadType) ); + fKeyword.insert( pair<string, EQOP_MthPtr>( "k_massfractionmin", & EQ_OneParameter::ReadEqMinFraction) ); + fKeyword.insert( pair<string, EQOP_MthPtr>( "k_massfractionmax", & EQ_OneParameter::ReadEqMaxFraction) ); + fKeyword.insert( pair<string, EQOP_MthPtr>( "k_list", & EQ_OneParameter::ReadList) ); + fKeyword.insert( pair<string, EQOP_MthPtr>( "k_specpower", & EQ_OneParameter::ReadSpecificPower) ); + if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_zainame", & EQ_OneParameter::ReadZAIName) ); + fKeyword.insert( pair<string, EQOP_MthPtr>( "k_maxburnup", & EQ_OneParameter::ReadMaxBurnUp) ); + if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_targetparameter", & EQ_OneParameter::ReadTargetParameter) ); + if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_predictortype", & EQ_OneParameter::ReadPredictorType) ); + if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_output", & EQ_OneParameter::ReadOutput) ); + fKeyword.insert( pair<string, EQOP_MthPtr>( "k_buffer", & EQ_OneParameter::ReadBuffer) ); + if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_modelparameter", & EQ_OneParameter::ReadModelParameter) ); + if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQOP_MthPtr>( "k_targetparameterstdev", & EQ_OneParameter::ReadTargetParameterStDev) ); + + DBGL +} +//________________________________________________________________________ +void EQ_OneParameter::ReadType(const string &line) +{ + DBGL + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_fuel" && keyword != "k_reactor" ) // Check the keyword + { + ERROR << " Bad keyword : " << keyword << " Not found !" << endl; + exit(1); + } + if( keyword == "k_fuel" ) + fDBFType = StringLine::NextWord(line, pos, ' '); + else if( keyword == "k_reactor" ) + fDBRType = StringLine::NextWord(line, pos, ' '); + + DBGL +} +//________________________________________________________________________ +void EQ_OneParameter::ReadZAIlimits(const string &line) +{ + DBGL + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_zail" ) // Check the keyword + { + ERROR << " Bad keyword : \"k_zail\" 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()); + + double downLimit = atof(StringLine::NextWord(line, pos, ' ').c_str()); + double upLimit = atof(StringLine::NextWord(line, pos, ' ').c_str()); + + if (upLimit < downLimit) + { + double tmp = upLimit; + upLimit = downLimit; + downLimit = tmp; + } + fZAILimits.insert(pair<ZAI, pair<double, double> >(ZAI(Z,A,I), pair<double,double>(downLimit, upLimit))); + DBGL +} +//________________________________________________________________________ +void EQ_OneParameter::ReadList(const string &line) +{ + DBGL + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_list" ) // Check the keyword + { + ERROR << " Bad keyword : \"k_list\" not found !" << endl; + exit(1); + } + string ListName= StringLine::NextWord(line, pos, ' '); + 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()); + double Q = atof(StringLine::NextWord(line, pos, ' ').c_str()); + fStreamList[ListName].Add(Z, A, I, Q); + + DBGL +} +//________________________________________________________________________ +void EQ_OneParameter::ReadEqMinFraction(const string &line) +{ + DBGL + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_massfractionmin" ) // Check the keyword + { + ERROR << " Bad keyword : \"k_massfractionmin\" not found !" << endl; + exit(1); + } + string ListName= StringLine::NextWord(line, pos, ' '); + double Q = atof(StringLine::NextWord(line, pos, ' ').c_str()); + fStreamListEqMMassFractionMin[ListName] = Q; + + DBGL +} + +//________________________________________________________________________ +void EQ_OneParameter::ReadEqMaxFraction(const string &line) +{ + DBGL + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_massfractionmax" ) // Check the keyword + { + ERROR << " Bad keyword : \"k_massfractionmax\" not found !" << endl; + exit(1); + } + string ListName= StringLine::NextWord(line, pos, ' '); + double Q = atof(StringLine::NextWord(line, pos, ' ').c_str()); + fStreamListEqMMassFractionMax[ListName] = Q; + + DBGL +} + +//________________________________________________________________________ +void EQ_OneParameter::ReadSpecificPower(const string &line) +{ + DBGL + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_specpower") // Check the keyword + { + ERROR << " Bad keyword : \"k_specpower\" Not found !" << endl; + exit(1); + } + + fSpecificPower = atof(StringLine::NextWord(line, pos, ' ').c_str()); + + DBGL +} +//________________________________________________________________________ +void EQ_OneParameter::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 EQ_OneParameter::ReadMaxBurnUp(const string &line) +{ + DBGL + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_maxburnup" ) // Check the keyword + { + ERROR << " Bad keyword : \"k_maxburnup\" not found !" << endl; + exit(1); + } + + fMaximalBU = atof(StringLine::NextWord(line, pos, ' ').c_str()); + + DBGL +} + +//________________________________________________________________________ +void EQ_OneParameter::ReadTargetParameter(const string &line) +{ + DBGL + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_targetparameter" ) // Check the keyword + { + ERROR << " Bad keyword : \"k_targetparameter\" not found !" << endl; + exit(1); + } + + fTargetParameter = StringLine::NextWord(line, pos, ' '); + + DBGL +} + +//________________________________________________________________________ +void EQ_OneParameter::ReadPredictorType(const string &line) +{ + DBGL + + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_predictortype" ) // Check the keyword + { + ERROR << " Bad keyword : \"k_predictortype\" not found !" << endl; + exit(1); + } + + fPredictorType = StringLine::NextWord(line, pos, ' '); + + DBGL +} +//________________________________________________________________________ +void EQ_OneParameter::ReadOutput(const string &line) +{ + DBGL + + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_output" ) // Check the keyword + { + ERROR << " Bad keyword : \"k_output\" not found !" << endl; + exit(1); + } + + fOutput = StringLine::NextWord(line, pos, ' '); + + DBGL +} +//________________________________________________________________________ +void EQ_OneParameter::ReadBuffer(const string &line) +{ + DBGL + + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_buffer" ) // Check the keyword + { + ERROR << " Bad keyword : \"k_buffer\" not found !" << endl; + exit(1); + } + + fBuffer = StringLine::NextWord(line, pos, ' '); + + DBGL +} +//________________________________________________________________________ +void EQ_OneParameter::ReadModelParameter(const string &line) +{ + DBGL + + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_modelparameter" ) // Check the keyword + { + ERROR << " Bad keyword : \"k_modelparameter\" not found !" << endl; + exit(1); + } + + keyword = StringLine::NextWord(line, pos, ' '); + + fModelParameter[keyword] = -1; + + DBGL +} + +//________________________________________________________________________ +void EQ_OneParameter::ReadTargetParameterStDev(const string &line) +{ + DBGL + + int pos = 0; + string keyword = tlc(StringLine::NextWord(line, pos, ' ')); + if( keyword != "k_targetparameterstdev" ) // Check the keyword + { + ERROR << " Bad keyword : \"k_targetparameterstdev\" not found !" << endl; + exit(1); + } + + fTargetParameterStDev = atof(StringLine::NextWord(line, pos, ' ').c_str());; + + DBGL +} + +//________________________________________________________________________ +void EQ_OneParameter::PrintInfo() +{ + INFO << "Reactor Type : "<< fDBRType << endl; + INFO << "Fuel Type : "<< fDBFType << endl; + INFO << "Specific Power [W/g]: "<< fSpecificPower << endl; + + map < string , IsotopicVector >::iterator it_s_IV; + map < string , double >::iterator it_s_D; + + for( it_s_IV = fStreamList.begin(); it_s_IV != fStreamList.end(); it_s_IV++) + { + INFO <<(* it_s_IV).first<<" (Z A I) :" << endl; + map<ZAI ,double >::iterator it1; + map<ZAI ,double > fMap1 = fStreamList[(* it_s_IV).first].GetIsotopicQuantity(); + for(it1 = fMap1.begin() ; it1 != fMap1.end() ; it1++) + INFO << (*it1).first.Z() <<" "<< (*it1).first.A() <<" "<< (*it1).first.I() << endl; + } + INFO<<"Minimum fraction in the fuel for each material : "<<endl; + for( it_s_D = fStreamListEqMMassFractionMin.begin(); it_s_D != fStreamListEqMMassFractionMin.end(); it_s_D++) + { + INFO <<(* it_s_D).first<<" "<<fStreamListEqMMassFractionMin[(* it_s_D).first]<<endl; + } + + INFO<<"Maximum fraction in the fuel for each material : "<<endl; + for( it_s_D = fStreamListEqMMassFractionMax.begin(); it_s_D != fStreamListEqMMassFractionMax.end(); it_s_D++) + { + INFO <<(* it_s_D).first<<" "<<fStreamListEqMMassFractionMax[(* it_s_D).first]<<endl; + } + + + INFO<<"ZAI limits (validity domain)[prop in fresh fuel] (Z A I min max) :"<<endl; + for (map< ZAI,pair<double,double> >::iterator Domain_it = fZAILimits.begin(); Domain_it != fZAILimits.end(); Domain_it++) + { + double ThatZAIMin = Domain_it->second.first; + double ThatZAIMax = Domain_it->second.second; + int Z = Domain_it->first.Z(); + int A = Domain_it->first.A(); + int I = Domain_it->first.I(); + + INFO <<ThatZAIMin<<" < ZAI ("<< Z<< " " << A << " " << I<<")"<<" < "<<ThatZAIMax<< endl; + } + +} diff --git a/source/Model/Equivalence/EQ_OneParameter.hxx b/source/Model/Equivalence/EQ_OneParameter.hxx new file mode 100644 index 0000000000000000000000000000000000000000..6922e56bb8baff759a96715afce03fef4261960c --- /dev/null +++ b/source/Model/Equivalence/EQ_OneParameter.hxx @@ -0,0 +1,285 @@ +#ifndef _EQONEPARAMETER_HXX_ +#define _EQONEPARAMETER_HXX_ + +/*! + \file + \brief Header file for EQ_OneParameter class. + + + @author BLG + @author BaM + @author FaC + @version 3.0 + */ + +#include "IsotopicVector.hxx" +#include <math.h> +#include "TTree.h" +#include <map> +#include "CLASSObject.hxx" + + +using namespace std; + +class EQ_OneParameter; +#ifndef __ROOTCLING__ +typedef void (EQ_OneParameter::*EQOP_MthPtr)( const string & ); +#endif + +//-----------------------------------------------------------------------------// + +//! Determines how to build a fresh fuel +/*! + Define an EQ_OneParameter. + The aim of these class is to gather all the commum properties of all + Equivalence Model. + + \warning + Never instantiate EQ_OneParameter in your CLASS input but it's derivated class + @see EQM_FBR_BakerRoss_MOX. + @see EQM_PWR_MLP_MOX + @see EQM_FBR_MLP_Keff.hxx + @see EQM_PWR_MLP_MOX_Am + @see EQM_FBR_MLP_Keff_BOUND + @see EQM_PWR_POL_UO2 + @see EQM_MLP_Kinf.hxx + @see EQM_PWR_QUAD_MOX + @see EQM_PWR_LIN_MOX + + @author BLG + @author BaM + @author FaC + + @version 3.0 + */ +//________________________________________________________________________ + + +class EQ_OneParameter : public EquivalenceModel +{ + public : + /*! + \name Constructor/Desctructor + */ + //@{ + EQ_OneParameter(string TMVAXMLFilePath, string TMVANFOFilePath); //!< Default constructor with path + EQ_OneParameter(CLASSLogger* log, string TMVAXMLFilePath, string TMVANFOFilePath); //!< Logger constructor with path + + EQ_OneParameter(string TMVANFOFilePath); //!< Default constructor without Eq Model + EQ_OneParameter(CLASSLogger* log, string TMVANFOFilePath); //!< Logger constructor Without Eq Model + + virtual ~EQ_OneParameter(); //!< Destructor + //@} + + /*! + \name Fuel Construction Method + */ + //@{ + + //{ + /// BuildFuel function. + /*! + Build the fuel following the equivalance model with the proper requierment in term of mass, burnup.... + \param double burnup reached by the fuel at the end of irradiation + \param double HMMass, Heavy metal mass needed + \param map < string , vector <IsotopicVector> > StreamArray, the string is the stream code (fissile fertile ,...) the IsotopicVector the fraction of each IV to take in the (fissile, fertile,..) stock . + */ + virtual map <string , vector<double> > BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray, map < string , double> StreamListMassFractionMin, map < string , double> StreamListMassFractionMax, map < int , string> StreamListPriority, map < string , bool> StreamListIsBuffer); + //} + + + //@} + TTree* CreateTMVAInputTree(IsotopicVector TheFreshfuel, double ThisTime) ; //!<Create input tmva tree to be read by ExecuteTMVA + double CalculateTargetParameter(IsotopicVector TheFuel, string TargetParameterName); //!<Get a fuel parameter associated to the fuel ---> ex : BurnUpMax, keffBOC, keffEOC, ... + double CalculateBurnUpMax(IsotopicVector TheFuel, map<string, double> ModelParameter);//!<Calculate the BU max associated to a fuel composition based on MLP prediction (suitable for PWR) + double CalculateKeffAtBOC(IsotopicVector TheFuel); //!<Calculate the keff at BOC associated to a fuel composition based on MLP prediction (suitable for SFR) + + string GetTMVAXMLFilePath() {return fTMVAXMLFilePath;} // Return the path to TMVA XML File path + string GetTMVANFOFilePath() {return fTMVANFOFilePath;} // Return the path to TMVA NFO File path + void SetTMVAXMLFilePath(string TMVAXMLFilePath) {fTMVAXMLFilePath = TMVAXMLFilePath;} // Set the path to TMVA XML File path + void SetTMVANFOFilePath(string TMVANFOFilePath) {fTMVANFOFilePath = TMVANFOFilePath;} // Set the path to TMVA NFO File path + + /*! + \name Get/Set Method + */ + //@{ + + int GetStreamListNumber(){return fStreamList.size();}; + int GetMaxIterration() const { return fMaxIterration; } //!< Max iterration in build fueld algorythm + double GetTargetParameterStDev(){return fTargetParameterStDev;}//!< Get the precision on fTargetParameterStDev + double GetStreamListEqMMassFractionMax(string keyword){return fStreamListEqMMassFractionMax[keyword] ;} + double GetStreamListEqMMassFractionMin(string keyword){return fStreamListEqMMassFractionMin[keyword] ;} + double GetPCMPrecision(){return fPCMprecision/1e5;}//!< Get the precision on @f$\langle k \rangle@f$ prediction []. Neural network predictor constructors + + void SetModelParameter(string sMP, double dMP) { fModelParameter[sMP] = dMP; } //!< Set Model Parameters precised in NFO file + map<string, double> GetModelParameter() { return fModelParameter; } //!< Get Model Parameters precised in NFO file + + void SetMaxIterration(int val) { fMaxIterration = val; } //!< Max iteration in build fuel algorithm + void SetTargetParameterStDev(double TPSD){fTargetParameterStDev = TPSD;} //!< Set the precision on Target Parameter + void SetStreamListEqMMassFractionMax(string keyword, double value){fStreamListEqMMassFractionMax[keyword] = value;} + void SetStreamListEqMMassFractionMin(string keyword, double value){fStreamListEqMMassFractionMin[keyword] = value;} + + void SetPCMPrecision(double pcm){fPCMprecision = pcm;} //!< Set the precision on @f$\langle k \rangle@f$ prediction [pcm]. Neural network predictor constructors + + /*! + \name Time <-> Burnup conversion + */ + //@{ + + double SecondToBurnup(double Second){return Second*fSpecificPower/(24*3.6e6);} + double BurnupToSecond(double BurnUp){return BurnUp/fSpecificPower*(24*3.6e6);} + + //@} + + //@} + + /*! + \name Reading NFO related Method + */ + //@{ + void ReadNFO(); + virtual void ReadLine(string line); + + + virtual void LoadKeyword(); + void ReadZAIlimits(const string &line); + void ReadType(const string &line); + + //{ + /// ReadZAIName : read the zai name in the TMWA MLP model + /*! + \param line : line suppossed to contain the ZAI name starts with "k_zainame" keyword + */ + void ReadZAIName(const string &line); + //} + + //{ + /// ReadMaxBurnUp : read a guessed (very overestimated) maximum burnup a fuel can reach (purpose : algorithm initialization) + /*! + \param line : line suppossed to contain the ZAI name starts with "k_maxburnup" keyword + */ + void ReadMaxBurnUp(const string &line); + //} + + //{ + /// ReadSpecificPower : read the Specific Power of the DataBase + /*! + \param line : line suppossed to contain the Specific Power information starts with "k_specpower" keyword + */ + void ReadSpecificPower(const string &line); + //} + + //{ + /// ReadTargetParameter : type of target parameter optimized in build fuel (ex. BUmax) + /*! + \param line : line suppossed to contain the Target Parameter information starts with "k_targetparameter" keyword + */ + void ReadTargetParameter(const string &line); + //} + + //{ + /// ReadOutput : read the output type of the predictor (ex : kinf) + /*! + \param line : line suppossed to contain the Specific Power information starts with "k_output" keyword + */ + void ReadOutput(const string &line); + //} + + //{ + /// ReadBuffer : read the Buffer material name in the fuel + /*! + \param line : line suppossed to contain the Buffer information starts with "k_buffer" keyword + */ + void ReadBuffer(const string &line); + //} + + //{ + /// ReadModelParameter : read the name of equivalence model parameter + /*! + \param line : line suppossed to contain the Buffer information starts with "k_modelparameter" keyword + */ + void ReadModelParameter(const string &line); + //} + + //{ + /// ReadPredictorType: read the type of predictor used (ex : MLP) + /*! + \param line : line suppossed to contain the Buffer information starts with "k_predictortype" keyword + */ + void ReadPredictorType(const string &line); + //} + + //{ + /// ReadTargetParameterStDev: read the target parameter standard deviation + /*! + \param line : line suppossed to contain the Buffer information starts with "k_targetparameterstdev" keyword + */ + void ReadTargetParameterStDev(const string &line); + //} + + void PrintInfo(); //Print the information red in the INFO stream + + //{ + /// ReadFissil : read the zai name in the TMWA MLP model starts with "k_fissil" keyword + /*! + \param line : line suppossed to contain the fissil list + */ + void ReadList(const string &line); + + void ReadEqMaxFraction(const string &line); + void ReadEqMinFraction(const string &line); + + IsotopicVector BuildFuelToTest(map < string, vector<double> >& lambda, map < string , vector <IsotopicVector> > const& StreamArray, double HMMass, map <string, bool> StreamListIsBuffer); //Build a fuel with the buffer according to fissile lambda + void CheckTargetParameterConsistency(map < int , string > StreamListPriority, map < string , double > TargetParameterMin, map < string , double > TargetParameterMax); + + protected : + + bool fUseTMVAPredictor; //!< Bool that says if we need a TMVA predictor. If not, fuel fraction isimpoased by the FP. + + map < string , double> fStreamListEqMMassFractionMax; //!< Map that contains lists of stream according to the EqModel with mass maximum fraction + map < string , double> fStreamListEqMMassFractionMin; //!< Map that contains lists of stream according to the EqModel with mass minimum fraction + + string fPredictorType ; //!< Type of predictor used in Equivalence Model (ex: MLP) + string fOutput ; //!< Type of output calculated by the predictor + string fBuffer ; //!< Name of material used as buffer in fuel + + map<string, double> fModelParameter ; ///< Map of equivalence model parameter + + map<ZAI, string> fMapOfTMVAVariableNames; //!< List of TMVA input variable names (read from fMLPInformationFile ) , name depends on the training step + + double fTargetParameterStDev; //!< Precision on target parameter calculation + + double fMaximalBU; //!< The Maximum burn-up of the model in GWd/t + string fTargetParameter; //!< Type of target parameter optimized in build fuel (ex. BUmax) + int fMaxIterration; //!< Max iterrations in build fueld algorithm + + string fTMVAXMLFilePath; //!<The weight needed by TMVA to construct and execute the multilayer perceptron + string fTMVANFOFilePath; //!<The weight needed by TMVA to construct and execute the multilayer perceptron + + +#ifndef __ROOTCLING__ + map<string, EQOP_MthPtr> fKeyword; +#endif + + map< ZAI, pair<double,double> > fZAILimits; //!< Fresh fuel range : map<ZAI<min edge ,max edge >> + + string fInformationFile; //!< file containing Reactor Type, Fuel type, HM mass, Power, time vector, and TMVA input variables names (looks the manual for format details) + string fDBFType; //!< Fuel Type (e.g MOX, UOX, ThU, ThPu ...) + string fDBRType; //!< Reactor Type (e.g PWR, FBR-Na, ADS..) + + private : + + double fPCMprecision; //!< precision on @f$\langle k \rangle@f$ prediction [pcm] + +}; + +#endif + + + + + + + + + diff --git a/source/Model/Equivalence/TMP/EQM_ADS_MLP_FixedRatioPuAM.cxx b/source/Model/Equivalence/OLD/EQM_ADS_MLP_FixedRatioPuAM.cxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_ADS_MLP_FixedRatioPuAM.cxx rename to source/Model/Equivalence/OLD/EQM_ADS_MLP_FixedRatioPuAM.cxx diff --git a/source/Model/Equivalence/TMP/EQM_ADS_MLP_FixedRatioPuAM.hxx b/source/Model/Equivalence/OLD/EQM_ADS_MLP_FixedRatioPuAM.hxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_ADS_MLP_FixedRatioPuAM.hxx rename to source/Model/Equivalence/OLD/EQM_ADS_MLP_FixedRatioPuAM.hxx diff --git a/source/Model/Equivalence/TMP/EQM_ADS_MLP_RatioPuAM.cxx b/source/Model/Equivalence/OLD/EQM_ADS_MLP_RatioPuAM.cxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_ADS_MLP_RatioPuAM.cxx rename to source/Model/Equivalence/OLD/EQM_ADS_MLP_RatioPuAM.cxx diff --git a/source/Model/Equivalence/TMP/EQM_ADS_MLP_RatioPuAM.hxx b/source/Model/Equivalence/OLD/EQM_ADS_MLP_RatioPuAM.hxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_ADS_MLP_RatioPuAM.hxx rename to source/Model/Equivalence/OLD/EQM_ADS_MLP_RatioPuAM.hxx diff --git a/source/Model/Equivalence/TMP/EQM_FBR_BakerRoss_MOX.cxx b/source/Model/Equivalence/OLD/EQM_FBR_BakerRoss_MOX.cxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_FBR_BakerRoss_MOX.cxx rename to source/Model/Equivalence/OLD/EQM_FBR_BakerRoss_MOX.cxx diff --git a/source/Model/Equivalence/TMP/EQM_FBR_BakerRoss_MOX.hxx b/source/Model/Equivalence/OLD/EQM_FBR_BakerRoss_MOX.hxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_FBR_BakerRoss_MOX.hxx rename to source/Model/Equivalence/OLD/EQM_FBR_BakerRoss_MOX.hxx diff --git a/source/Model/Equivalence/TMP/EQM_FBR_MLP_Keff.cxx b/source/Model/Equivalence/OLD/EQM_FBR_MLP_Keff.cxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_FBR_MLP_Keff.cxx rename to source/Model/Equivalence/OLD/EQM_FBR_MLP_Keff.cxx diff --git a/source/Model/Equivalence/TMP/EQM_FBR_MLP_Keff.hxx b/source/Model/Equivalence/OLD/EQM_FBR_MLP_Keff.hxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_FBR_MLP_Keff.hxx rename to source/Model/Equivalence/OLD/EQM_FBR_MLP_Keff.hxx diff --git a/source/Model/Equivalence/TMP/EQM_MLP_Kinf.cxx b/source/Model/Equivalence/OLD/EQM_MLP_Kinf.cxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_MLP_Kinf.cxx rename to source/Model/Equivalence/OLD/EQM_MLP_Kinf.cxx diff --git a/source/Model/Equivalence/TMP/EQM_MLP_Kinf.hxx b/source/Model/Equivalence/OLD/EQM_MLP_Kinf.hxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_MLP_Kinf.hxx rename to source/Model/Equivalence/OLD/EQM_MLP_Kinf.hxx diff --git a/source/Model/Equivalence/TMP/EQM_MLP_PWR_MOxEUS.cxx b/source/Model/Equivalence/OLD/EQM_MLP_PWR_MOxEUS.cxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_MLP_PWR_MOxEUS.cxx rename to source/Model/Equivalence/OLD/EQM_MLP_PWR_MOxEUS.cxx diff --git a/source/Model/Equivalence/TMP/EQM_MLP_PWR_MOxEUS.hxx b/source/Model/Equivalence/OLD/EQM_MLP_PWR_MOxEUS.hxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_MLP_PWR_MOxEUS.hxx rename to source/Model/Equivalence/OLD/EQM_MLP_PWR_MOxEUS.hxx diff --git a/source/Model/Equivalence/TMP/EQM_PWR_FixedContent.cxx b/source/Model/Equivalence/OLD/EQM_PWR_FixedContent.cxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_PWR_FixedContent.cxx rename to source/Model/Equivalence/OLD/EQM_PWR_FixedContent.cxx diff --git a/source/Model/Equivalence/TMP/EQM_PWR_FixedContent.hxx b/source/Model/Equivalence/OLD/EQM_PWR_FixedContent.hxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_PWR_FixedContent.hxx rename to source/Model/Equivalence/OLD/EQM_PWR_FixedContent.hxx diff --git a/source/Model/Equivalence/TMP/EQM_PWR_MLP_MOX.cxx b/source/Model/Equivalence/OLD/EQM_PWR_MLP_MOX.cxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_PWR_MLP_MOX.cxx rename to source/Model/Equivalence/OLD/EQM_PWR_MLP_MOX.cxx diff --git a/source/Model/Equivalence/TMP/EQM_PWR_MLP_MOX.hxx b/source/Model/Equivalence/OLD/EQM_PWR_MLP_MOX.hxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_PWR_MLP_MOX.hxx rename to source/Model/Equivalence/OLD/EQM_PWR_MLP_MOX.hxx diff --git a/source/Model/Equivalence/TMP/EQM_PWR_MLP_MOX_Am.cxx b/source/Model/Equivalence/OLD/EQM_PWR_MLP_MOX_Am.cxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_PWR_MLP_MOX_Am.cxx rename to source/Model/Equivalence/OLD/EQM_PWR_MLP_MOX_Am.cxx diff --git a/source/Model/Equivalence/TMP/EQM_PWR_MLP_MOX_Am.hxx b/source/Model/Equivalence/OLD/EQM_PWR_MLP_MOX_Am.hxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_PWR_MLP_MOX_Am.hxx rename to source/Model/Equivalence/OLD/EQM_PWR_MLP_MOX_Am.hxx diff --git a/source/Model/Equivalence/TMP/EQM_PWR_POL_UO2.cxx b/source/Model/Equivalence/OLD/EQM_PWR_POL_UO2.cxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_PWR_POL_UO2.cxx rename to source/Model/Equivalence/OLD/EQM_PWR_POL_UO2.cxx diff --git a/source/Model/Equivalence/TMP/EQM_PWR_POL_UO2.hxx b/source/Model/Equivalence/OLD/EQM_PWR_POL_UO2.hxx similarity index 100% rename from source/Model/Equivalence/TMP/EQM_PWR_POL_UO2.hxx rename to source/Model/Equivalence/OLD/EQM_PWR_POL_UO2.hxx diff --git a/source/Model/XS/CMakeLists.txt b/source/Model/XS/CMakeLists.txt index 08b512d34e9d11d2b6f86c28a7fdfa6b9947a30f..080ea25c8ef610375a6ad57b71356b326c5f4338 100644 --- a/source/Model/XS/CMakeLists.txt +++ b/source/Model/XS/CMakeLists.txt @@ -10,5 +10,4 @@ LIST( APPEND CLASS_SRC_FILE_XS "${CMAKE_CURRENT_SOURCE_DIR}/XSM_MLP.cxx") SET( CLASS_SRC_FILE_XS ${CLASS_SRC_FILE_XS} PARENT_SCOPE) -message("------- CLASS_SRC_FILE_XS VARIABLE -------") MESSAGE( STATUS "CLASS_SRC_FILE_XS=" ${CLASS_SRC_FILE_XS}) diff --git a/source/Model/XS/XSM_MLP.cxx b/source/Model/XS/XSM_MLP.cxx index 8000af6bbfe77965f7da95419918ca4eed499c07..5a347d687e27f66a08fbdb2fa85d5c8028908243 100644 --- a/source/Model/XS/XSM_MLP.cxx +++ b/source/Model/XS/XSM_MLP.cxx @@ -36,7 +36,6 @@ //________________________________________________________________________ XSM_MLP::XSM_MLP(string TMVA_Weight_Directory,string InformationFile, bool IsTimeStep):XSModel(new CLASSLogger("XSM_MLP.log")) { - fIsStepTime = IsTimeStep; fTMVAWeightFolder = TMVA_Weight_Directory; @@ -52,11 +51,9 @@ XSM_MLP::XSM_MLP(string TMVA_Weight_Directory,string InformationFile, bool IsTim ReadNFO(); } - //________________________________________________________________________ XSM_MLP::XSM_MLP(CLASSLogger* Log,string TMVA_Weight_Directory,string InformationFile, bool IsTimeStep):XSModel(Log) { - fIsStepTime = IsTimeStep; fTMVAWeightFolder = TMVA_Weight_Directory; @@ -70,9 +67,7 @@ XSM_MLP::XSM_MLP(CLASSLogger* Log,string TMVA_Weight_Directory,string Informatio LoadKeyword(); ReadNFO(); - } - //________________________________________________________________________ XSM_MLP::~XSM_MLP() { diff --git a/source/Model/XS/XSM_MLP.hxx b/source/Model/XS/XSM_MLP.hxx index f7c748f12b35093a111c89df6551dae09a0ecdfd..8ade4df29b8c519d41fe22fcc94d1116079b69e4 100644 --- a/source/Model/XS/XSM_MLP.hxx +++ b/source/Model/XS/XSM_MLP.hxx @@ -1,8 +1,6 @@ - #ifndef _XSM_MLP_HXX #define _XSM_MLP_HXX - /*! \file \brief Header file for XSM_MLP class. @@ -26,7 +24,7 @@ using namespace std; class XSM_MLP; -#ifndef __CINT__ +#ifndef __ROOTCLING__ typedef void (XSM_MLP::*XS_MLP_DMthPtr)( const string & ) ; #endif //-----------------------------------------------------------------------------// @@ -142,7 +140,7 @@ class XSM_MLP : public XSModel map<ZAI,string> fMapOfTMVAVariableNames;//!< List of TMVA input variable names (read from fMLPInformationFile ) , name depends on the training step -#ifndef __CINT__ +#ifndef __ROOTCLING__ map<string, XS_MLP_DMthPtr> fDKeyword; #endif }; diff --git a/source/include/EquivalenceModel.hxx b/source/include/EquivalenceModel.hxx index 3064bf23b5ca87d39ddb265bb2de887573fcd153..b9ec82cfd682562d1e87b8d5806f4c44a77abf2a 100644 --- a/source/include/EquivalenceModel.hxx +++ b/source/include/EquivalenceModel.hxx @@ -5,8 +5,8 @@ /*! \file \brief Header file for EquivalenceModel class. - - + + @author BLG @author BaM @author FaC @@ -29,28 +29,28 @@ typedef void (EquivalenceModel::*EQM_MthPtr)( const string & ); //-----------------------------------------------------------------------------// -//! Determines how to build a fresh fuel +//! Determines how to build a fresh fuel /*! Define an EquivalenceModel. The aim of these class is to gather all the commum properties of all Equivalence Model. - + \warning Never instantiate EquivalenceModel in your CLASS input but it's derivated class @see EQM_FBR_BakerRoss_MOX. @see EQM_PWR_MLP_MOX - @see EQM_FBR_MLP_Keff.hxx + @see EQM_FBR_MLP_Keff.hxx @see EQM_PWR_MLP_MOX_Am @see EQM_FBR_MLP_Keff_BOUND @see EQM_PWR_POL_UO2 - @see EQM_MLP_Kinf.hxx + @see EQM_MLP_Kinf.hxx @see EQM_PWR_QUAD_MOX @see EQM_PWR_LIN_MOX @author BLG @author BaM @author FaC - + @version 3.0 */ //________________________________________________________________________ @@ -58,244 +58,54 @@ typedef void (EquivalenceModel::*EQM_MthPtr)( const string & ); class EquivalenceModel : public CLASSObject { - public : +public : /*! \name Constructor/Desctructor */ //@{ - EquivalenceModel(string TMVAXMLFilePath, string TMVANFOFilePath); //!< Default constructor with path - EquivalenceModel(CLASSLogger* log, string TMVAXMLFilePath, string TMVANFOFilePath); //!< Logger constructor with path - - EquivalenceModel(string TMVANFOFilePath); //!< Default constructor without Eq Model - EquivalenceModel(CLASSLogger* log, string TMVANFOFilePath); //!< Logger constructor Without Eq Model + EquivalenceModel(); //!< Default constructor with path + EquivalenceModel(CLASSLogger* log); //!< Logger constructor with path virtual ~EquivalenceModel(); //!< Destructor //@} - - /*! - \name Fuel Construction Method - */ - //@{ - - //{ - /// BuildFuel function. - /*! - Build the fuel following the equivalance model with the proper requierment in term of mass, burnup.... - \param double burnup reached by the fuel at the end of irradiation - \param double HMMass, Heavy metal mass needed - \param map < string , vector <IsotopicVector> > StreamArray, the string is the stream code (fissile fertile ,...) the IsotopicVector the fraction of each IV to take in the (fissile, fertile,..) stock . - */ - virtual map <string , vector<double> > BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray, map < string , double> StreamListMassFractionMin, map < string , double> StreamListMassFractionMax, map < int , string> StreamListPriority, map < string , bool> StreamListIsBuffer); - //} - - - //@} - TTree* CreateTMVAInputTree(IsotopicVector TheFreshfuel, double ThisTime) ; //!<Create input tmva tree to be read by ExecuteTMVA - double CalculateTargetParameter(IsotopicVector TheFuel, string TargetParameterName); //!<Get a fuel parameter associated to the fuel ---> ex : BurnUpMax, keffBOC, keffEOC, ... - double CalculateBurnUpMax(IsotopicVector TheFuel, map<string, double> ModelParameter);//!<Calculate the BU max associated to a fuel composition based on MLP prediction (suitable for PWR) - double CalculateKeffAtBOC(IsotopicVector TheFuel); //!<Calculate the keff at BOC associated to a fuel composition based on MLP prediction (suitable for SFR) - - string GetTMVAXMLFilePath() {return fTMVAXMLFilePath;} // Return the path to TMVA XML File path - string GetTMVANFOFilePath() {return fTMVANFOFilePath;} // Return the path to TMVA NFO File path - void SetTMVAXMLFilePath(string TMVAXMLFilePath) {fTMVAXMLFilePath = TMVAXMLFilePath;} // Set the path to TMVA XML File path - void SetTMVANFOFilePath(string TMVANFOFilePath) {fTMVANFOFilePath = TMVANFOFilePath;} // Set the path to TMVA NFO File path - - /*! - \name Get/Set Method - */ - //@{ - - IsotopicVector GetStreamList(string keyword) {return fStreamList[keyword];} //!<return the list of ZAI of stream type keyword - map < string, IsotopicVector> GetAllStreamList() {return fStreamList;} //!<return all the lists - - int GetStreamListNumber(){return fStreamList.size();}; - int GetMaxIterration() const { return fMaxIterration; } //!< Max iterration in build fueld algorythm - double GetTargetParameterStDev(){return fTargetParameterStDev;}//!< Get the precision on fTargetParameterStDev - double GetStreamListEqMMassFractionMax(string keyword){return fStreamListEqMMassFractionMax[keyword] ;} - double GetStreamListEqMMassFractionMin(string keyword){return fStreamListEqMMassFractionMin[keyword] ;} - double GetPCMPrecision(){return fPCMprecision/1e5;}//!< Get the precision on @f$\langle k \rangle@f$ prediction []. Neural network predictor constructors - - void SetModelParameter(string sMP, double dMP) { fModelParameter[sMP] = dMP; } //!< Set Model Parameters precised in NFO file - map<string, double> GetModelParameter() { return fModelParameter; } //!< Get Model Parameters precised in NFO file - - void SetMaxIterration(int val) { fMaxIterration = val; } //!< Max iteration in build fuel algorithm - void SetTargetParameterStDev(double TPSD){fTargetParameterStDev = TPSD;} //!< Set the precision on Target Parameter - void SetStreamListEqMMassFractionMax(string keyword, double value){fStreamListEqMMassFractionMax[keyword] = value;} - void SetStreamListEqMMassFractionMin(string keyword, double value){fStreamListEqMMassFractionMin[keyword] = value;} - - void SetPCMPrecision(double pcm){fPCMprecision = pcm;} //!< Set the precision on @f$\langle k \rangle@f$ prediction [pcm]. Neural network predictor constructors - - /*! - \name Time <-> Burnup conversion - */ - //@{ - - double SecondToBurnup(double Second){return Second*fSpecificPower/(24*3.6e6);} - double BurnupToSecond(double BurnUp){return BurnUp/fSpecificPower*(24*3.6e6);} - - //@} - - //@} - - /*! - \name Reading NFO related Method - */ - //@{ - void ReadNFO(); - virtual void ReadLine(string line); - - - virtual void LoadKeyword(); - void ReadZAIlimits(const string &line); - void ReadType(const string &line); - - //{ - /// ReadZAIName : read the zai name in the TMWA MLP model - /*! - \param line : line suppossed to contain the ZAI name starts with "k_zainame" keyword - */ - void ReadZAIName(const string &line); - //} - - //{ - /// ReadMaxBurnUp : read a guessed (very overestimated) maximum burnup a fuel can reach (purpose : algorithm initialization) - /*! - \param line : line suppossed to contain the ZAI name starts with "k_maxburnup" keyword - */ - void ReadMaxBurnUp(const string &line); - //} - - //{ - /// ReadSpecificPower : read the Specific Power of the DataBase - /*! - \param line : line suppossed to contain the Specific Power information starts with "k_specpower" keyword - */ - void ReadSpecificPower(const string &line); - //} - //{ - /// ReadTargetParameter : type of target parameter optimized in build fuel (ex. BUmax) - /*! - \param line : line suppossed to contain the Target Parameter information starts with "k_targetparameter" keyword - */ - void ReadTargetParameter(const string &line); - //} + map < string, IsotopicVector> GetAllStreamList() {return fStreamList;} // return all the lists - //{ - /// ReadOutput : read the output type of the predictor (ex : kinf) - /*! - \param line : line suppossed to contain the Specific Power information starts with "k_output" keyword - */ - void ReadOutput(const string &line); - //} + virtual map <string , vector<double> > BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray, map < string , double> StreamListMassFractionMin, map < string , double> StreamListMassFractionMax, map < int , string> StreamListPriority, map < string , bool> StreamListIsBuffer); + + double SecondToBurnup(double Second) {return Second * fSpecificPower / (24 * 3.6e6);} + double BurnupToSecond(double BurnUp) {return BurnUp / fSpecificPower * (24 * 3.6e6);} - //{ - /// ReadBuffer : read the Buffer material name in the fuel - /*! - \param line : line suppossed to contain the Buffer information starts with "k_buffer" keyword - */ - void ReadBuffer(const string &line); - //} - - //{ - /// ReadModelParameter : read the name of equivalence model parameter - /*! - \param line : line suppossed to contain the Buffer information starts with "k_modelparameter" keyword - */ - void ReadModelParameter(const string &line); - //} - - //{ - /// ReadPredictorType: read the type of predictor used (ex : MLP) - /*! - \param line : line suppossed to contain the Buffer information starts with "k_predictortype" keyword - */ - void ReadPredictorType(const string &line); - //} - - //{ - /// ReadTargetParameterStDev: read the target parameter standard deviation - /*! - \param line : line suppossed to contain the Buffer information starts with "k_targetparameterstdev" keyword - */ - void ReadTargetParameterStDev(const string &line); - //} - - void PrintInfo(); //Print the information red in the INFO stream - - //{ - /// ReadFissil : read the zai name in the TMWA MLP model starts with "k_fissil" keyword - /*! - \param line : line suppossed to contain the fissil list - */ - void ReadList(const string &line); - - void ReadEqMaxFraction(const string &line); - void ReadEqMinFraction(const string &line); - bool isIVInDomain(IsotopicVector IV); void StocksTotalMassCalculation(map < string , vector <IsotopicVector> > const& Stocks); - void ConvertMassToLambdaVector(string MaterialDenomination, vector<double>& lambda, double MaterialMassNeeded, vector <IsotopicVector> Stocks); - IsotopicVector BuildFuelToTest(map < string, vector<double> >& lambda, map < string , vector <IsotopicVector> > const& StreamArray, double HMMass, map <string, bool> StreamListIsBuffer); //Build a fuel with the buffer according to fissile lambda - void CheckTargetParameterConsistency(map < int , string > StreamListPriority, map < string , double > TargetParameterMin, map < string , double > TargetParameterMax); - - protected : - - bool fUseTMVAPredictor; //!< Bool that says if we need a TMVA predictor. If not, fuel fraction isimpoased by the FP. - - map < string, IsotopicVector> fStreamList; //!< contains all lists of zai needed to build a fuel (example : 2 -> fissileList+fertileList) - //!< each list is identified by a keyword (example : -> "Fissile" & "Fertile") - map < string , double> fStreamListEqMMassFractionMax; //!< Map that contains lists of stream according to the EqModel with mass maximum fraction - map < string , double> fStreamListEqMMassFractionMin; //!< Map that contains lists of stream according to the EqModel with mass minimum fraction - - double fSpecificPower; //!< The specific power in W/gHM (HM: heavy Metal) - double fMaximalBU; //!< The Maximum burn-up of the model in GWd/t - string fTargetParameter; //!< Type of target parameter optimized in build fuel (ex. BUmax) - int fMaxIterration; //!< Max iterrations in build fueld algorithm + void ConvertMassToLambdaVector(string MaterialDenomination, vector<double>& lambda, double MaterialMassNeeded, vector <IsotopicVector> Stocks); - string fPredictorType ; //!< Type of predictor used in Equivalence Model (ex: MLP) - string fOutput ; //!< Type of output calculated by the predictor - string fBuffer ; //!< Name of material used as buffer in fuel - - map<string, double> fModelParameter ; ///< Map of equivalence model parameter - - map<ZAI,string> fMapOfTMVAVariableNames; //!< List of TMVA input variable names (read from fMLPInformationFile ) , name depends on the training step - - double fTargetParameterStDev; //!< Precision on target parameter calculation +protected : + map < string, IsotopicVector> fStreamList; //!< contains all lists of zai needed to build a fuel (example : 2 -> fissileList+fertileList) + //!< each list is identified by a keyword (example : -> "Fissile" & "Fertile") + double fSpecificPower; //!< The specific power in W/gHM (HM: heavy Metal) void SetLambdaToErrorCode(vector<double>& lambda); - - //@} - - string fTMVAXMLFilePath; //!<The weight needed by TMVA to construct and execute the multilayer perceptron - string fTMVANFOFilePath; //!<The weight needed by TMVA to construct and execute the multilayer perceptron - - #ifndef __ROOTCLING__ map<string, EQM_MthPtr> fKeyword; #endif - - bool freaded; - map< ZAI, pair<double,double> > fZAILimits; //!< Fresh fuel range : map<ZAI<min edge ,max edge >> - string fInformationFile; //!< file containing Reactor Type, Fuel type, HM mass, Power, time vector, and TMVA input variables names (looks the manual for format details) - string fDBFType; //!< Fuel Type (e.g MOX, UOX, ThU, ThPu ...) - string fDBRType; //!< Reactor Type (e.g PWR, FBR-Na, ADS..) + bool freaded; + + map< ZAI, pair<double, double> > fZAILimits; //!< Fresh fuel range : map<ZAI<min edge ,max edge >> /*! - \name Others + \name Others */ //@{ - map <string , double > fTotalMassInStocks; //!< Total mass in each vector of stock + map <string , double > fTotalMassInStocks; //!< Total mass in each vector of stock map <string , double > fLambdaMax; //!< Total lambda of available stocks //@} - private : - - double fPCMprecision; //!< precision on @f$\langle k \rangle@f$ prediction [pcm] +private : }; diff --git a/source/include/XSModel.hxx b/source/include/XSModel.hxx index bc517d3024e080136261f6b5671729effa3228c6..200b482a6f9302369c36386e382755010df44eec 100644 --- a/source/include/XSModel.hxx +++ b/source/include/XSModel.hxx @@ -57,9 +57,10 @@ class XSModel : public CLASSObject XSModel(); //!<Default constructor - XSModel(CLASSLogger* log); //!<Logger constructor + virtual ~XSModel(); //!< Destructor + //@} diff --git a/source/src/EquivalenceModel.cxx b/source/src/EquivalenceModel.cxx index 9f5777b5e238a4cfe4ec07d4fb8975a68e2551e1..9f6c3a68ccaede411ade43551919363dba6c622b 100644 --- a/source/src/EquivalenceModel.cxx +++ b/source/src/EquivalenceModel.cxx @@ -18,1293 +18,146 @@ //________________________________________________________________________ //________________________________________________________________________ -EquivalenceModel::EquivalenceModel(string TMVAXMLFilePath, string TMVANFOFilePath):CLASSObject() +EquivalenceModel::EquivalenceModel(): CLASSObject() { - fUseTMVAPredictor = true; - - fMaxIterration = 500; - freaded = false; - fPCMprecision = 10; - - fTMVAXMLFilePath = TMVAXMLFilePath; - fTMVANFOFilePath = TMVANFOFilePath; - - fDBRType = ""; - fDBFType = ""; - fSpecificPower = 0; - fMaximalBU = 0; - fTargetParameter = ""; - fTargetParameterStDev = 0; - fBuffer = ""; - fPredictorType = ""; - fOutput = ""; - - LoadKeyword(); // Load Key words defineds in NFO file - ReadNFO(); //Getting information from file NFO - - //Check if any information is missing in NFO file - if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); } - if(fMapOfTMVAVariableNames.empty()) { ERROR<<"Missing information for k_zainame in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fTargetParameter.empty()) { ERROR<<"Missing information for k_targetparameter in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(!fTargetParameterStDev) { ERROR<<"Missing information for fTargetParameterStDev in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fModelParameter.empty()) { ERROR<<"Missing information for k_modelparameter in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fPredictorType.empty()) { ERROR<<"Missing information for k_predictortype in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fOutput.empty()) { ERROR<<"Missing information for k_output in : "<<fTMVANFOFilePath<<endl; exit(1);} - - INFO << "__An equivalence model has been define__" << endl; - INFO << "\tThe TMVA weights file is :" << fTMVAXMLFilePath << endl; - INFO << "\tThe TMVA NFO file is :" << fTMVANFOFilePath << endl; - PrintInfo(); - + freaded = false; } //________________________________________________________________________ -EquivalenceModel::EquivalenceModel(CLASSLogger* log, string TMVAXMLFilePath, string TMVANFOFilePath):CLASSObject(log) -{ - fUseTMVAPredictor = true; - - fMaxIterration = 500; - freaded = false; - fPCMprecision = 10; - - fTMVAXMLFilePath = TMVAXMLFilePath; - fTMVANFOFilePath = TMVANFOFilePath; - - fDBRType = ""; - fDBFType = ""; - fSpecificPower = 0; - fMaximalBU = 0; - fTargetParameter = ""; - fTargetParameterStDev = 0; - fBuffer = ""; - fPredictorType = ""; - fOutput = ""; - - LoadKeyword(); // Load Key words defineds in NFO file - ReadNFO(); //Getting information from file NFO - - //Check if any information is missing in NFO file - if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); } - if(fMapOfTMVAVariableNames.empty()) { ERROR<<"Missing information for k_zainame in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fTargetParameter.empty()) { ERROR<<"Missing information for k_targetparameter in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(!fTargetParameterStDev) { ERROR<<"Missing information for fTargetParameterStDev in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fModelParameter.empty()) { ERROR<<"Missing information for k_modelparameter in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fPredictorType.empty()) { ERROR<<"Missing information for k_predictortype in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fOutput.empty()) { ERROR<<"Missing information for k_output in : "<<fTMVANFOFilePath<<endl; exit(1);} - - INFO << "__An equivalence model has been define__" << endl; - INFO << "\tThe TMVA weights file is :" << fTMVAXMLFilePath << endl; - INFO << "\tThe TMVA NFO file is :" << fTMVANFOFilePath << endl; - PrintInfo();} -//________________________________________________________________________ -EquivalenceModel::EquivalenceModel(string TMVANFOFilePath):CLASSObject() +EquivalenceModel::EquivalenceModel(CLASSLogger* log): CLASSObject(log) { - fUseTMVAPredictor = false; - - fTMVANFOFilePath = TMVANFOFilePath; - - fDBRType = ""; - fDBFType = ""; - fSpecificPower = 0; - fMaximalBU = 0; - fTargetParameter = ""; - fTargetParameterStDev = 0; - fBuffer = ""; - - LoadKeyword(); // Load Key words defineds in NFO file - ReadNFO(); //Getting information from file NFO - - //Check if any information is missing in NFO file - if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); } - if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);} - - INFO << "__An equivalence model without TMVA data has been define__" << endl; - INFO << "\tThe NFO file is :" << fTMVANFOFilePath << endl; - PrintInfo(); + freaded = false; } //________________________________________________________________________ -EquivalenceModel::EquivalenceModel(CLASSLogger* log, string TMVANFOFilePath):CLASSObject(log) -{ - fUseTMVAPredictor = false; - - fTMVANFOFilePath = TMVANFOFilePath; - - fDBRType = ""; - fDBFType = ""; - fSpecificPower = 0; - fMaximalBU = 0; - fBuffer = ""; - - LoadKeyword(); // Load Key words defineds in NFO file - ReadNFO(); //Getting information from file NFO - - //Check if any information is missing in NFO file - if(fZAILimits.empty()) {ERROR<<"Missing information for k_zail in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fDBRType.empty()) {ERROR<<"Missing information for k_reactor in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fDBFType.empty()) {ERROR<<"Missing information for k_fuel in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(!fSpecificPower) {ERROR<<"Missing information for k_specpower in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(!fMaximalBU) {ERROR<<"Missing information for k_maxburnup in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fStreamListEqMMassFractionMin.empty() || fStreamListEqMMassFractionMax.empty()) { ERROR<<"Missing information for k_massfractionmin and/or k_massfractionmax in : "<<fTMVANFOFilePath<<endl; exit(1);} - if(fStreamList.empty()) { ERROR<<"Missing information for k_list in : "<<fTMVANFOFilePath<<endl; exit(1); } - if(fBuffer.empty()) { ERROR<<"Missing information for k_buffer in : "<<fTMVANFOFilePath<<endl; exit(1);} - - INFO << "__An equivalence model has been define__" << endl; - INFO << "\tThe TMVA weights file is :" << fTMVAXMLFilePath << endl; - INFO << "\tThe TMVA NFO file is :" << fTMVANFOFilePath << endl; - PrintInfo();} -//________________________________________________________________________ EquivalenceModel::~EquivalenceModel() { - -} -//________________________________________________________________________ -void EquivalenceModel::SetLambdaToErrorCode(vector<double>& lambda) -{ - DBGL - if(lambda.size() == 0) //then we have to add an element to send the error code to the fab (case for no storage in stream) - { - lambda.push_back(-1); - } - - else // other errors (no enough material or too many steps) - { - for( int i=0; i < (int)lambda.size(); i++) - { - lambda[i] = -1; - } - } - DBGL } //________________________________________________________________________ -void EquivalenceModel::StocksTotalMassCalculation(map < string , vector <IsotopicVector> > const& Stocks) +map <string , vector<double> > EquivalenceModel::BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray, map < string , double> StreamListMassFractionMin, map < string , double> StreamListMassFractionMax, map < int , string> StreamListPriority, map < string , bool> StreamListIsBuffer) { - DBGL - // Calculating total mass of stock once and for all - double TotalMassInStocks = 0; - map < string , vector <IsotopicVector> >::const_iterator it_s_vIV; - - for( it_s_vIV = Stocks.begin(); it_s_vIV != Stocks.end(); it_s_vIV++) - { - fTotalMassInStocks[ it_s_vIV->first ] = 0; - fLambdaMax[ it_s_vIV->first ] = 0; - } - for( it_s_vIV = Stocks.begin(); it_s_vIV != Stocks.end(); it_s_vIV++) - { - TotalMassInStocks = 0; - for(int i=0; i < (int)Stocks.at((* it_s_vIV).first).size(); i++) - { - TotalMassInStocks += Stocks.at( it_s_vIV->first )[i].GetTotalMass(); - } - fLambdaMax[(*it_s_vIV).first] = Stocks.at( it_s_vIV->first ).size(); - fTotalMassInStocks[ it_s_vIV->first ] = TotalMassInStocks * 1e6; // in grams - } - DBGL } - //________________________________________________________________________ -void EquivalenceModel::ConvertMassToLambdaVector(string MaterialDenomination, vector<double>& lambda, double MaterialMassNeeded, vector <IsotopicVector> Stocks) -{ - DBGL - double Lambda_tot = 0; - - // Calculation of Lambda tot associated to the required mass MaterialMassNeeded - for( int i=0; i < (int)Stocks.size(); i++) - { - if( MaterialMassNeeded >= (Stocks[i].GetTotalMass()*1e6)) - { - Lambda_tot += 1; - MaterialMassNeeded -= (Stocks[i].GetTotalMass()*1e6); - } - else - { - Lambda_tot += MaterialMassNeeded/(Stocks[i].GetTotalMass()*1e6); - break; - } - } - - // Calculate lambda vector associated to the lambda tot - if(Lambda_tot > (int)lambda.size() ) - { - cout<<Lambda_tot<<" "<<lambda.size()<<endl; - ERROR << " FATAL ERROR " <<endl; - exit(0); - } - - for(int i = 0 ; i < (int)lambda.size() ; i++) //set to 0 all non touched value (to be sure) - lambda[i] = 0 ; - - int IntegerPart = floor( Lambda_tot ); - double DecimalPart = Lambda_tot - IntegerPart; - - for(int i=0 ; i < IntegerPart; i++ ) - lambda[i]=1; - - lambda[IntegerPart] = DecimalPart; - DBGL -} - -//________________________________________________________________________ -IsotopicVector EquivalenceModel::BuildFuelToTest(map < string, vector<double> >& lambda, map < string , vector <IsotopicVector> > const& StreamArray, double HMMass, map <string, bool> StreamListIsBuffer) -{ - DBGL - //Iterators declaration - map < string , vector <IsotopicVector> >::const_iterator it_s_vIV; - map < string , bool >::iterator it_s_B; - - //Find the buffer and set its lambda to 0 - string BufferDenomination =""; - for( it_s_B = StreamListIsBuffer.begin(); it_s_B != StreamListIsBuffer.end(); it_s_B++) - { - if((*it_s_B ).second==true){ BufferDenomination = (*it_s_B).first; } - } - - for(int i = 0; i< lambda[BufferDenomination].size(); i++) - { - lambda[BufferDenomination][i]=0; - } - - //Build an IV with all materials besides buffer to get the total mass of others materials - IsotopicVector IV; - for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++) - { - for(int i=0; i < (int)StreamArray.at( it_s_vIV->first ).size(); i++) - { - IV += lambda[(*it_s_vIV).first][i] * StreamArray.at( it_s_vIV->first )[i]; - } - } - - //Calculate MassBuffer - double MassBuffer = HMMass - IV.GetTotalMass()*1e06; - - //Set buffer lambda according to MassBuffer - ConvertMassToLambdaVector(BufferDenomination, lambda[BufferDenomination], MassBuffer, StreamArray.at(BufferDenomination)); - - IV.Clear(); - - //Build fuel with all materials - for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++) - { - for(int i=0; i < (int)StreamArray.at( it_s_vIV->first ).size(); i++) - { - IV += lambda[(*it_s_vIV).first][i] * StreamArray.at( it_s_vIV->first )[i]; - } - } - DBGL - return IV; - -} - -//________________________________________________________________________ -map <string , vector<double> > EquivalenceModel::BuildFuel(double BurnUp, double HMMass, map < string , vector <IsotopicVector> > StreamArray, map < string , double> StreamListFPMassFractionMin, map < string , double> StreamListFPMassFractionMax, map < int , string > StreamListPriority, map < string , bool> StreamListIsBuffer) +void EquivalenceModel::SetLambdaToErrorCode(vector<double>& lambda) { - DBGL - - HMMass *= 1e6; //Unit conversion : tons to gram - - map <string , vector<double> > lambda ; // map containing name of the list and associated vector of proportions taken from stocks - //Iterators declaration - map < string , vector <IsotopicVector> >::iterator it_s_vIV; - map < string , vector <double> >::iterator it_s_vD; - map < string , IsotopicVector >::iterator it_s_IV; - map < string , double >::iterator it_s_D; - map < int , string >::iterator it_i_s; - - // Initialize lambda to 0 // - for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++) - { - for(int i=0; i < (int)StreamArray[(*it_s_vIV).first].size(); i++) - { - lambda[(*it_s_vIV).first].push_back(0); - } + DBGL + if (lambda.size() == 0) //then we have to add an element to send the error code to the fab (case for no storage in stream) + { + lambda.push_back(-1); } - // Test if there is at least one stock available in each list, otherwise fuel is not built // - bool BreakReturnLambda = false; - for( it_s_vIV = StreamArray.begin(); it_s_vIV != StreamArray.end(); it_s_vIV++) + else // other errors (no enough material or too many steps) { - if(StreamArray[(*it_s_vIV).first].size() == 0) + for ( int i = 0; i < (int)lambda.size(); i++) { - WARNING << " No stock available for stream : "<< (*it_s_vIV).first <<". Fuel not built." << endl; - SetLambdaToErrorCode(lambda[(*it_s_vIV).first]); - BreakReturnLambda = true; + lambda[i] = -1; } } - if(BreakReturnLambda) { return lambda;} + DBGL +} +//________________________________________________________________________ +void EquivalenceModel::StocksTotalMassCalculation(map < string , vector <IsotopicVector> > const& Stocks) +{ + DBGL + // Calculating total mass of stock once and for all + double TotalMassInStocks = 0; + map < string , vector <IsotopicVector> >::const_iterator it_s_vIV; - // Check if the targeted burn-up is lower than maximum burn-up of model // - if(BurnUp > fMaximalBU) + for ( it_s_vIV = Stocks.begin(); it_s_vIV != Stocks.end(); it_s_vIV++) { - ERROR << " Targeted burn-up is higher than maximum burn-up defined in NFO file..."<< endl; - ERROR << " Targeted burn-up : "<<BurnUp<<" GWd/t"<<endl; - ERROR << " Maximum burn-up : "<<fMaximalBU<<" GWd/t"<<endl; - exit(1); + fTotalMassInStocks[ it_s_vIV->first ] = 0; + fLambdaMax[ it_s_vIV->first ] = 0; } - -// Fissile fraction calculation is needed. - if (fUseTMVAPredictor) + for ( it_s_vIV = Stocks.begin(); it_s_vIV != Stocks.end(); it_s_vIV++) { - // Check if EquivalenceModel->SetTMVAXMLFilePath() and/or EquivalenceModel->SetTMVANFOFilePath() have been defined - if (fTMVAXMLFilePath.empty() || fTMVANFOFilePath.empty()) + TotalMassInStocks = 0; + for (int i = 0; i < (int)Stocks.at((* it_s_vIV).first).size(); i++) { - ERROR << " TMVA XML and/or NFO File path are not defined..."<< endl; - ERROR << " You have to use EquivalenceModel->SetTMVAXMLFilePath() and/or EquivalenceModel->SetTMVANFOFilePath() methods."<<endl; - exit(1); + TotalMassInStocks += Stocks.at( it_s_vIV->first )[i].GetTotalMass(); } - - double TargetParameterValue = 0; - - for( it_s_D = fModelParameter.begin(); it_s_D != fModelParameter.end(); it_s_D++) - { - if(fModelParameter[(*it_s_D).first] == -1) - { - ERROR<< "Model parameter ( "<<fModelParameter[(*it_s_D).first] << " ) value is not defined in the input." <<endl; - ERROR<< "Use EqM->SetModelParameter( \" "<<(*it_s_D).first<<" \", value) to define it." <<endl; - exit(1); - } - } - - if(fTargetParameter=="BurnUpMax") {TargetParameterValue = BurnUp;} - else if (fTargetParameter=="keffBOC") {TargetParameterValue = fModelParameter["keffBOC"];} - else - { - ERROR<< "Target parameter defined in InformationFile ( "<<fTargetParameter<<" ) doesn't exist." <<endl; - ERROR<< "Possible target parameters for the moment are : "<< endl; - ERROR<< " - BurnUpMax - Used for PWR" <<endl; - ERROR<< " - keffBOC - Used for SFR" <<endl; - exit(1); - } - - /// Search for the minimum and maximum fraction of each material in fuel /// - map < string, double > StreamListMassFractionMin ; - map < string, double > StreamListMassFractionMax ; - for( it_s_D = StreamListFPMassFractionMin.begin(); it_s_D != StreamListFPMassFractionMin.end(); it_s_D++) - { - if(StreamListFPMassFractionMin[(*it_s_D).first] < fStreamListEqMMassFractionMin[(*it_s_D).first]) // if limits FP are lower than limits EqM - { - ERROR << " User mass fraction min requirement is lower than the model mass fraction min for list : "<<(*it_s_D).first << endl; - ERROR << " User mass fraction min requirement : "<<StreamListFPMassFractionMin[(*it_s_D).first]<<endl; - ERROR << " Model mass fraction min requirement : "<<fStreamListEqMMassFractionMin[(*it_s_D).first]<<endl; - exit(1); - } - else - { - StreamListMassFractionMin[(*it_s_D).first] = StreamListFPMassFractionMin[(*it_s_D).first]; - } - } - - for( it_s_D = StreamListFPMassFractionMax.begin(); it_s_D != StreamListFPMassFractionMax.end(); it_s_D++) - { - if(StreamListFPMassFractionMax[(*it_s_D).first] > fStreamListEqMMassFractionMax[(*it_s_D).first]) // if limits FP are higher than limits EqM - { - ERROR << " User mass fraction max requirement is higher than the model mass fraction max for list : "<<(*it_s_D).first << endl; - ERROR << " User mass fraction max requirement : "<<StreamListFPMassFractionMax[(*it_s_D).first]<<endl; - ERROR << " Model mass fraction max requirement : "<<fStreamListEqMMassFractionMax[(*it_s_D).first]<<endl; - exit(1); - } - else - { - StreamListMassFractionMax[(*it_s_D).first] = StreamListFPMassFractionMax[(*it_s_D).first]; - } - - } - - //Calculate Total mass in stock for each stream and fill fTotalMassInStocks - StocksTotalMassCalculation(StreamArray); - - // Check if there is enough material in stock to satisfy mass fraction min // - BreakReturnLambda = false; - for( it_s_D = StreamListMassFractionMin.begin(); it_s_D != StreamListMassFractionMin.end(); it_s_D++) - { - if(fTotalMassInStocks[(*it_s_D).first]< HMMass*StreamListMassFractionMin[(*it_s_D).first]) - { - WARNING << " Not enough material : "<< (*it_s_D).first << " in stocks to reach the build fuel lower limit of "<<StreamListMassFractionMin[(*it_s_D).first]<<" reactor mass. Fuel not built." << endl; - SetLambdaToErrorCode(lambda[(*it_s_D).first]); - BreakReturnLambda = true; - } - } - if(BreakReturnLambda) { return lambda;} - - // Check if there is enough material in stock to satisfy mass fraction max, if not mass fraction max is set to MassINStock/MassReactor// - for( it_s_D = StreamListMassFractionMax.begin(); it_s_D != StreamListMassFractionMax.end(); it_s_D++) - { - if(fTotalMassInStocks[(*it_s_D).first]< HMMass*StreamListMassFractionMax[(*it_s_D).first]) - { - StreamListMassFractionMax[(*it_s_D).first] = fTotalMassInStocks[(*it_s_D).first]/HMMass; - WARNING << " Not enough material : "<< (*it_s_D).first << " in stocks to reach the build fuel higher limit of "<<StreamListMassFractionMax[(*it_s_D).first]<<" reactor mass. " << endl; - WARNING << " Mass fraction max of material : "<< (*it_s_D).first << " is set to MassInStock/HMMassReactor : "<< StreamListMassFractionMax[(*it_s_D).first]<< endl; - } - } - - //Check if TargetParameter is inside [TargetParameterMin, TargetParameterMax] associated to fraction Min et Max// - - map < string , double > MassMin; - map < string , double > MassMax; - - map < string , double > TargetParameterMin; - map < string , double > TargetParameterMax; - - IsotopicVector FuelToTest; - - bool TargetParameterIncluded = false; - for( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++) - { - //Calculate TargetParameterMin for each possibility : min1 ; max1 + min2 ; max1 + max2 + min3 .... - MassMin[(*it_i_s ).second] = HMMass * StreamListMassFractionMin[(*it_i_s).second]; - ConvertMassToLambdaVector((*it_i_s ).second, lambda[(*it_i_s ).second], MassMin[(*it_i_s ).second], StreamArray[(*it_i_s ).second]); - FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer); - FuelToTest = FuelToTest/FuelToTest.GetSumOfAll(); - TargetParameterMin[(*it_i_s ).second] = CalculateTargetParameter(FuelToTest, fTargetParameter); - - //Check is TargetParameterMin < TargetParameter - if(TargetParameterMin[(*it_i_s ).second]>TargetParameterValue) - { - if((*it_i_s).first ==1) //Minimum of first material is too high - { - WARNING << "CRITICAL ! Minimum parameter value associated to the first priority material ( "<<(*it_i_s ).second <<" ) is higher than targeted parameter."<< endl; - WARNING << "Targeted parameter : "<<fTargetParameter<<" = "<<TargetParameterValue<<endl; - WARNING << "Minimum parameter value : " <<TargetParameterMin[(*it_i_s ).second]<<endl; - WARNING << "Try to increase targeted parameter." <<endl; - SetLambdaToErrorCode(lambda[(*it_i_s).second]); - return lambda; - DBGL - } - else if ((*it_i_s).first >1) //TargetParameter is located between max n-1 and min n - { - WARNING << "CRITICAL ! Targeted parameter value ( "<<fTargetParameter<<" ) is located between 2 materials. "<<endl; - it_i_s --; - WARNING << fTargetParameter <<" of max fraction of material : "<< (*it_i_s).second<<" ---> "<<TargetParameterMax[(*it_i_s ).second]<<endl; - it_i_s ++; - WARNING << fTargetParameter<< " of min fraction of material : "<< (*it_i_s ).second<<" ---> "<<TargetParameterMin[(*it_i_s ).second]<<endl; - WARNING << "Targeted "<<fTargetParameter<<" : " <<TargetParameterValue<<endl; - WARNING << "Try to decrease mimimum fraction of : "<< (*it_i_s ).second<<endl; - SetLambdaToErrorCode(lambda[(*it_i_s).second]); - return lambda; - } - } - FuelToTest.Clear(); - - //Calculate TargetParameter max for each possibility : max1 ; max1 + max2 ; max1 + max2 + max3 .... - MassMax[(*it_i_s ).second] = HMMass * StreamListMassFractionMax[(*it_i_s).second]; - ConvertMassToLambdaVector((*it_i_s ).second, lambda[(*it_i_s ).second], MassMax[(*it_i_s ).second], StreamArray[(*it_i_s ).second]); - FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer); - FuelToTest = FuelToTest/FuelToTest.GetSumOfAll(); - TargetParameterMax[(*it_i_s ).second] = CalculateTargetParameter(FuelToTest, fTargetParameter); - - if(TargetParameterMax[(*it_i_s ).second]>=TargetParameterValue) - { - TargetParameterIncluded = true ; - break; - } - } - - //Check if target parameter increases monotously with the material mass - CheckTargetParameterConsistency(StreamListPriority, TargetParameterMin, TargetParameterMax); - - if(!TargetParameterIncluded) - { - WARNING << "CRITICAL ! Maximum reachable "<<fTargetParameter<<" is lower than targeted "<< fTargetParameter<<". "<< endl; - WARNING << "Targeted "<<fTargetParameter<<" = "<<TargetParameterValue<<endl; - WARNING << "Maximum reachable "<<fTargetParameter<<" : "<<TargetParameterMax[(*--StreamListPriority.end()).second]<<endl; - WARNING << "Try to increase maximum fraction of materials, or decrease "<< fTargetParameter<<" ." <<endl; - SetLambdaToErrorCode(lambda[(*--StreamListPriority.end()).second]); - return lambda; - } - - //Search the TargetParameterValue location in the mass damain // - string MaterialToSearch = (*it_i_s ).second; - double CalculatedTargetParameter = TargetParameterMax[MaterialToSearch] ; //Algo start with maximum point - double MassToAdd = MassMax[MaterialToSearch]; //Algo start with maximum point - - double LastMassMinus = MassMin[MaterialToSearch]; //Used in bissection method - double LastMassPlus = MassMax[MaterialToSearch]; //Used in bissection method - - int count = 0; - - FuelToTest.Clear(); - - /* - if (fDBFType == "MOX") - { - cout<<"------------------------------------------------------"<<endl; - cout<<"START ALGO -> BU, Mass "<<BurnUp<<" "<<HMMass<<endl; - cout<<"------------------------------------------------------"<<endl; - double MassTest = MassMin[MaterialToSearch]; - cout<<MaterialToSearch<<" "<<MassMax[MaterialToSearch]<<" "<<MassMin[MaterialToSearch]<<" "<<endl; - do - { - ConvertMassToLambdaVector(MaterialToSearch, lambda[MaterialToSearch], MassTest, StreamArray[MaterialToSearch]); - FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer); - FuelToTest = FuelToTest/FuelToTest.GetSumOfAll(); - CalculatedTargetParameter = CalculateTargetParameter(FuelToTest, fTargetParameter); - - cout<<"Lambda vector : "<<MaterialToSearch<<" - "; for(int i=0; i < (int)lambda[MaterialToSearch].size(); i++) cout<<lambda[MaterialToSearch][i]<<" "; - cout<<endl; - - - MassTest += (MassMax[MaterialToSearch] - MassMin[MaterialToSearch])/100.; - - cout<<MassTest<<" "<<CalculatedTargetParameter<<endl; - - } while (MassTest <= MassMax[MaterialToSearch]); - cout<<"------------------------------------------------------"<<endl; - cout<<"STOP ALGO EXIT(1)..."<<endl; exit(1); - cout<<"------------------------------------------------------"<<endl; - } - */ - - do - { - if(count > fMaxIterration) - { - ERROR << "CRITICAL ! Can't manage to predict fissile content\nHint : Try to decrease the precision on the target parameter using :\nYourEquivalenceModel->SetTargetParameterStDev(Precision); " << endl; - ERROR << "Targeted "<<fTargetParameter<<" : "<<TargetParameterValue<<endl; - ERROR << "Last calculated "<<fTargetParameter<<" : "<<CalculatedTargetParameter<<endl; - ERROR << "Last Fresh fuel normalized composition : " <<endl; - ERROR << FuelToTest.sPrint()<<endl; - exit(1); - } - - if( (CalculatedTargetParameter - TargetParameterValue) < 0 ) //Need to add more fissile material in fuel - { - LastMassMinus = MassToAdd; - MassToAdd = MassToAdd + fabs(LastMassPlus - MassToAdd)/2.; - } - else if( (CalculatedTargetParameter - TargetParameterValue) > 0) //Need to add less fissile material in fuel - { - LastMassPlus = MassToAdd; - MassToAdd = MassToAdd - fabs(LastMassMinus - MassToAdd)/2.; - } - ConvertMassToLambdaVector(MaterialToSearch, lambda[MaterialToSearch], MassToAdd, StreamArray[MaterialToSearch]); - FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer); - FuelToTest = FuelToTest/FuelToTest.GetSumOfAll(); - CalculatedTargetParameter = CalculateTargetParameter(FuelToTest, fTargetParameter); - - count ++; - - }while(fabs(TargetParameterValue - CalculatedTargetParameter) > GetTargetParameterStDev()*TargetParameterValue); + fLambdaMax[(*it_s_vIV).first] = Stocks.at( it_s_vIV->first ).size(); + fTotalMassInStocks[ it_s_vIV->first ] = TotalMassInStocks * 1e6; // in grams } -// Fissile fraction is imposed by the FP -// No need to use algo -// Simplified fuel building - else - { - // Check if EquivalenceModel->SetTMVANFOFilePath() have been defined - if (fTMVANFOFilePath.empty()) - { - ERROR << " TMVA NFO File path is not defined..."<< endl; - ERROR << " You have to use EquivalenceModel->SetTMVANFOFilePath() methods."<<endl; - exit(1); - } - - /// Search for the fraction of each material in fuel /// - map < string, double > StreamListMassFraction; - for( it_s_D = StreamListFPMassFractionMin.begin(); it_s_D != StreamListFPMassFractionMin.end(); it_s_D++) - { - if(StreamListFPMassFractionMin[(*it_s_D).first] < fStreamListEqMMassFractionMin[(*it_s_D).first]) // if limits FP are lower than limits EqM - { - ERROR << " User mass fraction requirement is lower than the model mass fraction min for list : "<<(*it_s_D).first << endl; - ERROR << " User mass fraction requirement : "<<StreamListFPMassFractionMin[(*it_s_D).first]<<endl; - ERROR << " Model mass fraction min requirement : "<<fStreamListEqMMassFractionMin[(*it_s_D).first]<<endl; - exit(1); - } - else if(StreamListFPMassFractionMax[(*it_s_D).first] > fStreamListEqMMassFractionMax[(*it_s_D).first]) // if limits FP are higher than limits EqM - { - ERROR << " User mass fraction requirement is higher than the model mass fraction max for list : "<<(*it_s_D).first << endl; - ERROR << " User mass fraction requirement : "<<StreamListFPMassFractionMax[(*it_s_D).first]<<endl; - ERROR << " Model mass fraction max requirement : "<<fStreamListEqMMassFractionMax[(*it_s_D).first]<<endl; - exit(1); - } - else - { - StreamListMassFraction[(*it_s_D).first] = StreamListFPMassFractionMin[(*it_s_D).first]; // Because here, min = max - } - } - - //Calculate Total mass in stock for each stream and fill fTotalMassInStocks - StocksTotalMassCalculation(StreamArray); - - // Check if there is enough material in stock to satisfy requested mass fraction // - BreakReturnLambda = false; - for( it_s_D = StreamListMassFraction.begin(); it_s_D != StreamListMassFraction.end(); it_s_D++) - { - if(fTotalMassInStocks[(*it_s_D).first]< HMMass*StreamListMassFraction[(*it_s_D).first]) - { - WARNING << " Not enough material : "<< (*it_s_D).first << " in stocks to reach the build fuel limit of "<<StreamListMassFraction[(*it_s_D).first]<<" reactor mass. Fuel not built." << endl; - SetLambdaToErrorCode(lambda[(*it_s_D).first]); - BreakReturnLambda = true; - } - } - if(BreakReturnLambda) { return lambda;} + DBGL +} - IsotopicVector FuelToTest; - // Build Fuel - for( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++) - { - FuelToTest.Clear(); - ConvertMassToLambdaVector((*it_i_s).second, lambda[(*it_i_s).second], HMMass*StreamListMassFraction[(*it_i_s).second], StreamArray[(*it_i_s).second]); - FuelToTest = BuildFuelToTest(lambda, StreamArray, HMMass, StreamListIsBuffer); - } - } +//________________________________________________________________________ +void EquivalenceModel::ConvertMassToLambdaVector(string MaterialDenomination, vector<double>& lambda, double MaterialMassNeeded, vector <IsotopicVector> Stocks) +{ + DBGL + double Lambda_tot = 0; - //Final builded fuel - IsotopicVector IVStream; - for( it_s_vD = lambda.begin(); it_s_vD != lambda.end(); it_s_vD++) + // Calculation of Lambda tot associated to the required mass MaterialMassNeeded + for ( int i = 0; i < (int)Stocks.size(); i++) { - for(int i=0; i<(int)lambda[(*it_s_vD).first].size(); i++) + if ( MaterialMassNeeded >= (Stocks[i].GetTotalMass() * 1e6)) { - IVStream +=lambda[(*it_s_vD).first][i] * StreamArray[(*it_s_vD).first][i]; + Lambda_tot += 1; + MaterialMassNeeded -= (Stocks[i].GetTotalMass() * 1e6); } - } - - //Check if BuildedFuel is in Model isotopic bounds - (*this).isIVInDomain(IVStream); - - for( it_s_vD = lambda.begin(); it_s_vD != lambda.end(); it_s_vD++) - { - DBGV( "Lambda vector : "<<(*it_s_vD).first ); - - for(int i=0; i < (int)lambda[(*it_s_vD).first].size(); i++) + else { - DBGV(lambda[(*it_s_vD).first][i]); + Lambda_tot += MaterialMassNeeded / (Stocks[i].GetTotalMass() * 1e6); + break; } } - - DBGL - return lambda; -} - -//________________________________________________________________________ -TTree* EquivalenceModel::CreateTMVAInputTree(IsotopicVector TheFreshfuel, double ThisTime) -{ - /******Create Input data tree to be interpreted by TMVA::Reader***/ - TTree* InputTree = new TTree(fOutput.c_str(), fOutput.c_str()); - - vector<float> InputTMVA; - for(int i = 0 ; i< (int)fMapOfTMVAVariableNames.size() ; i++) - InputTMVA.push_back(0); - - float Time = 0; - - IsotopicVector IVInputTMVA; - map<ZAI ,string >::iterator it_ZAI_s; - int j = 0; - - for( it_ZAI_s = fMapOfTMVAVariableNames.begin() ; it_ZAI_s != fMapOfTMVAVariableNames.end() ; it_ZAI_s++) - { - InputTree->Branch( ((*it_ZAI_s).second).c_str(), &InputTMVA[j], ((*it_ZAI_s).second + "/F").c_str()); - IVInputTMVA+= ((*it_ZAI_s).first)*1; - j++; - } - - if(ThisTime != -1) - InputTree->Branch("Time" ,&Time ,"Time/F"); - - IsotopicVector IVAccordingToUserInfoFile = TheFreshfuel.GetThisComposition(IVInputTMVA); - double Ntot = IVAccordingToUserInfoFile.GetSumOfAll(); - IVAccordingToUserInfoFile = IVAccordingToUserInfoFile/Ntot; - - j = 0; - - for( it_ZAI_s = fMapOfTMVAVariableNames.begin() ; it_ZAI_s != fMapOfTMVAVariableNames.end() ; it_ZAI_s++) - { - InputTMVA[j] = IVAccordingToUserInfoFile.GetZAIIsotopicQuantity( (*it_ZAI_s).first ) ; - j++; - } - - Time = ThisTime; - InputTree->Fill(); - - return InputTree; - -} -//________________________________________________________________________ -void EquivalenceModel::CheckTargetParameterConsistency(map < int , string > StreamListPriority, map < string , double > TargetParameterMin, map < string , double > TargetParameterMax) -{ - map < int , string >::iterator it_i_s; + // Calculate lambda vector associated to the lambda tot + if (Lambda_tot > (int)lambda.size() ) + { + cout << Lambda_tot << " " << lambda.size() << endl; + ERROR << " FATAL ERROR " << endl; + exit(0); + } - //Loop on priority order to check if target parameter increases monotously with the material mass - for( it_i_s = StreamListPriority.begin(); it_i_s != StreamListPriority.end(); it_i_s++) - { - double TargetParameterUp = -1.0; //to be sure BUMin is > to BUmax even if BUmin is zero - double TargetParameterDown = 0.0; + for (int i = 0 ; i < (int)lambda.size() ; i++) //set to 0 all non touched value (to be sure) + lambda[i] = 0 ; - if(TargetParameterMin.find((*it_i_s).second) == TargetParameterMin.end()) - { - break; //if material is not in map, break the loop - } - TargetParameterDown = TargetParameterMin[(*it_i_s).second]; + int IntegerPart = floor( Lambda_tot ); + double DecimalPart = Lambda_tot - IntegerPart; - if (TargetParameterDown < 0.0 ) - { - ERROR<< "Target parameter evolution should always be positive." <<endl; - ERROR<< "TargetParameterDown = "<< TargetParameterDown<<" is negative "<<endl; - ERROR<< "Check the evolution..." <<endl; - exit(1); - } - TargetParameterUp = TargetParameterMax[(*it_i_s).second]; + for (int i = 0 ; i < IntegerPart; i++ ) + lambda[i] = 1; - if (TargetParameterDown > TargetParameterUp ) - { - ERROR<< "Target parameter evolution as a function of material mass is not monotonous." <<endl; - ERROR<< "TargetParameterDown = "<< TargetParameterDown<<" is greater than TargetParameterUp = "<< TargetParameterUp<<endl; - ERROR<< "Check the evolution..." <<endl; - exit(1); - } - } -} -//________________________________________________________________________ -double EquivalenceModel::CalculateTargetParameter(IsotopicVector TheFuel, string TargetParameterName) -{ - double ParameterToCalculate = 0; - if(TargetParameterName=="BurnUpMax") ParameterToCalculate = CalculateBurnUpMax(TheFuel, fModelParameter); - else if(TargetParameterName=="keffBOC") ParameterToCalculate = CalculateKeffAtBOC(TheFuel); - else - { - ERROR<< "Target parameter defined in InformationFile ( "<<TargetParameterName<<" ) doesn't exist" <<endl; - ERROR<< "Possible target parameters for the moment are : BurnUpMax and keffBOC." <<endl; - exit(1); - } - - return ParameterToCalculate ; -} -//________________________________________________________________________ -double EquivalenceModel::CalculateBurnUpMax(IsotopicVector TheFuel, map<string, double> ModelParameter) -{ - /**************************************************************************/ - //With a dichotomy, the maximal irradiation time (TheFinalTime) is calculated - //When average Kinf is very close (according "Precision") to the threshold - //then the corresponding irradiation time is convert in burnup and returned - /**************************************************************************/ - //Algorithm initialization - double KThreshold = fModelParameter["kThreshold"]; - int NumberOfBatch = (int)fModelParameter["NumberOfBatch"]; - double OldFinalTimeMinus = 0; - double MaximumBU = fMaximalBU; - double MinimumBU = 0 ; - double TheFinalTime = BurnupToSecond((MaximumBU-MinimumBU)/2.); - double OldFinalTimePlus = BurnupToSecond(MaximumBU); - double k_av = 0; //average kinf - double OldPredictedk_av = 0; - - CLASSReader * reader = new CLASSReader( fMapOfTMVAVariableNames ); - reader->AddVariable( "Time" ); - reader->BookMVA( "MLP method" , fTMVAXMLFilePath ); - - for(int b = 0;b<NumberOfBatch;b++) - { - float TheTime = (b+1)*TheFinalTime/NumberOfBatch; - - TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime); - reader->SetInputData( InputTree ); - - OldPredictedk_av += reader->EvaluateRegression( "MLP method" )[0]; - - delete InputTree; - } - OldPredictedk_av /= NumberOfBatch; - - //Algorithm control - int count = 0; - int MaximumLoopCount = 500; - do - { - if(count > MaximumLoopCount ) - { - ERROR << "CRITICAL ! Can't manage to predict burnup\nHint : Try to increase the precision on k effective using :\n YourEQM_MLP_Kinf->SetPCMPrecision(pcm); with pcm the precision in pcm (default 10) REDUCE IT\n If this message still appear mail to leniau@subatech.in2p3.fr\nor nicolas.thiolliere@subatech.in2p3.fr " << endl; - exit(1); - } - - if( (OldPredictedk_av-KThreshold) > 0) //The burnup can be increased - { - OldFinalTimeMinus = TheFinalTime; - TheFinalTime = TheFinalTime + fabs(OldFinalTimePlus - TheFinalTime)/2.; - - if(SecondToBurnup(TheFinalTime) >= (MaximumBU-MaximumBU*GetTargetParameterStDev() ) ) - { delete reader; return MaximumBU; } - } - - else if( (OldPredictedk_av-KThreshold) < 0)//The burnup is too high - { - OldFinalTimePlus = TheFinalTime; - TheFinalTime = TheFinalTime - fabs(OldFinalTimeMinus-TheFinalTime)/2.; - if( SecondToBurnup(TheFinalTime) < (MaximumBU-MinimumBU)/2.*GetTargetParameterStDev() ) - { delete reader; return 0; } - } - - k_av = 0; - for(int b = 0;b<NumberOfBatch;b++) - { - float TheTime = (b+1)*TheFinalTime/NumberOfBatch; - TTree* InputTree = CreateTMVAInputTree(TheFuel,TheTime); - reader->SetInputData( InputTree ); - - k_av += reader->EvaluateRegression("MLP method")[0]; - - delete InputTree; - } - k_av/= NumberOfBatch; - //cout<<SecondToBurnup(TheFinalTime)<<" "; - OldPredictedk_av = k_av; - count++; -//std::clog << "-> " << k_av << "\t\t(" << count << ") \t [" << TheFinalTime << "]" << "\t" << OldPredictedk_av-KThreshold << "\t" << GetPCMPrecision() << std::endl; - } while( fabs(OldPredictedk_av-KThreshold) > GetPCMPrecision() ) ; - - delete reader; - //cout<<endl; - return SecondToBurnup(TheFinalTime); + lambda[IntegerPart] = DecimalPart; + DBGL } - -//________________________________________________________________________ -double EquivalenceModel::CalculateKeffAtBOC(IsotopicVector FreshFuel) -{ - CLASSReader * reader = new CLASSReader( fMapOfTMVAVariableNames ); - reader->BookMVA( "MLP method" , fTMVAXMLFilePath ); - - TTree* InputTree = CreateTMVAInputTree(FreshFuel,-1) ; - reader->SetInputData( InputTree ); - - double keff = reader->EvaluateRegression( "MLP method" )[0]; - - delete InputTree; - - return keff; -} - //________________________________________________________________________ bool EquivalenceModel::isIVInDomain(IsotopicVector IV) { - DBGL - bool IsInDomain = true; - - if(fZAILimits.empty()) - { - WARNING << "Fresh Fuel variation domain is not set" << endl; - WARNING << "CLASS has no clue if the computed evolution for this fresh fuel is correct" << endl; - WARNING << "Proceed finger crossed !!" << endl; - return true; - } - - else - { - IsotopicVector IVNorm = IV /IV.GetSumOfAll(); - for (map< ZAI,pair<double,double> >::iterator Domain_it = fZAILimits.begin(); Domain_it != fZAILimits.end(); Domain_it++) - { - double ThatZAIProp = IVNorm.GetIsotopicQuantity()[Domain_it->first]; - double ThatZAIMin = Domain_it->second.first; - double ThatZAIMax = Domain_it->second.second; - if( (ThatZAIProp > ThatZAIMax) || (ThatZAIProp < ThatZAIMin) ) - { - IsInDomain = false; - - WARNING << "Fresh fuel out of model range" << endl; - WARNING << "\t AT LEAST this ZAI is accused to be outrange :" << endl; - WARNING << "\t\t" << Domain_it->first.Z() << " " << Domain_it->first.A() << " " << Domain_it->first.I() << endl; - WARNING << "\t\t min = " << ThatZAIMin << " value = " << ThatZAIProp << " max = " << ThatZAIMax << endl; - WARNING << "\t IV accused :" << endl << endl; - WARNING << IVNorm.sPrint() << endl; - break; - } - } - } - DBGL - return IsInDomain; - -} -//________________________________________________________________________ -void EquivalenceModel::ReadNFO() -{ - DBGL - ifstream NFO(fTMVANFOFilePath.c_str()); - - if(!NFO) - { - ERROR << "Can't find/open file " << fTMVANFOFilePath << endl; - exit(0); - } - - do - { - string line; - getline(NFO,line); - - EquivalenceModel::ReadLine(line); - - } while(!NFO.eof()); - - DBGL -} -//________________________________________________________________________ -void EquivalenceModel::ReadLine(string line) -{ - DBGL - - if (!freaded) - { - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - - map<string, EQM_MthPtr>::iterator it = fKeyword.find(keyword); - - if(it != fKeyword.end()) - (this->*(it->second))( line ); - - freaded = true; - ReadLine(line); - - } - - freaded = false; - - DBGL -} -//________________________________________________________________________ -void EquivalenceModel::LoadKeyword() -{ - DBGL - fKeyword.insert( pair<string, EQM_MthPtr>( "k_zail", & EquivalenceModel::ReadZAIlimits) ); - fKeyword.insert( pair<string, EQM_MthPtr>( "k_reactor", & EquivalenceModel::ReadType) ); - fKeyword.insert( pair<string, EQM_MthPtr>( "k_fuel", & EquivalenceModel::ReadType) ); - fKeyword.insert( pair<string, EQM_MthPtr>( "k_massfractionmin", & EquivalenceModel::ReadEqMinFraction) ); - fKeyword.insert( pair<string, EQM_MthPtr>( "k_massfractionmax", & EquivalenceModel::ReadEqMaxFraction) ); - fKeyword.insert( pair<string, EQM_MthPtr>( "k_list", & EquivalenceModel::ReadList) ); - fKeyword.insert( pair<string, EQM_MthPtr>( "k_specpower", & EquivalenceModel::ReadSpecificPower) ); - if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_zainame", & EquivalenceModel::ReadZAIName) ); - fKeyword.insert( pair<string, EQM_MthPtr>( "k_maxburnup", & EquivalenceModel::ReadMaxBurnUp) ); - if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_targetparameter", & EquivalenceModel::ReadTargetParameter) ); - if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_predictortype", & EquivalenceModel::ReadPredictorType) ); - if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_output", & EquivalenceModel::ReadOutput) ); - fKeyword.insert( pair<string, EQM_MthPtr>( "k_buffer", & EquivalenceModel::ReadBuffer) ); - if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_modelparameter", & EquivalenceModel::ReadModelParameter) ); - if (fUseTMVAPredictor) fKeyword.insert( pair<string, EQM_MthPtr>( "k_targetparameterstdev", & EquivalenceModel::ReadTargetParameterStDev) ); - - DBGL -} -//________________________________________________________________________ -void EquivalenceModel::ReadType(const string &line) -{ - DBGL - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_fuel" && keyword != "k_reactor" ) // Check the keyword - { - ERROR << " Bad keyword : " << keyword << " Not found !" << endl; - exit(1); - } - if( keyword == "k_fuel" ) - fDBFType = StringLine::NextWord(line, pos, ' '); - else if( keyword == "k_reactor" ) - fDBRType = StringLine::NextWord(line, pos, ' '); - - DBGL -} -//________________________________________________________________________ -void EquivalenceModel::ReadZAIlimits(const string &line) -{ - DBGL - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_zail" ) // Check the keyword - { - ERROR << " Bad keyword : \"k_zail\" 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()); - - double downLimit = atof(StringLine::NextWord(line, pos, ' ').c_str()); - double upLimit = atof(StringLine::NextWord(line, pos, ' ').c_str()); - - if (upLimit < downLimit) - { - double tmp = upLimit; - upLimit = downLimit; - downLimit = tmp; - } - fZAILimits.insert(pair<ZAI, pair<double, double> >(ZAI(Z,A,I), pair<double,double>(downLimit, upLimit))); - DBGL -} -//________________________________________________________________________ -void EquivalenceModel::ReadList(const string &line) -{ - DBGL - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_list" ) // Check the keyword - { - ERROR << " Bad keyword : \"k_list\" not found !" << endl; - exit(1); - } - string ListName= StringLine::NextWord(line, pos, ' '); - 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()); - double Q = atof(StringLine::NextWord(line, pos, ' ').c_str()); - fStreamList[ListName].Add(Z, A, I, Q); - - DBGL -} -//________________________________________________________________________ -void EquivalenceModel::ReadEqMinFraction(const string &line) -{ - DBGL - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_massfractionmin" ) // Check the keyword - { - ERROR << " Bad keyword : \"k_massfractionmin\" not found !" << endl; - exit(1); - } - string ListName= StringLine::NextWord(line, pos, ' '); - double Q = atof(StringLine::NextWord(line, pos, ' ').c_str()); - fStreamListEqMMassFractionMin[ListName] = Q; - - DBGL -} - -//________________________________________________________________________ -void EquivalenceModel::ReadEqMaxFraction(const string &line) -{ - DBGL - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_massfractionmax" ) // Check the keyword - { - ERROR << " Bad keyword : \"k_massfractionmax\" not found !" << endl; - exit(1); - } - string ListName= StringLine::NextWord(line, pos, ' '); - double Q = atof(StringLine::NextWord(line, pos, ' ').c_str()); - fStreamListEqMMassFractionMax[ListName] = Q; - - DBGL -} - -//________________________________________________________________________ -void EquivalenceModel::ReadSpecificPower(const string &line) -{ - DBGL - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_specpower") // Check the keyword - { - ERROR << " Bad keyword : \"k_specpower\" Not found !" << endl; - exit(1); - } - - fSpecificPower = atof(StringLine::NextWord(line, pos, ' ').c_str()); - - DBGL -} -//________________________________________________________________________ -void EquivalenceModel::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 EquivalenceModel::ReadMaxBurnUp(const string &line) -{ - DBGL - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_maxburnup" ) // Check the keyword - { - ERROR << " Bad keyword : \"k_maxburnup\" not found !" << endl; - exit(1); - } - - fMaximalBU = atof(StringLine::NextWord(line, pos, ' ').c_str()); - - DBGL -} - -//________________________________________________________________________ -void EquivalenceModel::ReadTargetParameter(const string &line) -{ - DBGL - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_targetparameter" ) // Check the keyword - { - ERROR << " Bad keyword : \"k_targetparameter\" not found !" << endl; - exit(1); - } - - fTargetParameter = StringLine::NextWord(line, pos, ' '); - - DBGL -} - -//________________________________________________________________________ -void EquivalenceModel::ReadPredictorType(const string &line) -{ - DBGL - - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_predictortype" ) // Check the keyword - { - ERROR << " Bad keyword : \"k_predictortype\" not found !" << endl; - exit(1); - } - - fPredictorType = StringLine::NextWord(line, pos, ' '); - - DBGL -} -//________________________________________________________________________ -void EquivalenceModel::ReadOutput(const string &line) -{ - DBGL - - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_output" ) // Check the keyword - { - ERROR << " Bad keyword : \"k_output\" not found !" << endl; - exit(1); - } - - fOutput = StringLine::NextWord(line, pos, ' '); - - DBGL -} -//________________________________________________________________________ -void EquivalenceModel::ReadBuffer(const string &line) -{ - DBGL - - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_buffer" ) // Check the keyword - { - ERROR << " Bad keyword : \"k_buffer\" not found !" << endl; - exit(1); - } - - fBuffer = StringLine::NextWord(line, pos, ' '); - - DBGL -} -//________________________________________________________________________ -void EquivalenceModel::ReadModelParameter(const string &line) -{ - DBGL - - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_modelparameter" ) // Check the keyword - { - ERROR << " Bad keyword : \"k_modelparameter\" not found !" << endl; - exit(1); - } - - keyword = StringLine::NextWord(line, pos, ' '); - - fModelParameter[keyword] = -1; - - DBGL -} - -//________________________________________________________________________ -void EquivalenceModel::ReadTargetParameterStDev(const string &line) -{ - DBGL - - int pos = 0; - string keyword = tlc(StringLine::NextWord(line, pos, ' ')); - if( keyword != "k_targetparameterstdev" ) // Check the keyword - { - ERROR << " Bad keyword : \"k_targetparameterstdev\" not found !" << endl; - exit(1); - } - - fTargetParameterStDev = atof(StringLine::NextWord(line, pos, ' ').c_str());; - - DBGL -} - -//________________________________________________________________________ -void EquivalenceModel::PrintInfo() -{ - INFO << "Reactor Type : "<< fDBRType << endl; - INFO << "Fuel Type : "<< fDBFType << endl; - INFO << "Specific Power [W/g]: "<< fSpecificPower << endl; - - map < string , IsotopicVector >::iterator it_s_IV; - map < string , double >::iterator it_s_D; - - for( it_s_IV = fStreamList.begin(); it_s_IV != fStreamList.end(); it_s_IV++) - { - INFO <<(* it_s_IV).first<<" (Z A I) :" << endl; - map<ZAI ,double >::iterator it1; - map<ZAI ,double > fMap1 = fStreamList[(* it_s_IV).first].GetIsotopicQuantity(); - for(it1 = fMap1.begin() ; it1 != fMap1.end() ; it1++) - INFO << (*it1).first.Z() <<" "<< (*it1).first.A() <<" "<< (*it1).first.I() << endl; - } - INFO<<"Minimum fraction in the fuel for each material : "<<endl; - for( it_s_D = fStreamListEqMMassFractionMin.begin(); it_s_D != fStreamListEqMMassFractionMin.end(); it_s_D++) - { - INFO <<(* it_s_D).first<<" "<<fStreamListEqMMassFractionMin[(* it_s_D).first]<<endl; - } - - INFO<<"Maximum fraction in the fuel for each material : "<<endl; - for( it_s_D = fStreamListEqMMassFractionMax.begin(); it_s_D != fStreamListEqMMassFractionMax.end(); it_s_D++) - { - INFO <<(* it_s_D).first<<" "<<fStreamListEqMMassFractionMax[(* it_s_D).first]<<endl; - } - + DBGL + bool IsInDomain = true; - INFO<<"ZAI limits (validity domain)[prop in fresh fuel] (Z A I min max) :"<<endl; - for (map< ZAI,pair<double,double> >::iterator Domain_it = fZAILimits.begin(); Domain_it != fZAILimits.end(); Domain_it++) - { - double ThatZAIMin = Domain_it->second.first; - double ThatZAIMax = Domain_it->second.second; - int Z = Domain_it->first.Z(); - int A = Domain_it->first.A(); - int I = Domain_it->first.I(); + if (fZAILimits.empty()) + { + WARNING << "Fresh Fuel variation domain is not set" << endl; + WARNING << "CLASS has no clue if the computed evolution for this fresh fuel is correct" << endl; + WARNING << "Proceed finger crossed !!" << endl; + return true; + } - INFO <<ThatZAIMin<<" < ZAI ("<< Z<< " " << A << " " << I<<")"<<" < "<<ThatZAIMax<< endl; - } + else + { + IsotopicVector IVNorm = IV / IV.GetSumOfAll(); + for (map< ZAI, pair<double, double> >::iterator Domain_it = fZAILimits.begin(); Domain_it != fZAILimits.end(); Domain_it++) + { + double ThatZAIProp = IVNorm.GetIsotopicQuantity()[Domain_it->first]; + double ThatZAIMin   = Domain_it->second.first; + double ThatZAIMax   = Domain_it->second.second; + if ( (ThatZAIProp > ThatZAIMax) || (ThatZAIProp <   ThatZAIMin) ) + { + IsInDomain = false; + + WARNING << "Fresh fuel out of model range" << endl; + WARNING << "\t AT LEAST this ZAI is accused to be outrange :" << endl; + WARNING << "\t\t" << Domain_it->first.Z() << " " << Domain_it->first.A() << " " << Domain_it->first.I() << endl; + WARNING << "\t\t min = " << ThatZAIMin   << " value = " << ThatZAIProp << " max = " << ThatZAIMax << endl; + WARNING << "\t IV accused :" << endl << endl; + WARNING << IVNorm.sPrint() << endl; + break; + } + } + } + DBGL + return IsInDomain; } diff --git a/source/src/XSModel.cxx b/source/src/XSModel.cxx index 096b4d154cc6cad3308d4bc3a910a88a5dce60fe..b0c947c0a03b66f7dc2861148a4926e96a2ffe0b 100644 --- a/source/src/XSModel.cxx +++ b/source/src/XSModel.cxx @@ -22,7 +22,6 @@ XSModel::XSModel():CLASSObject() freaded = false; } - //________________________________________________________________________ XSModel::XSModel(CLASSLogger* log):CLASSObject(log) @@ -31,7 +30,11 @@ XSModel::XSModel(CLASSLogger* log):CLASSObject(log) XSModel::LoadKeyword(); freaded = false; } +//________________________________________________________________________ +XSModel::~XSModel() +{ +} //________________________________________________________________________ void XSModel::ReadNFO() { diff --git a/test/main_test.cxx b/test/main_test.cxx deleted file mode 100644 index b2ad43c03e9992d45f87deeac82b45c6af412430..0000000000000000000000000000000000000000 --- a/test/main_test.cxx +++ /dev/null @@ -1,10 +0,0 @@ -#include <gtest/gtest.h> -#include "test_iv.inl" - -int main(int argc,char * argv[]) { - ::testing::InitGoogleTest(&argc,argv); - return RUN_ALL_TESTS(); -} - -// COMPILATION -// g++ -o CLASS_test main_test.cxx -I $CLASS_include -L $CLASS_lib -lCLASSpkg `root-config --cflags` `root-config --libs` -fopenmp -lgomp -lTMVA -lgtest