/*
Interface to Clermond-Ferrand artificial neural network
<p>
The CFMlpANN belong to the class of Multilayer Perceptrons (MLP), which are
feed-forward networks according to the following propagation schema:<br>
<center>
<img vspace=10 src="gif/tmva_mlp.gif" align="bottom" alt="Schema for artificial neural network">
</center>
The input layer contains as many neurons as input variables used in the MVA.
The output layer contains two neurons for the signal and background
event classes. In between the input and output layers are a variable number
of <i>k</i> hidden layers with arbitrary numbers of neurons. (While the
structure of the input and output layers is determined by the problem, the
hidden layers can be configured by the user through the option string
of the method booking.) <br>
As indicated in the sketch, all neuron inputs to a layer are linear
combinations of the neuron output of the previous layer. The transfer
from input to output within a neuron is performed by means of an "activation
function". In general, the activation function of a neuron can be
zero (deactivated), one (linear), or non-linear. The above example uses
a sigmoid activation function. The transfer function of the output layer
is usually linear. As a consequence: an ANN without hidden layer should
give identical discrimination power as a linear discriminant analysis (Fisher).
In case of one hidden layer, the ANN computes a linear combination of
sigmoid. <br>
The learning method used by the CFMlpANN is only stochastic.
*/
// End_Html
#include <string>
#include <cstdlib>
#include <iostream>
#include "TMatrix.h"
#include "TObjString.h"
#include "Riostream.h"
#include "TMath.h"
#include "TMVA/ClassifierFactory.h"
#include "TMVA/MethodCFMlpANN.h"
#include "TMVA/MethodCFMlpANN_def.h"
#include "TMVA/Tools.h"
REGISTER_METHOD(CFMlpANN)
ClassImp(TMVA::MethodCFMlpANN)
namespace TMVA {
Int_t MethodCFMlpANN_nsel = 0;
}
TMVA::MethodCFMlpANN* TMVA::MethodCFMlpANN::fgThis = 0;
TMVA::MethodCFMlpANN::MethodCFMlpANN( const TString& jobName,
const TString& methodTitle,
DataSetInfo& theData,
const TString& theOption,
TDirectory* theTargetDir ) :
TMVA::MethodBase( jobName, Types::kCFMlpANN, methodTitle, theData, theOption, theTargetDir ),
fNodes(0),
fYNN(0)
{
}
TMVA::MethodCFMlpANN::MethodCFMlpANN( DataSetInfo& theData,
const TString& theWeightFile,
TDirectory* theTargetDir ):
TMVA::MethodBase( Types::kCFMlpANN, theData, theWeightFile, theTargetDir ),
fNodes(0),
fYNN(0)
{
}
Bool_t TMVA::MethodCFMlpANN::HasAnalysisType( Types::EAnalysisType type, UInt_t numberClasses, UInt_t )
{
if (type == Types::kClassification && numberClasses == 2) return kTRUE;
return kFALSE;
}
void TMVA::MethodCFMlpANN::DeclareOptions()
{
DeclareOptionRef( fNcycles =3000, "NCycles", "Number of training cycles" );
DeclareOptionRef( fLayerSpec="N,N-1", "HiddenLayers", "Specification of hidden layer architecture" );
}
void TMVA::MethodCFMlpANN::ProcessOptions()
{
fNodes = new Int_t[20];
fNlayers = 2;
Int_t currentHiddenLayer = 1;
TString layerSpec(fLayerSpec);
while(layerSpec.Length()>0) {
TString sToAdd = "";
if (layerSpec.First(',')<0) {
sToAdd = layerSpec;
layerSpec = "";
}
else {
sToAdd = layerSpec(0,layerSpec.First(','));
layerSpec = layerSpec(layerSpec.First(',')+1,layerSpec.Length());
}
Int_t nNodes = 0;
if (sToAdd.BeginsWith("N") || sToAdd.BeginsWith("n")) { sToAdd.Remove(0,1); nNodes = GetNvar(); }
nNodes += atoi(sToAdd);
fNodes[currentHiddenLayer++] = nNodes;
fNlayers++;
}
fNodes[0] = GetNvar();
fNodes[fNlayers-1] = 2;
if (IgnoreEventsWithNegWeightsInTraining()) {
Log() << kFATAL << "Mechanism to ignore events with negative weights in training not yet available for method: "
<< GetMethodTypeName()
<< " --> please remove \"IgnoreNegWeightsInTraining\" option from booking string."
<< Endl;
}
Log() << kINFO << "Use configuration (nodes per layer): in=";
for (Int_t i=0; i<fNlayers-1; i++) Log() << kINFO << fNodes[i] << ":";
Log() << kINFO << fNodes[fNlayers-1] << "=out" << Endl;
Log() << "Use " << fNcycles << " training cycles" << Endl;
Int_t nEvtTrain = Data()->GetNTrainingEvents();
if (nEvtTrain>0) {
fData = new TMatrix( nEvtTrain, GetNvar() );
fClass = new vector<Int_t>( nEvtTrain );
UInt_t ivar;
for (Int_t ievt=0; ievt<nEvtTrain; ievt++) {
const Event * ev = GetEvent(ievt);
(*fClass)[ievt] = ev->IsSignal() ? 1 : 2;
for (ivar=0; ivar<GetNvar(); ivar++) {
(*fData)( ievt, ivar ) = ev->GetValue(ivar);
}
}
}
}
void TMVA::MethodCFMlpANN::Init( void )
{
SetNormalised( kTRUE );
fgThis = this;
TMVA::MethodCFMlpANN_nsel = 0;
}
TMVA::MethodCFMlpANN::~MethodCFMlpANN( void )
{
delete fData;
delete fClass;
delete[] fNodes;
if (fYNN!=0) {
for (Int_t i=0; i<fNlayers; i++) delete[] fYNN[i];
delete[] fYNN;
fYNN=0;
}
}
void TMVA::MethodCFMlpANN::Train( void )
{
Double_t dumDat(0);
Int_t ntrain(Data()->GetNTrainingEvents());
Int_t ntest(0);
Int_t nvar(GetNvar());
Int_t nlayers(fNlayers);
Int_t *nodes = new Int_t[nlayers];
Int_t ncycles(fNcycles);
for (Int_t i=0; i<nlayers; i++) nodes[i] = fNodes[i];
if (fYNN != 0) {
for (Int_t i=0; i<fNlayers; i++) delete[] fYNN[i];
delete[] fYNN;
fYNN = 0;
}
fYNN = new Double_t*[nlayers];
for (Int_t layer=0; layer<nlayers; layer++)
fYNN[layer] = new Double_t[fNodes[layer]];
#ifndef R__WIN32
Train_nn( &dumDat, &dumDat, &ntrain, &ntest, &nvar, &nlayers, nodes, &ncycles );
#else
Log() << kWARNING << "<Train> sorry CFMlpANN does not run on Windows" << Endl;
#endif
delete [] nodes;
}
Double_t TMVA::MethodCFMlpANN::GetMvaValue( Double_t* err )
{
Bool_t isOK = kTRUE;
const Event* ev = GetEvent();
vector<Double_t> inputVec( GetNvar() );
for (UInt_t ivar=0; ivar<GetNvar(); ivar++) inputVec[ivar] = ev->GetValue(ivar);
Double_t myMVA = EvalANN( inputVec, isOK );
if (!isOK) Log() << kFATAL << "EvalANN returns (!isOK) for event " << Endl;
if (err != 0) *err = -1;
return myMVA;
}
Double_t TMVA::MethodCFMlpANN::EvalANN( vector<Double_t>& inVar, Bool_t& isOK )
{
Double_t* xeev = new Double_t[GetNvar()];
for (UInt_t ivar=0; ivar<GetNvar(); ivar++) xeev[ivar] = inVar[ivar];
isOK = kTRUE;
for (UInt_t jvar=0; jvar<GetNvar(); jvar++) {
if (fVarn_1.xmax[jvar] < xeev[jvar]) xeev[jvar] = fVarn_1.xmax[jvar];
if (fVarn_1.xmin[jvar] > xeev[jvar]) xeev[jvar] = fVarn_1.xmin[jvar];
if (fVarn_1.xmax[jvar] == fVarn_1.xmin[jvar]) {
isOK = kFALSE;
xeev[jvar] = 0;
}
else {
xeev[jvar] = xeev[jvar] - ((fVarn_1.xmax[jvar] + fVarn_1.xmin[jvar])/2);
xeev[jvar] = xeev[jvar] / ((fVarn_1.xmax[jvar] - fVarn_1.xmin[jvar])/2);
}
}
NN_ava( xeev );
Double_t retval = 0.5*(1.0 + fYNN[fParam_1.layerm-1][0]);
delete [] xeev;
return retval;
}
void TMVA::MethodCFMlpANN::NN_ava( Double_t* xeev )
{
for (Int_t ivar=0; ivar<fNeur_1.neuron[0]; ivar++) fYNN[0][ivar] = xeev[ivar];
for (Int_t layer=1; layer<fParam_1.layerm; layer++) {
for (Int_t j=1; j<=fNeur_1.neuron[layer]; j++) {
Double_t x = Ww_ref(fNeur_1.ww, layer+1,j);
for (Int_t k=1; k<=fNeur_1.neuron[layer-1]; k++) {
x += fYNN[layer-1][k-1]*W_ref(fNeur_1.w, layer+1, j, k);
}
fYNN[layer][j-1] = NN_fonc( layer, x );
}
}
}
Double_t TMVA::MethodCFMlpANN::NN_fonc( Int_t i, Double_t u ) const
{
Double_t f(0);
if (u/fDel_1.temp[i] > 170) f = +1;
else if (u/fDel_1.temp[i] < -170) f = -1;
else {
Double_t yy = TMath::Exp(-u/fDel_1.temp[i]);
f = (1 - yy)/(1 + yy);
}
return f;
}
void TMVA::MethodCFMlpANN::ReadWeightsFromStream( istream & istr )
{
TString var;
UInt_t nva(0), lclass(0);
istr >> nva >> lclass;
if (GetNvar() != nva)
Log() << kFATAL << "<ReadWeightsFromFile> mismatch in number of variables" << Endl;
if (lclass != 2)
Log() << kFATAL << "<ReadWeightsFromFile> mismatch in number of classes" << Endl;
if (istr.eof( ))
Log() << kFATAL << "<ReadWeightsFromStream> reached EOF prematurely " << Endl;
for (UInt_t ivar=0; ivar<GetNvar(); ivar++)
istr >> fVarn_1.xmax[ivar] >> fVarn_1.xmin[ivar];
istr >> fParam_1.layerm;
if (fYNN != 0) {
for (Int_t i=0; i<fNlayers; i++) delete[] fYNN[i];
delete[] fYNN;
fYNN = 0;
}
fYNN = new Double_t*[fParam_1.layerm];
for (Int_t layer=0; layer<fParam_1.layerm; layer++) {
istr >> fNeur_1.neuron[layer];
fYNN[layer] = new Double_t[fNeur_1.neuron[layer]];
}
const Int_t nchar( 100 );
char* dumchar = new char[nchar];
for (Int_t layer=1; layer<=fParam_1.layerm-1; layer++) {
Int_t nq = fNeur_1.neuron[layer]/10;
Int_t nr = fNeur_1.neuron[layer] - nq*10;
Int_t kk(0);
if (nr==0) kk = nq;
else kk = nq+1;
for (Int_t k=1; k<=kk; k++) {
Int_t jmin = 10*k - 9;
Int_t jmax = 10*k;
if (fNeur_1.neuron[layer]<jmax) jmax = fNeur_1.neuron[layer];
for (Int_t j=jmin; j<=jmax; j++) {
istr >> Ww_ref(fNeur_1.ww, layer+1, j);
}
for (Int_t i=1; i<=fNeur_1.neuron[layer-1]; i++) {
for (Int_t j=jmin; j<=jmax; j++) {
istr >> W_ref(fNeur_1.w, layer+1, j, i);
}
}
istr.getline( dumchar, nchar );
}
}
for (Int_t layer=0; layer<fParam_1.layerm; layer++) {
istr.getline( dumchar, nchar );
istr.getline( dumchar, nchar );
istr >> fDel_1.temp[layer];
}
if ((Int_t)GetNvar() != fNeur_1.neuron[0]) {
Log() << kFATAL << "<ReadWeightsFromFile> mismatch in zeroth layer:"
<< GetNvar() << " " << fNeur_1.neuron[0] << Endl;
}
fNlayers = fParam_1.layerm;
delete dumchar;
}
Int_t TMVA::MethodCFMlpANN::DataInterface( Double_t* , Double_t* ,
Int_t* , Int_t* ,
Int_t* , Int_t* nvar,
Double_t* xpg, Int_t* iclass, Int_t* ikend )
{
*ikend = 0;
TMVA::MethodCFMlpANN* opt = TMVA::MethodCFMlpANN::This();
if (0 == xpg) {
Log() << kFATAL << "ERROR in MethodCFMlpANN_DataInterface zero pointer xpg" << Endl;
}
if (*nvar != (Int_t)opt->GetNvar()) {
Log() << kFATAL << "ERROR in MethodCFMlpANN_DataInterface mismatch in num of variables: "
<< *nvar << " " << opt->GetNvar() << Endl;
}
*iclass = (int)opt->GetClass( TMVA::MethodCFMlpANN_nsel );
for (UInt_t ivar=0; ivar<opt->GetNvar(); ivar++)
xpg[ivar] = (double)opt->GetData( TMVA::MethodCFMlpANN_nsel, ivar );
++TMVA::MethodCFMlpANN_nsel;
return 0;
}
void TMVA::MethodCFMlpANN::AddWeightsXMLTo( void* parent ) const
{
void *wght = gTools().xmlengine().NewChild(parent, 0, "Weights");
gTools().AddAttr(wght,"NVars",fParam_1.nvar);
gTools().AddAttr(wght,"NClasses",fParam_1.lclass);
gTools().AddAttr(wght,"NLayers",fParam_1.layerm);
void* minmaxnode = gTools().xmlengine().NewChild(wght, 0, "VarMinMax");
stringstream s;
s.precision( 16 );
for (Int_t ivar=0; ivar<fParam_1.nvar; ivar++)
s << std::scientific << fVarn_1.xmin[ivar] << " " << fVarn_1.xmax[ivar] << " ";
gTools().xmlengine().AddRawLine( minmaxnode, s.str().c_str() );
void* neurons = gTools().xmlengine().NewChild(wght, 0, "NNeurons");
stringstream n;
n.precision( 16 );
for (Int_t layer=0; layer<fParam_1.layerm; layer++)
n << std::scientific << fNeur_1.neuron[layer] << " ";
gTools().xmlengine().AddRawLine( neurons, n.str().c_str() );
for (Int_t layer=1; layer<fParam_1.layerm; layer++) {
void* layernode = gTools().xmlengine().NewChild(wght, 0, "Layer"+gTools().StringFromInt(layer));
gTools().AddAttr(layernode,"NNeurons",fNeur_1.neuron[layer]);
void* neuronnode=NULL;
for (Int_t neuron=0; neuron<fNeur_1.neuron[layer]; neuron++) {
neuronnode = gTools().xmlengine().NewChild(layernode,0,"Neuron"+gTools().StringFromInt(neuron));
stringstream weights;
weights.precision( 16 );
weights << std::scientific << Ww_ref(fNeur_1.ww, layer+1, neuron+1);
for (Int_t i=0; i<fNeur_1.neuron[layer-1]; i++) {
weights << " " << std::scientific << W_ref(fNeur_1.w, layer+1, neuron+1, i+1);
}
gTools().xmlengine().AddRawLine( neuronnode, weights.str().c_str() );
}
}
void* tempnode = gTools().xmlengine().NewChild(wght, 0, "LayerTemp");
stringstream temp;
temp.precision( 16 );
for (Int_t layer=0; layer<fParam_1.layerm; layer++) {
temp << std::scientific << fDel_1.temp[layer] << " ";
}
gTools().xmlengine().AddRawLine(tempnode, temp.str().c_str() );
}
void TMVA::MethodCFMlpANN::ReadWeightsFromXML( void* wghtnode )
{
gTools().ReadAttr( wghtnode, "NLayers",fParam_1.layerm );
void* minmaxnode = gTools().xmlengine().GetChild(wghtnode);
const char* minmaxcontent = gTools().xmlengine().GetNodeContent(minmaxnode);
std::stringstream content(minmaxcontent);
for (UInt_t ivar=0; ivar<GetNvar(); ivar++)
content >> fVarn_1.xmin[ivar] >> fVarn_1.xmax[ivar];
if (fYNN != 0) {
for (Int_t i=0; i<fNlayers; i++) delete[] fYNN[i];
delete[] fYNN;
fYNN = 0;
}
fYNN = new Double_t*[fParam_1.layerm];
void *layernode=gTools().xmlengine().GetNext(minmaxnode);
const char* neuronscontent = gTools().xmlengine().GetNodeContent(layernode);
stringstream ncontent(neuronscontent);
for (Int_t layer=0; layer<fParam_1.layerm; layer++) {
ncontent >> fNeur_1.neuron[layer];
fYNN[layer] = new Double_t[fNeur_1.neuron[layer]];
}
for (Int_t layer=1; layer<fParam_1.layerm; layer++) {
layernode=gTools().xmlengine().GetNext(layernode);
void* neuronnode=NULL;
neuronnode = gTools().xmlengine().GetChild(layernode);
for (Int_t neuron=0; neuron<fNeur_1.neuron[layer]; neuron++) {
const char* neuronweights = gTools().xmlengine().GetNodeContent(neuronnode);
stringstream weights(neuronweights);
weights >> Ww_ref(fNeur_1.ww, layer+1, neuron+1);
for (Int_t i=0; i<fNeur_1.neuron[layer-1]; i++) {
weights >> W_ref(fNeur_1.w, layer+1, neuron+1, i+1);
}
neuronnode=gTools().xmlengine().GetNext(neuronnode);
}
}
void* tempnode=gTools().xmlengine().GetNext(layernode);
const char* temp = gTools().xmlengine().GetNodeContent(tempnode);
stringstream t(temp);
for (Int_t layer=0; layer<fParam_1.layerm; layer++) {
t >> fDel_1.temp[layer];
}
fNlayers = fParam_1.layerm;
}
void TMVA::MethodCFMlpANN::PrintWeights( std::ostream & o ) const
{
o << "Number of vars " << fParam_1.nvar << endl;
o << "Output nodes " << fParam_1.lclass << endl;
for (Int_t ivar=0; ivar<fParam_1.nvar; ivar++)
o << "Var " << ivar << " [" << fVarn_1.xmin[ivar] << " - " << fVarn_1.xmax[ivar] << "]" << endl;
o << "Number of layers " << fParam_1.layerm << endl;
o << "Nodes per layer ";
for (Int_t layer=0; layer<fParam_1.layerm; layer++)
o << fNeur_1.neuron[layer] << " ";
o << endl;
for (Int_t layer=1; layer<=fParam_1.layerm-1; layer++) {
Int_t nq = fNeur_1.neuron[layer]/10;
Int_t nr = fNeur_1.neuron[layer] - nq*10;
Int_t kk(0);
if (nr==0) kk = nq;
else kk = nq+1;
for (Int_t k=1; k<=kk; k++) {
Int_t jmin = 10*k - 9;
Int_t jmax = 10*k;
Int_t i, j;
if (fNeur_1.neuron[layer]<jmax) jmax = fNeur_1.neuron[layer];
for (j=jmin; j<=jmax; j++) {
o << Ww_ref(fNeur_1.ww, layer+1, j) << " ";
}
o << endl;
for (i=1; i<=fNeur_1.neuron[layer-1]; i++) {
for (j=jmin; j<=jmax; j++) {
o << W_ref(fNeur_1.w, layer+1, j, i) << " ";
}
o << endl;
}
o << endl;
}
}
for (Int_t layer=0; layer<fParam_1.layerm; layer++) {
o << "Del.temp in layer " << layer << " : " << fDel_1.temp[layer] << endl;
}
}
void TMVA::MethodCFMlpANN::MakeClassSpecific( std::ostream& fout, const TString& className ) const
{
fout << " // not implemented for class: \"" << className << "\"" << endl;
fout << "};" << endl;
}
void TMVA::MethodCFMlpANN::MakeClassSpecificHeader( std::ostream& , const TString& ) const
{
}
void TMVA::MethodCFMlpANN::GetHelpMessage() const
{
Log() << Endl;
Log() << gTools().Color("bold") << "--- Short description:" << gTools().Color("reset") << Endl;
Log() << Endl;
Log() << "<None>" << Endl;
Log() << Endl;
Log() << gTools().Color("bold") << "--- Performance optimisation:" << gTools().Color("reset") << Endl;
Log() << Endl;
Log() << "<None>" << Endl;
Log() << Endl;
Log() << gTools().Color("bold") << "--- Performance tuning via configuration options:" << gTools().Color("reset") << Endl;
Log() << Endl;
Log() << "<None>" << Endl;
}