doc/hackathon/TMVA__CNN__Classification_8C.html

/***


    # TMVA Classification Example Using a Convolutional Neural Network


**/


/// Helper function to create input images data

/// we create a signal and background 2D histograms from 2d gaussians

/// with a location (means in X and Y)  different for each event

/// The difference between signal and background is in the gaussian width.

/// The width for the background gaussian is slightly larger than the signal width by few % values

///

///

void MakeImagesTree(int n, int nh, int nw)

{


   // image size (nh x nw)

   const int ntot = nh * nw;

   const TString fileOutName = TString::Format("images_data_%dx%d.root", nh, nw);

   TFile f(fileOutName, "RECREATE");


   const int nRndmEvts = 10000; // number of events we use to fill each image

   double delta_sigma = 0.1;    // 5% difference in the sigma

   double pixelNoise = 5;


   double sX1 = 3;

   double sY1 = 3;

   double sX2 = sX1 + delta_sigma;

   double sY2 = sY1 - delta_sigma;


   TH2D h1("h1", "h1", nh, 0, 10, nw, 0, 10);

   TH2D h2("h2", "h2", nh, 0, 10, nw, 0, 10);


   TF2 f1("f1", "xygaus");

   TF2 f2("f2", "xygaus");


   TTree sgn("sig_tree", "signal_tree");

   TTree bkg("bkg_tree", "background_tree");


   std::vector<float> x1(ntot);

   std::vector<float> x2(ntot);


   // create signal and background trees with a single branch

   // an std::vector<float> of size nh x nw containing the image data


   std::vector<float> *px1 = &x1;

   std::vector<float> *px2 = &x2;


   bkg.Branch("vars", "std::vector<float>", &px1);

   sgn.Branch("vars", "std::vector<float>", &px2);


   // std::cout << "create tree " << std::endl;


   sgn.SetDirectory(&f);

   bkg.SetDirectory(&f);


   f1.SetParameters(1, 5, sX1, 5, sY1);

   f2.SetParameters(1, 5, sX2, 5, sY2);

   gRandom->SetSeed(0);

   std::cout << "Filling ROOT tree " << std::endl;

   for (int i = 0; i < n; ++i) {

      if (i % 1000 == 0)

         std::cout << "Generating image event ... " << i << std::endl;

      h1.Reset();

      h2.Reset();

      // generate random means in range [3,7] to be not too much on the border

      f1.SetParameter(1, gRandom->Uniform(3, 7));

      f1.SetParameter(3, gRandom->Uniform(3, 7));

      f2.SetParameter(1, gRandom->Uniform(3, 7));

      f2.SetParameter(3, gRandom->Uniform(3, 7));


      h1.FillRandom("f1", nRndmEvts);

      h2.FillRandom("f2", nRndmEvts);


      for (int k = 0; k < nh; ++k) {

         for (int l = 0; l < nw; ++l) {

            int m = k * nw + l;

            // add some noise in each bin

            x1[m] = h1.GetBinContent(k + 1, l + 1) + gRandom->Gaus(0, pixelNoise);

            x2[m] = h2.GetBinContent(k + 1, l + 1) + gRandom->Gaus(0, pixelNoise);

         }

      }

      sgn.Fill();

      bkg.Fill();

   }

   sgn.Write();

   bkg.Write();


   Info("MakeImagesTree", "Signal and background tree with images data written to the file %s", f.GetName());

   sgn.Print();

   bkg.Print();

   f.Close();

}


/// @brief Run the TMVA CNN Classification example

/// @param nevts : number of signal/background events. Use by default a low value (1000)

///                but increase to at least 5000 to get a good result

/// @param opt :   vector of bool with method used (default all on if available). The order is:

///                   - TMVA CNN

///                   - Keras CNN

///                   - TMVA DNN

///                   - TMVA BDT

///                   - PyTorch CNN

void TMVA_CNN_Classification(int nevts = 1000, std::vector<bool> opt = {1, 1, 1, 1, 1})

{


   int imgSize = 16 * 16;

   TString inputFileName = "images_data_16x16.root";


   bool fileExist = !gSystem->AccessPathName(inputFileName);


   // if file does not exists create it

   if (!fileExist) {

      MakeImagesTree(nevts, 16, 16);

   }


   bool useTMVACNN = (opt.size() > 0) ? opt[0] : false;

   bool useKerasCNN = (opt.size() > 1) ? opt[1] : false;

   bool useTMVADNN = (opt.size() > 2) ? opt[2] : false;

   bool useTMVABDT = (opt.size() > 3) ? opt[3] : false;

   bool usePyTorchCNN = (opt.size() > 4) ? opt[4] : false;

#ifndef R__HAS_TMVACPU

#ifndef R__HAS_TMVAGPU

   Warning("TMVA_CNN_Classification",

           "TMVA is not build with GPU or CPU multi-thread support. Cannot use TMVA Deep Learning for CNN");

   useTMVACNN = false;

#endif

#endif


   bool writeOutputFile = true;


#ifdef R__USE_IMT

   int num_threads = 4;  // use by default 4 threads if value is not set before

   // switch off MT in OpenBLAS to avoid conflict with tbb

   gSystem->Setenv("OMP_NUM_THREADS", "1");


   // do enable MT running

   if (num_threads >= 0) {

      ROOT::EnableImplicitMT(num_threads);

   }

#endif


   TMVA::Tools::Instance();


   std::cout << "Running with nthreads  = " << ROOT::GetThreadPoolSize() << std::endl;


#ifdef R__HAS_PYMVA

   gSystem->Setenv("KERAS_BACKEND", "tensorflow");

   // for using Keras

   TMVA::PyMethodBase::PyInitialize();

#else

   useKerasCNN = false;

   usePyTorchCNN = false;

#endif


   TFile *outputFile = nullptr;

   if (writeOutputFile)

      outputFile = TFile::Open("TMVA_CNN_ClassificationOutput.root", "RECREATE");


   /***

       ## Create TMVA Factory


    Create the Factory class. Later you can choose the methods

    whose performance you'd like to investigate.


    The factory is the major TMVA object you have to interact with. Here is the list of parameters you need to pass


    - The first argument is the base of the name of all the output

    weight files in the directory weight/ that will be created with the

    method parameters


    - The second argument is the output file for the training results


    - The third argument is a string option defining some general configuration for the TMVA session.

      For example all TMVA output can be suppressed by removing the "!" (not) in front of the "Silent" argument in the

   option string


    - note that we disable any pre-transformation of the input variables and we avoid computing correlations between

   input variables

   ***/


   TMVA::Factory factory(

      "TMVA_CNN_Classification", outputFile,

      "!V:ROC:!Silent:Color:AnalysisType=Classification:Transformations=None:!Correlations");


   /***


       ## Declare DataLoader(s)


       The next step is to declare the DataLoader class that deals with input variables


       Define the input variables that shall be used for the MVA training

       note that you may also use variable expressions, which can be parsed by TTree::Draw( "expression" )]


       In this case the input data consists of an image of 16x16 pixels. Each single pixel is a branch in a ROOT TTree


   **/


   TMVA::DataLoader loader("dataset");


   /***


       ## Setup Dataset(s)


       Define input data file and signal and background trees


   **/


   std::unique_ptr<TFile> inputFile{TFile::Open(inputFileName)};

   if (!inputFile) {

      Error("TMVA_CNN_Classification", "Error opening input file %s - exit", inputFileName.Data());

      return;

   }


   // --- Register the training and test trees


   auto signalTree = inputFile->Get<TTree>("sig_tree");

   auto backgroundTree = inputFile->Get<TTree>("bkg_tree");


   if (!signalTree) {

      Error("TMVA_CNN_Classification", "Could not find signal tree in file '%s'", inputFileName.Data());

      return;

   }

   if (!backgroundTree) {

      Error("TMVA_CNN_Classification", "Could not find background tree in file '%s'", inputFileName.Data());

      return;

   }


   int nEventsSig = signalTree->GetEntries();

   int nEventsBkg = backgroundTree->GetEntries();


   // global event weights per tree (see below for setting event-wise weights)

   Double_t signalWeight = 1.0;

   Double_t backgroundWeight = 1.0;


   // You can add an arbitrary number of signal or background trees

   loader.AddSignalTree(signalTree, signalWeight);

   loader.AddBackgroundTree(backgroundTree, backgroundWeight);


   /// add event variables (image)

   /// use new method (from ROOT 6.20 to add a variable array for all image data)

   loader.AddVariablesArray("vars", imgSize);


   // Set individual event weights (the variables must exist in the original TTree)

   //    for signal    : factory->SetSignalWeightExpression    ("weight1*weight2");

   //    for background: factory->SetBackgroundWeightExpression("weight1*weight2");

   // loader.SetBackgroundWeightExpression( "weight" );


   // Apply additional cuts on the signal and background samples (can be different)

   TCut mycuts = ""; // for example: TCut mycuts = "abs(var1)<0.5 && abs(var2-0.5)<1";

   TCut mycutb = ""; // for example: TCut mycutb = "abs(var1)<0.5";


   // Tell the factory how to use the training and testing events

   //

   // If no numbers of events are given, half of the events in the tree are used

   // for training, and the other half for testing:

   //    loader.PrepareTrainingAndTestTree( mycut, "SplitMode=random:!V" );

   // It is possible also to specify the number of training and testing events,

   // note we disable the computation of the correlation matrix of the input variables


   int nTrainSig = 0.8 * nEventsSig;

   int nTrainBkg = 0.8 * nEventsBkg;


   // build the string options for DataLoader::PrepareTrainingAndTestTree

   TString prepareOptions = TString::Format(

      "nTrain_Signal=%d:nTrain_Background=%d:SplitMode=Random:SplitSeed=100:NormMode=NumEvents:!V:!CalcCorrelations",

      nTrainSig, nTrainBkg);


   loader.PrepareTrainingAndTestTree(mycuts, mycutb, prepareOptions);


   /***


       DataSetInfo              : [dataset] : Added class "Signal"

       : Add Tree sig_tree of type Signal with 10000 events

       DataSetInfo              : [dataset] : Added class "Background"

       : Add Tree bkg_tree of type Background with 10000 events


   **/


   /****

        # Booking Methods


        Here we book the TMVA methods. We book a Boosted Decision Tree method (BDT)


   **/


   // Boosted Decision Trees

   if (useTMVABDT) {

      factory.BookMethod(&loader, TMVA::Types::kBDT, "BDT",

                         "!V:NTrees=200:MinNodeSize=2.5%:MaxDepth=2:BoostType=AdaBoost:AdaBoostBeta=0.5:"

                         "UseBaggedBoost:BaggedSampleFraction=0.5:SeparationType=GiniIndex:nCuts=20");

   }

   /**


      #### Booking Deep Neural Network


      Here we book the DNN of TMVA. See the example TMVA_Higgs_Classification.C for a detailed description of the

   options


   **/


   if (useTMVADNN) {


      TString layoutString(

         "Layout=DENSE|100|RELU,BNORM,DENSE|100|RELU,BNORM,DENSE|100|RELU,BNORM,DENSE|100|RELU,DENSE|1|LINEAR");


      // Training strategies

      // one can catenate several training strings with different parameters (e.g. learning rates or regularizations

      // parameters) The training string must be concatenates with the `|` delimiter

      TString trainingString1("LearningRate=1e-3,Momentum=0.9,Repetitions=1,"

                              "ConvergenceSteps=5,BatchSize=100,TestRepetitions=1,"

                              "MaxEpochs=10,WeightDecay=1e-4,Regularization=None,"

                              "Optimizer=ADAM,DropConfig=0.0+0.0+0.0+0.");


      TString trainingStrategyString("TrainingStrategy=");

      trainingStrategyString += trainingString1; // + "|" + trainingString2 + ....


      // Build now the full DNN Option string


      TString dnnOptions("!H:V:ErrorStrategy=CROSSENTROPY:VarTransform=None:"

                         "WeightInitialization=XAVIER");

      dnnOptions.Append(":");

      dnnOptions.Append(layoutString);

      dnnOptions.Append(":");

      dnnOptions.Append(trainingStrategyString);


      TString dnnMethodName = "TMVA_DNN_CPU";

// use GPU if available

#ifdef R__HAS_TMVAGPU

      dnnOptions += ":Architecture=GPU";

      dnnMethodName = "TMVA_DNN_GPU";

#elif defined(R__HAS_TMVACPU)

      dnnOptions += ":Architecture=CPU";

#endif


      factory.BookMethod(&loader, TMVA::Types::kDL, dnnMethodName, dnnOptions);

   }


   /***

    ### Book Convolutional Neural Network in TMVA


    For building a CNN one needs to define


    -  Input Layout :  number of channels (in this case = 1)  | image height | image width

    -  Batch Layout :  batch size | number of channels | image size = (height*width)


    Then one add Convolutional layers and MaxPool layers.


    -  For Convolutional layer the option string has to be:

       - CONV | number of units | filter height | filter width | stride height | stride width | padding height | paddig

   width | activation function


       - note in this case we are using a filer 3x3 and padding=1 and stride=1 so we get the output dimension of the

   conv layer equal to the input


      - note we use after the first convolutional layer a batch normalization layer. This seems to help significantly the

   convergence


     - For the MaxPool layer:

        - MAXPOOL  | pool height | pool width | stride height | stride width


    The RESHAPE layer is needed to flatten the output before the Dense layer


    Note that to run the CNN is required to have CPU  or GPU support


   ***/


   if (useTMVACNN) {


      TString inputLayoutString("InputLayout=1|16|16");


      // Batch Layout

      TString layoutString("Layout=CONV|10|3|3|1|1|1|1|RELU,BNORM,CONV|10|3|3|1|1|1|1|RELU,MAXPOOL|2|2|1|1,"

                           "RESHAPE|FLAT,DENSE|100|RELU,DENSE|1|LINEAR");


      // Training strategies.

      TString trainingString1("LearningRate=1e-3,Momentum=0.9,Repetitions=1,"

                              "ConvergenceSteps=5,BatchSize=100,TestRepetitions=1,"

                              "MaxEpochs=10,WeightDecay=1e-4,Regularization=None,"

                              "Optimizer=ADAM,DropConfig=0.0+0.0+0.0+0.0");


      TString trainingStrategyString("TrainingStrategy=");

      trainingStrategyString +=

         trainingString1; // + "|" + trainingString2 + "|" + trainingString3; for concatenating more training strings


      // Build full CNN Options.

      TString cnnOptions("!H:V:ErrorStrategy=CROSSENTROPY:VarTransform=None:"

                         "WeightInitialization=XAVIER");


      cnnOptions.Append(":");

      cnnOptions.Append(inputLayoutString);

      cnnOptions.Append(":");

      cnnOptions.Append(layoutString);

      cnnOptions.Append(":");

      cnnOptions.Append(trainingStrategyString);


      //// New DL (CNN)

      TString cnnMethodName = "TMVA_CNN_CPU";

// use GPU if available

#ifdef R__HAS_TMVAGPU

      cnnOptions += ":Architecture=GPU";

      cnnMethodName = "TMVA_CNN_GPU";

#else

      cnnOptions += ":Architecture=CPU";

      cnnMethodName = "TMVA_CNN_CPU";

#endif


      factory.BookMethod(&loader, TMVA::Types::kDL, cnnMethodName, cnnOptions);

   }


   /**

      ### Book Convolutional Neural Network in Keras using a generated model


   **/


#ifdef R__HAS_PYMVA

   // The next section uses Python packages, execute it only if PyMVA is available

   TString tmva_python_exe{TMVA::Python_Executable()};

   TString python_exe = tmva_python_exe.IsNull() ? "python" : tmva_python_exe;


   if (useKerasCNN) {


      Info("TMVA_CNN_Classification", "Building convolutional keras model");

      // create python script which can be executed

      // create 2 conv2d layer + maxpool + dense

      TMacro m;

      m.AddLine("import tensorflow");

      m.AddLine("from tensorflow.keras.models import Sequential");

      m.AddLine("from tensorflow.keras.optimizers import Adam");

      m.AddLine(

         "from tensorflow.keras.layers import Input, Dense, Dropout, Flatten, Conv2D, MaxPooling2D, Reshape, BatchNormalization");

      m.AddLine("");

      m.AddLine("model = Sequential() ");

      m.AddLine("model.add(Reshape((16, 16, 1), input_shape = (256, )))");

      m.AddLine("model.add(Conv2D(10, kernel_size = (3, 3), kernel_initializer = 'glorot_normal',activation = "

                "'relu', padding = 'same'))");

      m.AddLine("model.add(BatchNormalization())");

      m.AddLine("model.add(Conv2D(10, kernel_size = (3, 3), kernel_initializer = 'glorot_normal',activation = "

                "'relu', padding = 'same'))");

      // m.AddLine("model.add(BatchNormalization())");

      m.AddLine("model.add(MaxPooling2D(pool_size = (2, 2), strides = (1,1))) ");

      m.AddLine("model.add(Flatten())");

      m.AddLine("model.add(Dense(256, activation = 'relu')) ");

      m.AddLine("model.add(Dense(2, activation = 'sigmoid')) ");

      m.AddLine("model.compile(loss = 'binary_crossentropy', optimizer = Adam(learning_rate = 0.001), weighted_metrics = ['accuracy'])");

      m.AddLine("model.save('model_cnn.keras')");

      m.AddLine("model.summary()");


      m.SaveSource("make_cnn_model.py");

      // execute

      gSystem->Exec(python_exe + " make_cnn_model.py");


      if (gSystem->AccessPathName("model_cnn.keras")) {

         Warning("TMVA_CNN_Classification", "Error creating Keras model file - skip using Keras");

      } else {

         // book PyKeras method only if Keras model could be created

         Info("TMVA_CNN_Classification", "Booking tf.Keras CNN model");

         factory.BookMethod(

            &loader, TMVA::Types::kPyKeras, "PyKeras",

            "H:!V:VarTransform=None:FilenameModel=model_cnn.keras:tf.keras:"

            "FilenameTrainedModel=trained_model_cnn.keras:NumEpochs=10:BatchSize=100:");

      }

   }


   if (usePyTorchCNN) {


      Info("TMVA_CNN_Classification", "Using Convolutional PyTorch Model");

      TString pyTorchFileName = gROOT->GetTutorialDir() + TString("/machine_learning/PyTorch_Generate_CNN_Model.py");

      // check that pytorch can be imported and file defining the model and used later when booking the method is

      // existing

      if (gSystem->Exec(python_exe + " -c 'import torch'") || gSystem->AccessPathName(pyTorchFileName)) {

         Warning("TMVA_CNN_Classification", "PyTorch is not installed or model building file is not existing - skip using PyTorch");

      } else {

         // book PyTorch method only if PyTorch model could be created

         Info("TMVA_CNN_Classification", "Booking PyTorch CNN model");

         TString methodOpt = "H:!V:VarTransform=None:FilenameModel=PyTorchModelCNN.pt:"

                             "FilenameTrainedModel=PyTorchTrainedModelCNN.pt:NumEpochs=10:BatchSize=100";

         methodOpt += TString(":UserCode=") + pyTorchFileName;

         factory.BookMethod(&loader, TMVA::Types::kPyTorch, "PyTorch", methodOpt);

      }

   }

#endif


   ////  ## Train Methods


   factory.TrainAllMethods();


   /// ## Test and Evaluate Methods


   factory.TestAllMethods();


   factory.EvaluateAllMethods();


   /// ## Plot ROC Curve


   auto c1 = factory.GetROCCurve(&loader);

   c1->Draw();


   // close outputfile to save output file

   outputFile->Close();

}

f
#define f(i)
Definition RSha256.hxx:104

Double_t
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73

Error
Error("WriteTObject","The current directory (%s) is not associated with a file. The object (%s) has not been written.", GetName(), objname)

Info
void Info(const char *location, const char *msgfmt,...)
Use this function for informational messages.
Definition TError.cxx:241

Warning
void Warning(const char *location, const char *msgfmt,...)
Use this function in warning situations.
Definition TError.cxx:252

x2
Option_t Option_t TPoint TPoint const char x2
Definition TGWin32VirtualXProxy.cxx:70

x1
Option_t Option_t TPoint TPoint const char x1
Definition TGWin32VirtualXProxy.cxx:70

gROOT
#define gROOT
Definition TROOT.h:417

gRandom
externTRandom * gRandom
Definition TRandom.h:62

gSystem
externTSystem * gSystem
Definition TSystem.h:582

TCut
A specialized string object used for TTree selections.
Definition TCut.h:25

TF2
Definition TF2.h:29

TFile
Definition TFile.h:130

TFile::Open
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)

TFile::Close
void Close(Option_t *option="") override
Delete all objects from memory and directory structure itself.

TH2D
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:400

TMVA::DataLoader
Definition DataLoader.h:50

TMVA::Factory
This is the main MVA steering class.
Definition Factory.h:80

TMVA::PyMethodBase::PyInitialize
static void PyInitialize()
Initialize Python interpreter.
Definition PyMethodBase.cxx:152

TMVA::Tools::Instance
static Tools & Instance()
Definition Tools.cxx:72

TMVA::Types::kPyKeras
@ kPyKeras
Definition Types.h:103

TMVA::Types::kBDT
@ kBDT
Definition Types.h:86

TMVA::Types::kDL
@ kDL
Definition Types.h:99

TMVA::Types::kPyTorch
@ kPyTorch
Definition Types.h:104

TMacro
Class supporting a collection of lines with C++ code.
Definition TMacro.h:31

TString
Basic string class.
Definition TString.h:138

TString::Data
const char * Data() const
Definition TString.h:384

TString::IsNull
Bool_t IsNull() const
Definition TString.h:422

TString::Format
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2385

TTree
A TTree represents a columnar dataset.
Definition TTree.h:89

TTree::GetEntries
virtual Long64_t GetEntries() const
Definition TTree.h:510

c1
return c1
Definition legend1.C:41

n
const Int_t n
Definition legend1.C:16

h1
TH1F * h1
Definition legend1.C:5

f1
TF1 * f1
Definition legend1.C:11

ROOT::EnableImplicitMT
void EnableImplicitMT(UInt_t numthreads=0)
Enable ROOT's implicit multi-threading for all objects and methods that provide an internal paralleli...
Definition TROOT.cxx:613

ROOT::GetThreadPoolSize
UInt_t GetThreadPoolSize()
Returns the size of ROOT's thread pool.
Definition TROOT.cxx:676

TMVA_CNN_Classification
Definition TMVA_CNN_Classification.py:1

TMVA::Python_Executable
TString Python_Executable()
Function to find current Python executable used by ROOT If "Python3" is installed,...
Definition PyMethodBase.cxx:43

m
TMarker m
Definition textangle.C:8

l
TLine l
Definition textangle.C:4