Logo ROOT   6.21/01
Reference Guide
TMultiLayerPerceptron.cxx
Go to the documentation of this file.
1 // @(#)root/mlp:$Id$
2 // Author: Christophe.Delaere@cern.ch 20/07/03
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2003, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /** \class TMultiLayerPerceptron
13 
14 
15 This class describes a neural network.
16 There are facilities to train the network and use the output.
17 
18 The input layer is made of inactive neurons (returning the
19 optionally normalized input) and output neurons are linear.
20 The type of hidden neurons is free, the default being sigmoids.
21 (One should still try to pass normalized inputs, e.g. between [0.,1])
22 
23 The basic input is a TTree and two (training and test) TEventLists.
24 Input and output neurons are assigned a value computed for each event
25 with the same possibilities as for TTree::Draw().
26 Events may be weighted individually or via TTree::SetWeight().
27 6 learning methods are available: kStochastic, kBatch,
28 kSteepestDescent, kRibierePolak, kFletcherReeves and kBFGS.
29 
30 This implementation, written by C. Delaere, is *inspired* from
31 the mlpfit package from J.Schwindling et al. with some extensions:
32 
33  - the algorithms are globally the same
34  - in TMultilayerPerceptron, there is no limitation on the number of
35  layers/neurons, while MLPFIT was limited to 2 hidden layers
36  - TMultilayerPerceptron allows you to save the network in a root file, and
37  provides more export functionalities
38  - TMultilayerPerceptron gives more flexibility regarding the normalization of
39  inputs/outputs
40  - TMultilayerPerceptron provides, thanks to Andrea Bocci, the possibility to
41  use cross-entropy errors, which allows to train a network for pattern
42  classification based on Bayesian posterior probability.
43 
44 ### Introduction
45 
46 Neural Networks are more and more used in various fields for data
47 analysis and classification, both for research and commercial
48 institutions. Some randomly chosen examples are:
49 
50  - image analysis
51  - financial movements predictions and analysis
52  - sales forecast and product shipping optimisation
53  - in particles physics: mainly for classification tasks (signal
54  over background discrimination)
55 
56 More than 50% of neural networks are multilayer perceptrons. This
57 implementation of multilayer perceptrons is inspired from the
58 <A HREF="http://schwind.home.cern.ch/schwind/MLPfit.html">MLPfit
59 package</A> originally written by Jerome Schwindling. MLPfit remains
60 one of the fastest tool for neural networks studies, and this ROOT
61 add-on will not try to compete on that. A clear and flexible Object
62 Oriented implementation has been chosen over a faster but more
63 difficult to maintain code. Nevertheless, the time penalty does not
64 exceed a factor 2.
65 
66 ### The MLP
67 
68 The multilayer perceptron is a simple feed-forward network with
69 the following structure:
70 
71 \image html mlp.png
72 
73 It is made of neurons characterized by a bias and weighted links
74 between them (let's call those links synapses). The input neurons
75 receive the inputs, normalize them and forward them to the first
76 hidden layer.
77 
78 Each neuron in any subsequent layer first computes a linear
79 combination of the outputs of the previous layer. The output of the
80 neuron is then function of that combination with <I>f</I> being
81 linear for output neurons or a sigmoid for hidden layers. This is
82 useful because of two theorems:
83 
84  1. A linear combination of sigmoids can approximate any
85  continuous function.
86  2. Trained with output = 1 for the signal and 0 for the
87  background, the approximated function of inputs X is the probability
88  of signal, knowing X.
89 
90 ### Learning methods
91 
92 The aim of all learning methods is to minimize the total error on
93 a set of weighted examples. The error is defined as the sum in
94 quadrature, divided by two, of the error on each individual output
95 neuron.
96 In all methods implemented, one needs to compute
97 the first derivative of that error with respect to the weights.
98 Exploiting the well-known properties of the derivative, especially the
99 derivative of compound functions, one can write:
100 
101  - for a neuron: product of the local derivative with the
102  weighted sum on the outputs of the derivatives.
103  - for a synapse: product of the input with the local derivative
104  of the output neuron.
105 
106 This computation is called back-propagation of the errors. A
107 loop over all examples is called an epoch.
108 Six learning methods are implemented.
109 
110 #### Stochastic minimization:
111 
112 is the most trivial learning method. This is the Robbins-Monro
113 stochastic approximation applied to multilayer perceptrons. The
114 weights are updated after each example according to the formula:
115 \f$w_{ij}(t+1) = w_{ij}(t) + \Delta w_{ij}(t)\f$
116 
117 with
118 
119 \f$\Delta w_{ij}(t) = - \eta(d e_p / d w_{ij} + \delta) + \epsilon \Delta w_{ij}(t-1)\f$
120 
121 The parameters for this method are Eta, EtaDecay, Delta and
122 Epsilon.
123 
124 #### Steepest descent with fixed step size (batch learning):
125 
126 It is the same as the stochastic
127 minimization, but the weights are updated after considering all the
128 examples, with the total derivative dEdw. The parameters for this
129 method are Eta, EtaDecay, Delta and Epsilon.
130 
131 #### Steepest descent algorithm:
132 
133 Weights are set to the minimum along the line defined by the gradient. The
134 only parameter for this method is Tau. Lower tau = higher precision =
135 slower search. A value Tau = 3 seems reasonable.
136 
137 #### Conjugate gradients with the Polak-Ribiere updating formula:
138 
139 Weights are set to the minimum along the line defined by the conjugate gradient.
140 Parameters are Tau and Reset, which defines the epochs where the direction is
141 reset to the steepest descent.
142 
143 #### Conjugate gradients with the Fletcher-Reeves updating formula:
144 
145 Weights are set to the minimum along the line defined by the conjugate gradient. Parameters
146 are Tau and Reset, which defines the epochs where the direction is
147 reset to the steepest descent.
148 
149 #### Broyden, Fletcher, Goldfarb, Shanno (BFGS) method:
150 
151  Implies the computation of a NxN matrix
152 computation, but seems more powerful at least for less than 300
153 weights. Parameters are Tau and Reset, which defines the epochs where
154 the direction is reset to the steepest descent.
155 
156 ### How to use it...
157 
158 TMLP is build from 3 classes: TNeuron, TSynapse and
159 TMultiLayerPerceptron. Only TMultiLayerPerceptron should be used
160 explicitly by the user.
161 
162 TMultiLayerPerceptron will take examples from a TTree
163 given in the constructor. The network is described by a simple
164 string: The input/output layers are defined by giving the expression for
165 each neuron, separated by comas. Hidden layers are just described
166 by the number of neurons. The layers are separated by colons.
167 In addition, input/output layer formulas can be preceded by '@' (e.g "@out")
168 if one wants to also normalize the data from the TTree.
169 Input and outputs are taken from the TTree given as second argument.
170 Expressions are evaluated as for TTree::Draw(), arrays are expended in
171 distinct neurons, one for each index.
172 This can only be done for fixed-size arrays.
173 If the formula ends with "!", softmax functions are used for the output layer.
174 One defines the training and test datasets by TEventLists.
175 
176 Example:
177 ~~~ {.cpp}
178 TMultiLayerPerceptron("x,y:10:5:f",inputTree);
179 ~~~
180 
181 Both the TTree and the TEventLists can be defined in
182 the constructor, or later with the suited setter method. The lists
183 used for training and test can be defined either explicitly, or via
184 a string containing the formula to be used to define them, exactly as
185 for a TCut.
186 
187 The learning method is defined using the TMultiLayerPerceptron::SetLearningMethod() .
188 Learning methods are :
189 
190  - TMultiLayerPerceptron::kStochastic,
191  - TMultiLayerPerceptron::kBatch,
192  - TMultiLayerPerceptron::kSteepestDescent,
193  - TMultiLayerPerceptron::kRibierePolak,
194  - TMultiLayerPerceptron::kFletcherReeves,
195  - TMultiLayerPerceptron::kBFGS
196 
197 A weight can be assigned to events, either in the constructor, either
198 with TMultiLayerPerceptron::SetEventWeight(). In addition, the TTree weight
199 is taken into account.
200 
201 Finally, one starts the training with
202 TMultiLayerPerceptron::Train(Int_t nepoch, Option_t* options). The
203 first argument is the number of epochs while option is a string that
204 can contain: "text" (simple text output) , "graph"
205 (evoluting graphical training curves), "update=X" (step for
206 the text/graph output update) or "+" (will skip the
207 randomisation and start from the previous values). All combinations
208 are available.
209 
210 Example:
211 ~~~ {.cpp}
212 net.Train(100,"text, graph, update=10");
213 ~~~
214 
215 When the neural net is trained, it can be used
216 directly ( TMultiLayerPerceptron::Evaluate() ) or exported to a
217 standalone C++ code ( TMultiLayerPerceptron::Export() ).
218 
219 Finally, note that even if this implementation is inspired from the mlpfit code,
220 the feature lists are not exactly matching:
221 
222  - mlpfit hybrid learning method is not implemented
223  - output neurons can be normalized, this is not the case for mlpfit
224  - the neural net is exported in C++, FORTRAN or PYTHON
225  - the drawResult() method allows a fast check of the learning procedure
226 
227 In addition, the paw version of mlpfit had additional limitations on the number of
228 neurons, hidden layers and inputs/outputs that does not apply to TMultiLayerPerceptron.
229 */
230 
231 
232 #include "TMultiLayerPerceptron.h"
233 #include "TSynapse.h"
234 #include "TNeuron.h"
235 #include "TClass.h"
236 #include "TTree.h"
237 #include "TEventList.h"
238 #include "TRandom3.h"
239 #include "TTimeStamp.h"
240 #include "TRegexp.h"
241 #include "TCanvas.h"
242 #include "TH2.h"
243 #include "TGraph.h"
244 #include "TLegend.h"
245 #include "TMultiGraph.h"
246 #include "TDirectory.h"
247 #include "TSystem.h"
248 #include "Riostream.h"
249 #include "TMath.h"
250 #include "TTreeFormula.h"
251 #include "TTreeFormulaManager.h"
252 #include "TMarker.h"
253 #include "TLine.h"
254 #include "TText.h"
255 #include "TObjString.h"
256 #include <stdlib.h>
257 
259 
260 ////////////////////////////////////////////////////////////////////////////////
261 /// Default constructor
262 
264 {
265  if(!TClass::GetClass("TTreePlayer")) gSystem->Load("libTreePlayer");
266  fNetwork.SetOwner(true);
267  fFirstLayer.SetOwner(false);
268  fLastLayer.SetOwner(false);
269  fSynapses.SetOwner(true);
270  fData = 0;
271  fCurrentTree = -1;
272  fCurrentTreeWeight = 1;
273  fStructure = "";
274  fWeight = "1";
275  fTraining = 0;
276  fTrainingOwner = false;
277  fTest = 0;
278  fTestOwner = false;
279  fEventWeight = 0;
280  fManager = 0;
282  fEta = .1;
283  fEtaDecay = 1;
284  fDelta = 0;
285  fEpsilon = 0;
286  fTau = 3;
287  fLastAlpha = 0;
288  fReset = 50;
291  fextF = "";
292  fextD = "";
293 }
294 
295 ////////////////////////////////////////////////////////////////////////////////
296 /// The network is described by a simple string:
297 /// The input/output layers are defined by giving
298 /// the branch names separated by comas.
299 /// Hidden layers are just described by the number of neurons.
300 /// The layers are separated by colons.
301 ///
302 /// Ex: "x,y:10:5:f"
303 ///
304 /// The output can be prepended by '@' if the variable has to be
305 /// normalized.
306 /// The output can be followed by '!' to use Softmax neurons for the
307 /// output layer only.
308 ///
309 /// Ex: "x,y:10:5:c1,c2,c3!"
310 ///
311 /// Input and outputs are taken from the TTree given as second argument.
312 /// training and test are the two TEventLists defining events
313 /// to be used during the neural net training.
314 /// Both the TTree and the TEventLists can be defined in the constructor,
315 /// or later with the suited setter method.
316 
318  TEventList * training,
319  TEventList * test,
321  const char* extF, const char* extD)
322 {
323  if(!TClass::GetClass("TTreePlayer")) gSystem->Load("libTreePlayer");
324  fNetwork.SetOwner(true);
325  fFirstLayer.SetOwner(false);
326  fLastLayer.SetOwner(false);
327  fSynapses.SetOwner(true);
328  fStructure = layout;
329  fData = data;
330  fCurrentTree = -1;
331  fCurrentTreeWeight = 1;
332  fTraining = training;
333  fTrainingOwner = false;
334  fTest = test;
335  fTestOwner = false;
336  fWeight = "1";
337  fType = type;
339  fextF = extF;
340  fextD = extD;
341  fEventWeight = 0;
342  fManager = 0;
343  if (data) {
344  BuildNetwork();
345  AttachData();
346  }
348  fEta = .1;
349  fEpsilon = 0;
350  fDelta = 0;
351  fEtaDecay = 1;
352  fTau = 3;
353  fLastAlpha = 0;
354  fReset = 50;
355 }
356 
357 ////////////////////////////////////////////////////////////////////////////////
358 /// The network is described by a simple string:
359 /// The input/output layers are defined by giving
360 /// the branch names separated by comas.
361 /// Hidden layers are just described by the number of neurons.
362 /// The layers are separated by colons.
363 ///
364 /// Ex: "x,y:10:5:f"
365 ///
366 /// The output can be prepended by '@' if the variable has to be
367 /// normalized.
368 /// The output can be followed by '!' to use Softmax neurons for the
369 /// output layer only.
370 ///
371 /// Ex: "x,y:10:5:c1,c2,c3!"
372 ///
373 /// Input and outputs are taken from the TTree given as second argument.
374 /// training and test are the two TEventLists defining events
375 /// to be used during the neural net training.
376 /// Both the TTree and the TEventLists can be defined in the constructor,
377 /// or later with the suited setter method.
378 
380  const char * weight, TTree * data,
381  TEventList * training,
382  TEventList * test,
384  const char* extF, const char* extD)
385 {
386  if(!TClass::GetClass("TTreePlayer")) gSystem->Load("libTreePlayer");
387  fNetwork.SetOwner(true);
388  fFirstLayer.SetOwner(false);
389  fLastLayer.SetOwner(false);
390  fSynapses.SetOwner(true);
391  fStructure = layout;
392  fData = data;
393  fCurrentTree = -1;
394  fCurrentTreeWeight = 1;
395  fTraining = training;
396  fTrainingOwner = false;
397  fTest = test;
398  fTestOwner = false;
399  fWeight = weight;
400  fType = type;
402  fextF = extF;
403  fextD = extD;
404  fEventWeight = 0;
405  fManager = 0;
406  if (data) {
407  BuildNetwork();
408  AttachData();
409  }
411  fEta = .1;
412  fEtaDecay = 1;
413  fDelta = 0;
414  fEpsilon = 0;
415  fTau = 3;
416  fLastAlpha = 0;
417  fReset = 50;
418 }
419 
420 ////////////////////////////////////////////////////////////////////////////////
421 /// The network is described by a simple string:
422 /// The input/output layers are defined by giving
423 /// the branch names separated by comas.
424 /// Hidden layers are just described by the number of neurons.
425 /// The layers are separated by colons.
426 ///
427 /// Ex: "x,y:10:5:f"
428 ///
429 /// The output can be prepended by '@' if the variable has to be
430 /// normalized.
431 /// The output can be followed by '!' to use Softmax neurons for the
432 /// output layer only.
433 ///
434 /// Ex: "x,y:10:5:c1,c2,c3!"
435 ///
436 /// Input and outputs are taken from the TTree given as second argument.
437 /// training and test are two cuts (see TTreeFormula) defining events
438 /// to be used during the neural net training and testing.
439 ///
440 /// Example: "Entry$%2", "(Entry$+1)%2".
441 ///
442 /// Both the TTree and the cut can be defined in the constructor,
443 /// or later with the suited setter method.
444 
446  const char * training,
447  const char * test,
449  const char* extF, const char* extD)
450 {
451  if(!TClass::GetClass("TTreePlayer")) gSystem->Load("libTreePlayer");
452  fNetwork.SetOwner(true);
453  fFirstLayer.SetOwner(false);
454  fLastLayer.SetOwner(false);
455  fSynapses.SetOwner(true);
456  fStructure = layout;
457  fData = data;
458  fCurrentTree = -1;
459  fCurrentTreeWeight = 1;
460  fTraining = new TEventList(Form("fTrainingList_%lu",(ULong_t)this));
461  fTrainingOwner = true;
462  fTest = new TEventList(Form("fTestList_%lu",(ULong_t)this));
463  fTestOwner = true;
464  fWeight = "1";
465  TString testcut = test;
466  if(testcut=="") testcut = Form("!(%s)",training);
467  fType = type;
469  fextF = extF;
470  fextD = extD;
471  fEventWeight = 0;
472  fManager = 0;
473  if (data) {
474  BuildNetwork();
475  data->Draw(Form(">>fTrainingList_%lu",(ULong_t)this),training,"goff");
476  data->Draw(Form(">>fTestList_%lu",(ULong_t)this),(const char *)testcut,"goff");
477  AttachData();
478  }
479  else {
480  Warning("TMultiLayerPerceptron::TMultiLayerPerceptron","Data not set. Cannot define datasets");
481  }
483  fEta = .1;
484  fEtaDecay = 1;
485  fDelta = 0;
486  fEpsilon = 0;
487  fTau = 3;
488  fLastAlpha = 0;
489  fReset = 50;
490 }
491 
492 ////////////////////////////////////////////////////////////////////////////////
493 /// The network is described by a simple string:
494 /// The input/output layers are defined by giving
495 /// the branch names separated by comas.
496 /// Hidden layers are just described by the number of neurons.
497 /// The layers are separated by colons.
498 ///
499 /// Ex: "x,y:10:5:f"
500 ///
501 /// The output can be prepended by '@' if the variable has to be
502 /// normalized.
503 /// The output can be followed by '!' to use Softmax neurons for the
504 /// output layer only.
505 ///
506 /// Ex: "x,y:10:5:c1,c2,c3!"
507 ///
508 /// Input and outputs are taken from the TTree given as second argument.
509 /// training and test are two cuts (see TTreeFormula) defining events
510 /// to be used during the neural net training and testing.
511 ///
512 /// Example: "Entry$%2", "(Entry$+1)%2".
513 ///
514 /// Both the TTree and the cut can be defined in the constructor,
515 /// or later with the suited setter method.
516 
518  const char * weight, TTree * data,
519  const char * training,
520  const char * test,
522  const char* extF, const char* extD)
523 {
524  if(!TClass::GetClass("TTreePlayer")) gSystem->Load("libTreePlayer");
525  fNetwork.SetOwner(true);
526  fFirstLayer.SetOwner(false);
527  fLastLayer.SetOwner(false);
528  fSynapses.SetOwner(true);
529  fStructure = layout;
530  fData = data;
531  fCurrentTree = -1;
532  fCurrentTreeWeight = 1;
533  fTraining = new TEventList(Form("fTrainingList_%lu",(ULong_t)this));
534  fTrainingOwner = true;
535  fTest = new TEventList(Form("fTestList_%lu",(ULong_t)this));
536  fTestOwner = true;
537  fWeight = weight;
538  TString testcut = test;
539  if(testcut=="") testcut = Form("!(%s)",training);
540  fType = type;
542  fextF = extF;
543  fextD = extD;
544  fEventWeight = 0;
545  fManager = 0;
546  if (data) {
547  BuildNetwork();
548  data->Draw(Form(">>fTrainingList_%lu",(ULong_t)this),training,"goff");
549  data->Draw(Form(">>fTestList_%lu",(ULong_t)this),(const char *)testcut,"goff");
550  AttachData();
551  }
552  else {
553  Warning("TMultiLayerPerceptron::TMultiLayerPerceptron","Data not set. Cannot define datasets");
554  }
556  fEta = .1;
557  fEtaDecay = 1;
558  fDelta = 0;
559  fEpsilon = 0;
560  fTau = 3;
561  fLastAlpha = 0;
562  fReset = 50;
563 }
564 
565 ////////////////////////////////////////////////////////////////////////////////
566 /// Destructor
567 
569 {
570  if(fTraining && fTrainingOwner) delete fTraining;
571  if(fTest && fTestOwner) delete fTest;
572 }
573 
574 ////////////////////////////////////////////////////////////////////////////////
575 /// Set the data source
576 
578 {
579  if (fData) {
580  std::cerr << "Error: data already defined." << std::endl;
581  return;
582  }
583  fData = data;
584  if (data) {
585  BuildNetwork();
586  AttachData();
587  }
588 }
589 
590 ////////////////////////////////////////////////////////////////////////////////
591 /// Set the event weight
592 
593 void TMultiLayerPerceptron::SetEventWeight(const char * branch)
594 {
595  fWeight=branch;
596  if (fData) {
597  if (fEventWeight) {
599  delete fEventWeight;
600  }
601  fManager->Add((fEventWeight = new TTreeFormula("NNweight",fWeight.Data(),fData)));
602  }
603 }
604 
605 ////////////////////////////////////////////////////////////////////////////////
606 /// Sets the Training dataset.
607 /// Those events will be used for the minimization
608 
610 {
611  if(fTraining && fTrainingOwner) delete fTraining;
612  fTraining = train;
613  fTrainingOwner = false;
614 }
615 
616 ////////////////////////////////////////////////////////////////////////////////
617 /// Sets the Test dataset.
618 /// Those events will not be used for the minimization but for control
619 
621 {
622  if(fTest && fTestOwner) delete fTest;
623  fTest = test;
624  fTestOwner = false;
625 }
626 
627 ////////////////////////////////////////////////////////////////////////////////
628 /// Sets the Training dataset.
629 /// Those events will be used for the minimization.
630 /// Note that the tree must be already defined.
631 
633 {
634  if(fTraining && fTrainingOwner) delete fTraining;
635  fTraining = new TEventList(Form("fTrainingList_%lu",(ULong_t)this));
636  fTrainingOwner = true;
637  if (fData) {
638  fData->Draw(Form(">>fTrainingList_%lu",(ULong_t)this),train,"goff");
639  }
640  else {
641  Warning("TMultiLayerPerceptron::TMultiLayerPerceptron","Data not set. Cannot define datasets");
642  }
643 }
644 
645 ////////////////////////////////////////////////////////////////////////////////
646 /// Sets the Test dataset.
647 /// Those events will not be used for the minimization but for control.
648 /// Note that the tree must be already defined.
649 
651 {
652  if(fTest && fTestOwner) {delete fTest; fTest=0;}
653  if(fTest) if(strncmp(fTest->GetName(),Form("fTestList_%lu",(ULong_t)this),10)) delete fTest;
654  fTest = new TEventList(Form("fTestList_%lu",(ULong_t)this));
655  fTestOwner = true;
656  if (fData) {
657  fData->Draw(Form(">>fTestList_%lu",(ULong_t)this),test,"goff");
658  }
659  else {
660  Warning("TMultiLayerPerceptron::TMultiLayerPerceptron","Data not set. Cannot define datasets");
661  }
662 }
663 
664 ////////////////////////////////////////////////////////////////////////////////
665 /// Sets the learning method.
666 /// Available methods are: kStochastic, kBatch,
667 /// kSteepestDescent, kRibierePolak, kFletcherReeves and kBFGS.
668 /// (look at the constructor for the complete description
669 /// of learning methods and parameters)
670 
672 {
673  fLearningMethod = method;
674 }
675 
676 ////////////////////////////////////////////////////////////////////////////////
677 /// Sets Eta - used in stochastic minimisation
678 /// (look at the constructor for the complete description
679 /// of learning methods and parameters)
680 
682 {
683  fEta = eta;
684 }
685 
686 ////////////////////////////////////////////////////////////////////////////////
687 /// Sets Epsilon - used in stochastic minimisation
688 /// (look at the constructor for the complete description
689 /// of learning methods and parameters)
690 
692 {
693  fEpsilon = eps;
694 }
695 
696 ////////////////////////////////////////////////////////////////////////////////
697 /// Sets Delta - used in stochastic minimisation
698 /// (look at the constructor for the complete description
699 /// of learning methods and parameters)
700 
702 {
703  fDelta = delta;
704 }
705 
706 ////////////////////////////////////////////////////////////////////////////////
707 /// Sets EtaDecay - Eta *= EtaDecay at each epoch
708 /// (look at the constructor for the complete description
709 /// of learning methods and parameters)
710 
712 {
713  fEtaDecay = ed;
714 }
715 
716 ////////////////////////////////////////////////////////////////////////////////
717 /// Sets Tau - used in line search
718 /// (look at the constructor for the complete description
719 /// of learning methods and parameters)
720 
722 {
723  fTau = tau;
724 }
725 
726 ////////////////////////////////////////////////////////////////////////////////
727 /// Sets number of epochs between two resets of the
728 /// search direction to the steepest descent.
729 /// (look at the constructor for the complete description
730 /// of learning methods and parameters)
731 
733 {
734  fReset = reset;
735 }
736 
737 ////////////////////////////////////////////////////////////////////////////////
738 /// Load an entry into the network
739 
741 {
742  if (!fData) return;
743  fData->GetEntry(entry);
744  if (fData->GetTreeNumber() != fCurrentTree) {
745  ((TMultiLayerPerceptron*)this)->fCurrentTree = fData->GetTreeNumber();
746  fManager->Notify();
747  ((TMultiLayerPerceptron*)this)->fCurrentTreeWeight = fData->GetWeight();
748  }
750  for (Int_t i=0;i<nentries;i++) {
751  TNeuron *neuron = (TNeuron *)fNetwork.UncheckedAt(i);
752  neuron->SetNewEvent();
753  }
754 }
755 
756 ////////////////////////////////////////////////////////////////////////////////
757 /// Train the network.
758 /// nEpoch is the number of iterations.
759 /// option can contain:
760 /// - "text" (simple text output)
761 /// - "graph" (evoluting graphical training curves)
762 /// - "update=X" (step for the text/graph output update)
763 /// - "+" will skip the randomisation and start from the previous values.
764 /// - "current" (draw in the current canvas)
765 /// - "minErrorTrain" (stop when NN error on the training sample gets below minE
766 /// - "minErrorTest" (stop when NN error on the test sample gets below minE
767 /// All combinations are available.
768 
770 {
771  Int_t i;
772  TString opt = option;
773  opt.ToLower();
774  // Decode options and prepare training.
775  Int_t verbosity = 0;
776  Bool_t newCanvas = true;
777  Bool_t minE_Train = false;
778  Bool_t minE_Test = false;
779  if (opt.Contains("text"))
780  verbosity += 1;
781  if (opt.Contains("graph"))
782  verbosity += 2;
783  Int_t displayStepping = 1;
784  if (opt.Contains("update=")) {
785  TRegexp reg("update=[0-9]*");
786  TString out = opt(reg);
787  displayStepping = atoi(out.Data() + 7);
788  }
789  if (opt.Contains("current"))
790  newCanvas = false;
791  if (opt.Contains("minerrortrain"))
792  minE_Train = true;
793  if (opt.Contains("minerrortest"))
794  minE_Test = true;
795  TVirtualPad *canvas = 0;
796  TMultiGraph *residual_plot = 0;
797  TGraph *train_residual_plot = 0;
798  TGraph *test_residual_plot = 0;
799  if ((!fData) || (!fTraining) || (!fTest)) {
800  Error("Train","Training/Test samples still not defined. Cannot train the neural network");
801  return;
802  }
803  Info("Train","Using %d train and %d test entries.",
804  fTraining->GetN(), fTest->GetN());
805  // Text and Graph outputs
806  if (verbosity % 2)
807  std::cout << "Training the Neural Network" << std::endl;
808  if (verbosity / 2) {
809  residual_plot = new TMultiGraph;
810  if(newCanvas)
811  canvas = new TCanvas("NNtraining", "Neural Net training");
812  else {
813  canvas = gPad;
814  if(!canvas) canvas = new TCanvas("NNtraining", "Neural Net training");
815  }
816  train_residual_plot = new TGraph(nEpoch);
817  test_residual_plot = new TGraph(nEpoch);
818  canvas->SetLeftMargin(0.14);
819  train_residual_plot->SetLineColor(4);
820  test_residual_plot->SetLineColor(2);
821  residual_plot->Add(train_residual_plot);
822  residual_plot->Add(test_residual_plot);
823  residual_plot->Draw("LA");
824  if (residual_plot->GetXaxis()) residual_plot->GetXaxis()->SetTitle("Epoch");
825  if (residual_plot->GetYaxis()) residual_plot->GetYaxis()->SetTitle("Error");
826  }
827  // If the option "+" is not set, one has to randomize the weights first
828  if (!opt.Contains("+"))
829  Randomize();
830  // Initialisation
831  fLastAlpha = 0;
833  Double_t *buffer = new Double_t[els];
834  Double_t *dir = new Double_t[els];
835  for (i = 0; i < els; i++)
836  buffer[i] = 0;
837  Int_t matrix_size = fLearningMethod==TMultiLayerPerceptron::kBFGS ? els : 1;
838  TMatrixD bfgsh(matrix_size, matrix_size);
839  TMatrixD gamma(matrix_size, 1);
840  TMatrixD delta(matrix_size, 1);
841  // Epoch loop. Here is the training itself.
842  Double_t training_E = 1e10;
843  Double_t test_E = 1e10;
844  for (Int_t iepoch = 0; (iepoch < nEpoch) && (!minE_Train || training_E>minE) && (!minE_Test || test_E>minE) ; iepoch++) {
845  switch (fLearningMethod) {
847  {
848  MLP_Stochastic(buffer);
849  break;
850  }
852  {
853  ComputeDEDw();
854  MLP_Batch(buffer);
855  break;
856  }
858  {
859  ComputeDEDw();
860  SteepestDir(dir);
861  if (LineSearch(dir, buffer))
862  MLP_Batch(buffer);
863  break;
864  }
866  {
867  ComputeDEDw();
868  if (!(iepoch % fReset)) {
869  SteepestDir(dir);
870  } else {
871  Double_t norm = 0;
872  Double_t onorm = 0;
873  for (i = 0; i < els; i++)
874  onorm += dir[i] * dir[i];
875  Double_t prod = 0;
876  Int_t idx = 0;
877  TNeuron *neuron = 0;
878  TSynapse *synapse = 0;
880  for (i=0;i<nentries;i++) {
881  neuron = (TNeuron *) fNetwork.UncheckedAt(i);
882  prod -= dir[idx++] * neuron->GetDEDw();
883  norm += neuron->GetDEDw() * neuron->GetDEDw();
884  }
886  for (i=0;i<nentries;i++) {
887  synapse = (TSynapse *) fSynapses.UncheckedAt(i);
888  prod -= dir[idx++] * synapse->GetDEDw();
889  norm += synapse->GetDEDw() * synapse->GetDEDw();
890  }
891  ConjugateGradientsDir(dir, (norm - prod) / onorm);
892  }
893  if (LineSearch(dir, buffer))
894  MLP_Batch(buffer);
895  break;
896  }
898  {
899  ComputeDEDw();
900  if (!(iepoch % fReset)) {
901  SteepestDir(dir);
902  } else {
903  Double_t norm = 0;
904  Double_t onorm = 0;
905  for (i = 0; i < els; i++)
906  onorm += dir[i] * dir[i];
907  TNeuron *neuron = 0;
908  TSynapse *synapse = 0;
910  for (i=0;i<nentries;i++) {
911  neuron = (TNeuron *) fNetwork.UncheckedAt(i);
912  norm += neuron->GetDEDw() * neuron->GetDEDw();
913  }
915  for (i=0;i<nentries;i++) {
916  synapse = (TSynapse *) fSynapses.UncheckedAt(i);
917  norm += synapse->GetDEDw() * synapse->GetDEDw();
918  }
919  ConjugateGradientsDir(dir, norm / onorm);
920  }
921  if (LineSearch(dir, buffer))
922  MLP_Batch(buffer);
923  break;
924  }
926  {
927  SetGammaDelta(gamma, delta, buffer);
928  if (!(iepoch % fReset)) {
929  SteepestDir(dir);
930  bfgsh.UnitMatrix();
931  } else {
932  if (GetBFGSH(bfgsh, gamma, delta)) {
933  SteepestDir(dir);
934  bfgsh.UnitMatrix();
935  } else {
936  BFGSDir(bfgsh, dir);
937  }
938  }
939  if (DerivDir(dir) > 0) {
940  SteepestDir(dir);
941  bfgsh.UnitMatrix();
942  }
943  if (LineSearch(dir, buffer)) {
944  bfgsh.UnitMatrix();
945  SteepestDir(dir);
946  if (LineSearch(dir, buffer)) {
947  Error("TMultiLayerPerceptron::Train()","Line search fail");
948  iepoch = nEpoch;
949  }
950  }
951  break;
952  }
953  }
954  // Security: would the learning lead to non real numbers,
955  // the learning should stop now.
957  Error("TMultiLayerPerceptron::Train()","Stop.");
958  iepoch = nEpoch;
959  }
960  // Process other ROOT events. Time penalty is less than
961  // 1/1000 sec/evt on a mobile AMD Athlon(tm) XP 1500+
965  // Intermediate graph and text output
966  if ((verbosity % 2) && ((!(iepoch % displayStepping)) || (iepoch == nEpoch - 1))) {
967  std::cout << "Epoch: " << iepoch
968  << " learn=" << training_E
969  << " test=" << test_E
970  << std::endl;
971  }
972  if (verbosity / 2) {
973  train_residual_plot->SetPoint(iepoch, iepoch,training_E);
974  test_residual_plot->SetPoint(iepoch, iepoch,test_E);
975  if (!iepoch) {
976  Double_t trp = train_residual_plot->GetY()[iepoch];
977  Double_t tep = test_residual_plot->GetY()[iepoch];
978  for (i = 1; i < nEpoch; i++) {
979  train_residual_plot->SetPoint(i, i, trp);
980  test_residual_plot->SetPoint(i, i, tep);
981  }
982  }
983  if ((!(iepoch % displayStepping)) || (iepoch == nEpoch - 1)) {
984  if (residual_plot->GetYaxis()) {
985  residual_plot->GetYaxis()->UnZoom();
986  residual_plot->GetYaxis()->SetTitleOffset(1.4);
987  residual_plot->GetYaxis()->SetDecimals();
988  }
989  canvas->Modified();
990  canvas->Update();
991  }
992  }
993  }
994  // Cleaning
995  delete [] buffer;
996  delete [] dir;
997  // Final Text and Graph outputs
998  if (verbosity % 2)
999  std::cout << "Training done." << std::endl;
1000  if (verbosity / 2) {
1001  TLegend *legend = new TLegend(.75, .80, .95, .95);
1002  legend->AddEntry(residual_plot->GetListOfGraphs()->At(0),
1003  "Training sample", "L");
1004  legend->AddEntry(residual_plot->GetListOfGraphs()->At(1),
1005  "Test sample", "L");
1006  legend->Draw();
1007  }
1008 }
1009 
1010 ////////////////////////////////////////////////////////////////////////////////
1011 /// Computes the output for a given event.
1012 /// Look at the output neuron designed by index.
1013 
1015 {
1016  GetEntry(event);
1017  TNeuron *out = (TNeuron *) (fLastLayer.At(index));
1018  if (out)
1019  return out->GetValue();
1020  else
1021  return 0;
1022 }
1023 
1024 ////////////////////////////////////////////////////////////////////////////////
1025 /// Error on the output for a given event
1026 
1028 {
1029  GetEntry(event);
1030  Double_t error = 0;
1031  // look at 1st output neuron to determine type and error function
1032  Int_t nEntries = fLastLayer.GetEntriesFast();
1033  if (nEntries == 0) return 0.0;
1034  switch (fOutType) {
1035  case (TNeuron::kSigmoid):
1036  error = GetCrossEntropyBinary();
1037  break;
1038  case (TNeuron::kSoftmax):
1039  error = GetCrossEntropy();
1040  break;
1041  case (TNeuron::kLinear):
1042  error = GetSumSquareError();
1043  break;
1044  default:
1045  // default to sum-of-squares error
1046  error = GetSumSquareError();
1047  }
1048  error *= fEventWeight->EvalInstance();
1049  error *= fCurrentTreeWeight;
1050  return error;
1051 }
1052 
1053 ////////////////////////////////////////////////////////////////////////////////
1054 /// Error on the whole dataset
1055 
1057 {
1058  TEventList *list =
1060  Double_t error = 0;
1061  Int_t i;
1062  if (list) {
1063  Int_t nEvents = list->GetN();
1064  for (i = 0; i < nEvents; i++) {
1065  error += GetError(list->GetEntry(i));
1066  }
1067  } else if (fData) {
1068  Int_t nEvents = (Int_t) fData->GetEntries();
1069  for (i = 0; i < nEvents; i++) {
1070  error += GetError(i);
1071  }
1072  }
1073  return error;
1074 }
1075 
1076 ////////////////////////////////////////////////////////////////////////////////
1077 /// Error on the output for a given event
1078 
1080 {
1081  Double_t error = 0;
1082  for (Int_t i = 0; i < fLastLayer.GetEntriesFast(); i++) {
1083  TNeuron *neuron = (TNeuron *) fLastLayer[i];
1084  error += neuron->GetError() * neuron->GetError();
1085  }
1086  return (error / 2.);
1087 }
1088 
1089 ////////////////////////////////////////////////////////////////////////////////
1090 /// Cross entropy error for sigmoid output neurons, for a given event
1091 
1093 {
1094  Double_t error = 0;
1095  for (Int_t i = 0; i < fLastLayer.GetEntriesFast(); i++) {
1096  TNeuron *neuron = (TNeuron *) fLastLayer[i];
1097  Double_t output = neuron->GetValue(); // sigmoid output and target
1098  Double_t target = neuron->GetTarget(); // values lie in [0,1]
1099  if (target < DBL_EPSILON) {
1100  if (output == 1.0)
1101  error = DBL_MAX;
1102  else
1103  error -= TMath::Log(1 - output);
1104  } else
1105  if ((1 - target) < DBL_EPSILON) {
1106  if (output == 0.0)
1107  error = DBL_MAX;
1108  else
1109  error -= TMath::Log(output);
1110  } else {
1111  if (output == 0.0 || output == 1.0)
1112  error = DBL_MAX;
1113  else
1114  error -= target * TMath::Log(output / target) + (1-target) * TMath::Log((1 - output)/(1 - target));
1115  }
1116  }
1117  return error;
1118 }
1119 
1120 ////////////////////////////////////////////////////////////////////////////////
1121 /// Cross entropy error for a softmax output neuron, for a given event
1122 
1124 {
1125  Double_t error = 0;
1126  for (Int_t i = 0; i < fLastLayer.GetEntriesFast(); i++) {
1127  TNeuron *neuron = (TNeuron *) fLastLayer[i];
1128  Double_t output = neuron->GetValue(); // softmax output and target
1129  Double_t target = neuron->GetTarget(); // values lie in [0,1]
1130  if (target > DBL_EPSILON) { // (target == 0) => dE = 0
1131  if (output == 0.0)
1132  error = DBL_MAX;
1133  else
1134  error -= target * TMath::Log(output / target);
1135  }
1136  }
1137  return error;
1138 }
1139 
1140 ////////////////////////////////////////////////////////////////////////////////
1141 /// Compute the DEDw = sum on all training events of dedw for each weight
1142 /// normalized by the number of events.
1143 
1145 {
1146  Int_t i,j;
1148  TSynapse *synapse;
1149  for (i=0;i<nentries;i++) {
1150  synapse = (TSynapse *) fSynapses.UncheckedAt(i);
1151  synapse->SetDEDw(0.);
1152  }
1153  TNeuron *neuron;
1155  for (i=0;i<nentries;i++) {
1156  neuron = (TNeuron *) fNetwork.UncheckedAt(i);
1157  neuron->SetDEDw(0.);
1158  }
1159  Double_t eventWeight = 1.;
1160  if (fTraining) {
1161  Int_t nEvents = fTraining->GetN();
1162  for (i = 0; i < nEvents; i++) {
1164  eventWeight = fEventWeight->EvalInstance();
1165  eventWeight *= fCurrentTreeWeight;
1167  for (j=0;j<nentries;j++) {
1168  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
1169  synapse->SetDEDw(synapse->GetDEDw() + (synapse->GetDeDw()*eventWeight));
1170  }
1172  for (j=0;j<nentries;j++) {
1173  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
1174  neuron->SetDEDw(neuron->GetDEDw() + (neuron->GetDeDw()*eventWeight));
1175  }
1176  }
1178  for (j=0;j<nentries;j++) {
1179  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
1180  synapse->SetDEDw(synapse->GetDEDw() / (Double_t) nEvents);
1181  }
1183  for (j=0;j<nentries;j++) {
1184  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
1185  neuron->SetDEDw(neuron->GetDEDw() / (Double_t) nEvents);
1186  }
1187  } else if (fData) {
1188  Int_t nEvents = (Int_t) fData->GetEntries();
1189  for (i = 0; i < nEvents; i++) {
1190  GetEntry(i);
1191  eventWeight = fEventWeight->EvalInstance();
1192  eventWeight *= fCurrentTreeWeight;
1194  for (j=0;j<nentries;j++) {
1195  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
1196  synapse->SetDEDw(synapse->GetDEDw() + (synapse->GetDeDw()*eventWeight));
1197  }
1199  for (j=0;j<nentries;j++) {
1200  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
1201  neuron->SetDEDw(neuron->GetDEDw() + (neuron->GetDeDw()*eventWeight));
1202  }
1203  }
1205  for (j=0;j<nentries;j++) {
1206  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
1207  synapse->SetDEDw(synapse->GetDEDw() / (Double_t) nEvents);
1208  }
1210  for (j=0;j<nentries;j++) {
1211  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
1212  neuron->SetDEDw(neuron->GetDEDw() / (Double_t) nEvents);
1213  }
1214  }
1215 }
1216 
1217 ////////////////////////////////////////////////////////////////////////////////
1218 /// Randomize the weights
1219 
1221 {
1223  Int_t j;
1224  TSynapse *synapse;
1225  TNeuron *neuron;
1226  TTimeStamp ts;
1227  TRandom3 gen(ts.GetSec());
1228  for (j=0;j<nentries;j++) {
1229  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
1230  synapse->SetWeight(gen.Rndm() - 0.5);
1231  }
1233  for (j=0;j<nentries;j++) {
1234  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
1235  neuron->SetWeight(gen.Rndm() - 0.5);
1236  }
1237 }
1238 
1239 ////////////////////////////////////////////////////////////////////////////////
1240 /// Connects the TTree to Neurons in input and output
1241 /// layers. The formulas associated to each neuron are created
1242 /// and reported to the network formula manager.
1243 /// By default, the branch is not normalised since this would degrade
1244 /// performance for classification jobs.
1245 /// Normalisation can be requested by putting '@' in front of the formula.
1246 
1248 {
1249  Int_t j = 0;
1250  TNeuron *neuron = 0;
1251  Bool_t normalize = false;
1253 
1254  // Set the size of the internal array of parameters of the formula
1255  Int_t maxop, maxpar, maxconst;
1256  ROOT::v5::TFormula::GetMaxima(maxop, maxpar, maxconst);
1257  ROOT::v5::TFormula::SetMaxima(10, 10, 10);
1258 
1259  //first layer
1260  const TString input = TString(fStructure(0, fStructure.First(':')));
1261  const TObjArray *inpL = input.Tokenize(", ");
1263  // make sure nentries == entries in inpL
1264  R__ASSERT(nentries == inpL->GetLast()+1);
1265  for (j=0;j<nentries;j++) {
1266  normalize = false;
1267  const TString brName = ((TObjString *)inpL->At(j))->GetString();
1268  neuron = (TNeuron *) fFirstLayer.UncheckedAt(j);
1269  if (brName[0]=='@')
1270  normalize = true;
1271  fManager->Add(neuron->UseBranch(fData,brName.Data() + (normalize?1:0)));
1272  if(!normalize) neuron->SetNormalisation(0., 1.);
1273  }
1274  delete inpL;
1275 
1276  // last layer
1278  fStructure(fStructure.Last(':') + 1,
1279  fStructure.Length() - fStructure.Last(':')));
1280  const TObjArray *outL = output.Tokenize(", ");
1282  // make sure nentries == entries in outL
1283  R__ASSERT(nentries == outL->GetLast()+1);
1284  for (j=0;j<nentries;j++) {
1285  normalize = false;
1286  const TString brName = ((TObjString *)outL->At(j))->GetString();
1287  neuron = (TNeuron *) fLastLayer.UncheckedAt(j);
1288  if (brName[0]=='@')
1289  normalize = true;
1290  fManager->Add(neuron->UseBranch(fData,brName.Data() + (normalize?1:0)));
1291  if(!normalize) neuron->SetNormalisation(0., 1.);
1292  }
1293  delete outL;
1294 
1295  fManager->Add((fEventWeight = new TTreeFormula("NNweight",fWeight.Data(),fData)));
1296  //fManager->Sync();
1297 
1298  // Set the old values
1299  ROOT::v5::TFormula::SetMaxima(maxop, maxpar, maxconst);
1300 }
1301 
1302 ////////////////////////////////////////////////////////////////////////////////
1303 /// Expand the structure of the first layer
1304 
1306 {
1307  TString input = TString(fStructure(0, fStructure.First(':')));
1308  const TObjArray *inpL = input.Tokenize(", ");
1309  Int_t nneurons = inpL->GetLast()+1;
1310 
1311  TString hiddenAndOutput = TString(
1312  fStructure(fStructure.First(':') + 1,
1313  fStructure.Length() - fStructure.First(':')));
1314  TString newInput;
1315  Int_t i = 0;
1316  // loop on input neurons
1317  for (i = 0; i<nneurons; i++) {
1318  const TString name = ((TObjString *)inpL->At(i))->GetString();
1319  TTreeFormula f("sizeTestFormula",name,fData);
1320  // Variable size arrays are unreliable
1321  if(f.GetMultiplicity()==1 && f.GetNdata()>1) {
1322  Warning("TMultiLayerPerceptron::ExpandStructure()","Variable size arrays cannot be used to build implicitly an input layer. The index 0 will be assumed.");
1323  }
1324  // Check if we are coping with an array... then expand
1325  // The array operator used is {}. It is detected in TNeuron, and
1326  // passed directly as instance index of the TTreeFormula,
1327  // so that complex compounds made of arrays can be used without
1328  // parsing the details.
1329  else if(f.GetNdata()>1) {
1330  for(Int_t j=0; j<f.GetNdata(); j++) {
1331  if(i||j) newInput += ",";
1332  newInput += name;
1333  newInput += "{";
1334  newInput += j;
1335  newInput += "}";
1336  }
1337  continue;
1338  }
1339  if(i) newInput += ",";
1340  newInput += name;
1341  }
1342  delete inpL;
1343 
1344  // Save the result
1345  fStructure = newInput + ":" + hiddenAndOutput;
1346 }
1347 
1348 ////////////////////////////////////////////////////////////////////////////////
1349 /// Instantiates the network from the description
1350 
1352 {
1353  ExpandStructure();
1354  TString input = TString(fStructure(0, fStructure.First(':')));
1355  TString hidden = TString(
1356  fStructure(fStructure.First(':') + 1,
1357  fStructure.Last(':') - fStructure.First(':') - 1));
1359  fStructure(fStructure.Last(':') + 1,
1360  fStructure.Length() - fStructure.Last(':')));
1361  Int_t bll = atoi(TString(
1362  hidden(hidden.Last(':') + 1,
1363  hidden.Length() - (hidden.Last(':') + 1))).Data());
1364  if (input.Length() == 0) {
1365  Error("BuildNetwork()","malformed structure. No input layer.");
1366  return;
1367  }
1368  if (output.Length() == 0) {
1369  Error("BuildNetwork()","malformed structure. No output layer.");
1370  return;
1371  }
1372  BuildFirstLayer(input);
1373  BuildHiddenLayers(hidden);
1374  BuildLastLayer(output, bll);
1375 }
1376 
1377 ////////////////////////////////////////////////////////////////////////////////
1378 /// Instantiates the neurons in input
1379 /// Inputs are normalised and the type is set to kOff
1380 /// (simple forward of the formula value)
1381 
1383 {
1384  const TObjArray *inpL = input.Tokenize(", ");
1385  const Int_t nneurons =inpL->GetLast()+1;
1386  TNeuron *neuron = 0;
1387  Int_t i = 0;
1388  for (i = 0; i<nneurons; i++) {
1389  const TString name = ((TObjString *)inpL->At(i))->GetString();
1390  neuron = new TNeuron(TNeuron::kOff, name);
1391  fFirstLayer.AddLast(neuron);
1392  fNetwork.AddLast(neuron);
1393  }
1394  delete inpL;
1395 }
1396 
1397 ////////////////////////////////////////////////////////////////////////////////
1398 /// Builds hidden layers.
1399 
1401 {
1402  Int_t beg = 0;
1403  Int_t end = hidden.Index(":", beg + 1);
1404  Int_t prevStart = 0;
1405  Int_t prevStop = fNetwork.GetEntriesFast();
1406  Int_t layer = 1;
1407  while (end != -1) {
1408  BuildOneHiddenLayer(hidden(beg, end - beg), layer, prevStart, prevStop, false);
1409  beg = end + 1;
1410  end = hidden.Index(":", beg + 1);
1411  }
1412 
1413  BuildOneHiddenLayer(hidden(beg, hidden.Length() - beg), layer, prevStart, prevStop, true);
1414 }
1415 
1416 ////////////////////////////////////////////////////////////////////////////////
1417 /// Builds a hidden layer, updates the number of layers.
1418 
1420  Int_t& prevStart, Int_t& prevStop,
1421  Bool_t lastLayer)
1422 {
1423  TNeuron *neuron = 0;
1424  TSynapse *synapse = 0;
1425  TString name;
1426  if (!sNumNodes.IsAlnum() || sNumNodes.IsAlpha()) {
1427  Error("BuildOneHiddenLayer",
1428  "The specification '%s' for hidden layer %d must contain only numbers!",
1429  sNumNodes.Data(), layer - 1);
1430  } else {
1431  Int_t num = atoi(sNumNodes.Data());
1432  for (Int_t i = 0; i < num; i++) {
1433  name.Form("HiddenL%d:N%d",layer,i);
1434  neuron = new TNeuron(fType, name, "", (const char*)fextF, (const char*)fextD);
1435  fNetwork.AddLast(neuron);
1436  for (Int_t j = prevStart; j < prevStop; j++) {
1437  synapse = new TSynapse((TNeuron *) fNetwork[j], neuron);
1438  fSynapses.AddLast(synapse);
1439  }
1440  }
1441 
1442  if (!lastLayer) {
1443  // tell each neuron which ones are in its own layer (for Softmax)
1444  Int_t nEntries = fNetwork.GetEntriesFast();
1445  for (Int_t i = prevStop; i < nEntries; i++) {
1446  neuron = (TNeuron *) fNetwork[i];
1447  for (Int_t j = prevStop; j < nEntries; j++)
1448  neuron->AddInLayer((TNeuron *) fNetwork[j]);
1449  }
1450  }
1451 
1452  prevStart = prevStop;
1453  prevStop = fNetwork.GetEntriesFast();
1454  layer++;
1455  }
1456 }
1457 
1458 ////////////////////////////////////////////////////////////////////////////////
1459 /// Builds the output layer
1460 /// Neurons are linear combinations of input, by default.
1461 /// If the structure ends with "!", neurons are set up for classification,
1462 /// ie. with a sigmoid (1 neuron) or softmax (more neurons) activation function.
1463 
1465 {
1466  Int_t nneurons = output.CountChar(',')+1;
1467  if (fStructure.EndsWith("!")) {
1468  fStructure = TString(fStructure(0, fStructure.Length() - 1)); // remove "!"
1469  if (nneurons == 1)
1471  else
1473  }
1474  Int_t prevStop = fNetwork.GetEntriesFast();
1475  Int_t prevStart = prevStop - prev;
1476  Ssiz_t pos = 0;
1477  TNeuron *neuron;
1478  TSynapse *synapse;
1479  TString name;
1480  Int_t i,j;
1481  for (i = 0; i<nneurons; i++) {
1482  Ssiz_t nextpos=output.Index(",",pos);
1483  if (nextpos!=kNPOS)
1484  name=output(pos,nextpos-pos);
1485  else name=output(pos,output.Length());
1486  pos+=nextpos+1;
1487  neuron = new TNeuron(fOutType, name);
1488  for (j = prevStart; j < prevStop; j++) {
1489  synapse = new TSynapse((TNeuron *) fNetwork[j], neuron);
1490  fSynapses.AddLast(synapse);
1491  }
1492  fLastLayer.AddLast(neuron);
1493  fNetwork.AddLast(neuron);
1494  }
1495  // tell each neuron which ones are in its own layer (for Softmax)
1496  Int_t nEntries = fNetwork.GetEntriesFast();
1497  for (i = prevStop; i < nEntries; i++) {
1498  neuron = (TNeuron *) fNetwork[i];
1499  for (j = prevStop; j < nEntries; j++)
1500  neuron->AddInLayer((TNeuron *) fNetwork[j]);
1501  }
1502 
1503 }
1504 
1505 ////////////////////////////////////////////////////////////////////////////////
1506 /// Draws the neural net output
1507 /// It produces an histogram with the output for the two datasets.
1508 /// Index is the number of the desired output neuron.
1509 /// "option" can contain:
1510 /// - test or train to select a dataset
1511 /// - comp to produce a X-Y comparison plot
1512 /// - nocanv to not create a new TCanvas for the plot
1513 
1515 {
1516  TString opt = option;
1517  opt.ToLower();
1518  TNeuron *out = (TNeuron *) (fLastLayer.At(index));
1519  if (!out) {
1520  Error("DrawResult()","no such output.");
1521  return;
1522  }
1523  //TCanvas *canvas = new TCanvas("NNresult", "Neural Net output");
1524  if (!opt.Contains("nocanv"))
1525  new TCanvas("NNresult", "Neural Net output");
1526  const Double_t *norm = out->GetNormalisation();
1527  TEventList *events = 0;
1528  TString setname;
1529  Int_t i;
1530  if (opt.Contains("train")) {
1531  events = fTraining;
1532  setname = Form("train%d",index);
1533  } else if (opt.Contains("test")) {
1534  events = fTest;
1535  setname = Form("test%d",index);
1536  }
1537  if ((!fData) || (!events)) {
1538  Error("DrawResult()","no dataset.");
1539  return;
1540  }
1541  if (opt.Contains("comp")) {
1542  //comparison plot
1543  TString title = "Neural Net Output control. ";
1544  title += setname;
1545  setname = "MLP_" + setname + "_comp";
1546  TH2D *hist = ((TH2D *) gDirectory->Get(setname.Data()));
1547  if (!hist)
1548  hist = new TH2D(setname.Data(), title.Data(), 50, -1, 1, 50, -1, 1);
1549  hist->Reset();
1550  Int_t nEvents = events->GetN();
1551  for (i = 0; i < nEvents; i++) {
1552  GetEntry(events->GetEntry(i));
1553  hist->Fill(out->GetValue(), (out->GetBranch() - norm[1]) / norm[0]);
1554  }
1555  hist->Draw();
1556  } else {
1557  //output plot
1558  TString title = "Neural Net Output. ";
1559  title += setname;
1560  setname = "MLP_" + setname;
1561  TH1D *hist = ((TH1D *) gDirectory->Get(setname.Data()));
1562  if (!hist)
1563  hist = new TH1D(setname, title, 50, 1, -1);
1564  hist->Reset();
1565  Int_t nEvents = events->GetN();
1566  for (i = 0; i < nEvents; i++)
1567  hist->Fill(Result(events->GetEntry(i), index));
1568  hist->Draw();
1569  if (opt.Contains("train") && opt.Contains("test")) {
1570  events = fTraining;
1571  setname = "train";
1572  hist = ((TH1D *) gDirectory->Get("MLP_test"));
1573  if (!hist)
1574  hist = new TH1D(setname, title, 50, 1, -1);
1575  hist->Reset();
1576  nEvents = events->GetN();
1577  for (i = 0; i < nEvents; i++)
1578  hist->Fill(Result(events->GetEntry(i), index));
1579  hist->Draw("same");
1580  }
1581  }
1582 }
1583 
1584 ////////////////////////////////////////////////////////////////////////////////
1585 /// Dumps the weights to a text file.
1586 /// Set filename to "-" (default) to dump to the standard output
1587 
1589 {
1590  TString filen = filename;
1591  std::ostream * output;
1592  if (filen == "") {
1593  Error("TMultiLayerPerceptron::DumpWeights()","Invalid file name");
1594  return kFALSE;
1595  }
1596  if (filen == "-")
1597  output = &std::cout;
1598  else
1599  output = new std::ofstream(filen.Data());
1600  TNeuron *neuron = 0;
1601  *output << "#input normalization" << std::endl;
1603  Int_t j=0;
1604  for (j=0;j<nentries;j++) {
1605  neuron = (TNeuron *) fFirstLayer.UncheckedAt(j);
1606  *output << neuron->GetNormalisation()[0] << " "
1607  << neuron->GetNormalisation()[1] << std::endl;
1608  }
1609  *output << "#output normalization" << std::endl;
1611  for (j=0;j<nentries;j++) {
1612  neuron = (TNeuron *) fLastLayer.UncheckedAt(j);
1613  *output << neuron->GetNormalisation()[0] << " "
1614  << neuron->GetNormalisation()[1] << std::endl;
1615  }
1616  *output << "#neurons weights" << std::endl;
1618  while ((neuron = (TNeuron *) it->Next()))
1619  *output << neuron->GetWeight() << std::endl;
1620  delete it;
1622  TSynapse *synapse = 0;
1623  *output << "#synapses weights" << std::endl;
1624  while ((synapse = (TSynapse *) it->Next()))
1625  *output << synapse->GetWeight() << std::endl;
1626  delete it;
1627  if (filen != "-") {
1628  ((std::ofstream *) output)->close();
1629  delete output;
1630  }
1631  return kTRUE;
1632 }
1633 
1634 ////////////////////////////////////////////////////////////////////////////////
1635 /// Loads the weights from a text file conforming to the format
1636 /// defined by DumpWeights.
1637 
1639 {
1640  TString filen = filename;
1641  Double_t w;
1642  if (filen == "") {
1643  Error("TMultiLayerPerceptron::LoadWeights()","Invalid file name");
1644  return kFALSE;
1645  }
1646  char *buff = new char[100];
1647  std::ifstream input(filen.Data());
1648  // input normalzation
1649  input.getline(buff, 100);
1651  Float_t n1,n2;
1652  TNeuron *neuron = 0;
1653  while ((neuron = (TNeuron *) it->Next())) {
1654  input >> n1 >> n2;
1655  neuron->SetNormalisation(n2,n1);
1656  }
1657  input.getline(buff, 100);
1658  // output normalization
1659  input.getline(buff, 100);
1660  delete it;
1662  while ((neuron = (TNeuron *) it->Next())) {
1663  input >> n1 >> n2;
1664  neuron->SetNormalisation(n2,n1);
1665  }
1666  input.getline(buff, 100);
1667  // neuron weights
1668  input.getline(buff, 100);
1669  delete it;
1670  it = (TObjArrayIter *) fNetwork.MakeIterator();
1671  while ((neuron = (TNeuron *) it->Next())) {
1672  input >> w;
1673  neuron->SetWeight(w);
1674  }
1675  delete it;
1676  input.getline(buff, 100);
1677  // synapse weights
1678  input.getline(buff, 100);
1680  TSynapse *synapse = 0;
1681  while ((synapse = (TSynapse *) it->Next())) {
1682  input >> w;
1683  synapse->SetWeight(w);
1684  }
1685  delete it;
1686  delete[] buff;
1687  return kTRUE;
1688 }
1689 
1690 ////////////////////////////////////////////////////////////////////////////////
1691 /// Returns the Neural Net for a given set of input parameters
1692 /// #parameters must equal #input neurons
1693 
1695 {
1697  TNeuron *neuron;
1698  while ((neuron = (TNeuron *) it->Next()))
1699  neuron->SetNewEvent();
1700  delete it;
1702  Int_t i=0;
1703  while ((neuron = (TNeuron *) it->Next()))
1704  neuron->ForceExternalValue(params[i++]);
1705  delete it;
1706  TNeuron *out = (TNeuron *) (fLastLayer.At(index));
1707  if (out)
1708  return out->GetValue();
1709  else
1710  return 0;
1711 }
1712 
1713 ////////////////////////////////////////////////////////////////////////////////
1714 /// Exports the NN as a function for any non-ROOT-dependant code
1715 /// Supported languages are: only C++ , FORTRAN and Python (yet)
1716 /// This feature is also useful if you want to plot the NN as
1717 /// a function (TF1 or TF2).
1718 
1719 void TMultiLayerPerceptron::Export(Option_t * filename, Option_t * language) const
1720 {
1721  TString lg = language;
1722  lg.ToUpper();
1723  Int_t i;
1724  if(GetType()==TNeuron::kExternal) {
1725  Warning("TMultiLayerPerceptron::Export","Request to export a network using an external function");
1726  }
1727  if (lg == "C++") {
1728  TString basefilename = filename;
1729  Int_t slash = basefilename.Last('/')+1;
1730  if (slash) basefilename = TString(basefilename(slash, basefilename.Length()-slash));
1731 
1732  TString classname = basefilename;
1733  TString header = filename;
1734  header += ".h";
1735  TString source = filename;
1736  source += ".cxx";
1737  std::ofstream headerfile(header);
1738  std::ofstream sourcefile(source);
1739  headerfile << "#ifndef " << basefilename << "_h" << std::endl;
1740  headerfile << "#define " << basefilename << "_h" << std::endl << std::endl;
1741  headerfile << "class " << classname << " { " << std::endl;
1742  headerfile << "public:" << std::endl;
1743  headerfile << " " << classname << "() {}" << std::endl;
1744  headerfile << " ~" << classname << "() {}" << std::endl;
1745  sourcefile << "#include \"" << header << "\"" << std::endl;
1746  sourcefile << "#include <cmath>" << std::endl << std::endl;
1747  headerfile << " double Value(int index";
1748  sourcefile << "double " << classname << "::Value(int index";
1749  for (i = 0; i < fFirstLayer.GetEntriesFast(); i++) {
1750  headerfile << ",double in" << i;
1751  sourcefile << ",double in" << i;
1752  }
1753  headerfile << ");" << std::endl;
1754  sourcefile << ") {" << std::endl;
1755  for (i = 0; i < fFirstLayer.GetEntriesFast(); i++)
1756  sourcefile << " input" << i << " = (in" << i << " - "
1757  << ((TNeuron *) fFirstLayer[i])->GetNormalisation()[1] << ")/"
1758  << ((TNeuron *) fFirstLayer[i])->GetNormalisation()[0] << ";"
1759  << std::endl;
1760  sourcefile << " switch(index) {" << std::endl;
1761  TNeuron *neuron;
1763  Int_t idx = 0;
1764  while ((neuron = (TNeuron *) it->Next()))
1765  sourcefile << " case " << idx++ << ":" << std::endl
1766  << " return neuron" << neuron << "();" << std::endl;
1767  sourcefile << " default:" << std::endl
1768  << " return 0.;" << std::endl << " }"
1769  << std::endl;
1770  sourcefile << "}" << std::endl << std::endl;
1771  headerfile << " double Value(int index, double* input);" << std::endl;
1772  sourcefile << "double " << classname << "::Value(int index, double* input) {" << std::endl;
1773  for (i = 0; i < fFirstLayer.GetEntriesFast(); i++)
1774  sourcefile << " input" << i << " = (input[" << i << "] - "
1775  << ((TNeuron *) fFirstLayer[i])->GetNormalisation()[1] << ")/"
1776  << ((TNeuron *) fFirstLayer[i])->GetNormalisation()[0] << ";"
1777  << std::endl;
1778  sourcefile << " switch(index) {" << std::endl;
1779  delete it;
1781  idx = 0;
1782  while ((neuron = (TNeuron *) it->Next()))
1783  sourcefile << " case " << idx++ << ":" << std::endl
1784  << " return neuron" << neuron << "();" << std::endl;
1785  sourcefile << " default:" << std::endl
1786  << " return 0.;" << std::endl << " }"
1787  << std::endl;
1788  sourcefile << "}" << std::endl << std::endl;
1789  headerfile << "private:" << std::endl;
1790  for (i = 0; i < fFirstLayer.GetEntriesFast(); i++)
1791  headerfile << " double input" << i << ";" << std::endl;
1792  delete it;
1793  it = (TObjArrayIter *) fNetwork.MakeIterator();
1794  idx = 0;
1795  while ((neuron = (TNeuron *) it->Next())) {
1796  if (!neuron->GetPre(0)) {
1797  headerfile << " double neuron" << neuron << "();" << std::endl;
1798  sourcefile << "double " << classname << "::neuron" << neuron
1799  << "() {" << std::endl;
1800  sourcefile << " return input" << idx++ << ";" << std::endl;
1801  sourcefile << "}" << std::endl << std::endl;
1802  } else {
1803  headerfile << " double input" << neuron << "();" << std::endl;
1804  sourcefile << "double " << classname << "::input" << neuron
1805  << "() {" << std::endl;
1806  sourcefile << " double input = " << neuron->GetWeight()
1807  << ";" << std::endl;
1808  TSynapse *syn = 0;
1809  Int_t n = 0;
1810  while ((syn = neuron->GetPre(n++))) {
1811  sourcefile << " input += synapse" << syn << "();" << std::endl;
1812  }
1813  sourcefile << " return input;" << std::endl;
1814  sourcefile << "}" << std::endl << std::endl;
1815 
1816  headerfile << " double neuron" << neuron << "();" << std::endl;
1817  sourcefile << "double " << classname << "::neuron" << neuron << "() {" << std::endl;
1818  sourcefile << " double input = input" << neuron << "();" << std::endl;
1819  switch(neuron->GetType()) {
1820  case (TNeuron::kSigmoid):
1821  {
1822  sourcefile << " return ((input < -709. ? 0. : (1/(1+exp(-input)))) * ";
1823  break;
1824  }
1825  case (TNeuron::kLinear):
1826  {
1827  sourcefile << " return (input * ";
1828  break;
1829  }
1830  case (TNeuron::kTanh):
1831  {
1832  sourcefile << " return (tanh(input) * ";
1833  break;
1834  }
1835  case (TNeuron::kGauss):
1836  {
1837  sourcefile << " return (exp(-input*input) * ";
1838  break;
1839  }
1840  case (TNeuron::kSoftmax):
1841  {
1842  sourcefile << " return (exp(input) / (";
1843  Int_t nn = 0;
1844  TNeuron* side = neuron->GetInLayer(nn++);
1845  sourcefile << "exp(input" << side << "())";
1846  while ((side = neuron->GetInLayer(nn++)))
1847  sourcefile << " + exp(input" << side << "())";
1848  sourcefile << ") * ";
1849  break;
1850  }
1851  default:
1852  {
1853  sourcefile << " return (0.0 * ";
1854  }
1855  }
1856  sourcefile << neuron->GetNormalisation()[0] << ")+" ;
1857  sourcefile << neuron->GetNormalisation()[1] << ";" << std::endl;
1858  sourcefile << "}" << std::endl << std::endl;
1859  }
1860  }
1861  delete it;
1862  TSynapse *synapse = 0;
1864  while ((synapse = (TSynapse *) it->Next())) {
1865  headerfile << " double synapse" << synapse << "();" << std::endl;
1866  sourcefile << "double " << classname << "::synapse"
1867  << synapse << "() {" << std::endl;
1868  sourcefile << " return (neuron" << synapse->GetPre()
1869  << "()*" << synapse->GetWeight() << ");" << std::endl;
1870  sourcefile << "}" << std::endl << std::endl;
1871  }
1872  delete it;
1873  headerfile << "};" << std::endl << std::endl;
1874  headerfile << "#endif // " << basefilename << "_h" << std::endl << std::endl;
1875  headerfile.close();
1876  sourcefile.close();
1877  std::cout << header << " and " << source << " created." << std::endl;
1878  }
1879  else if(lg == "FORTRAN") {
1880  TString implicit = " implicit double precision (a-h,n-z)\n";
1881  std::ofstream sigmoid("sigmoid.f");
1882  sigmoid << " double precision FUNCTION SIGMOID(X)" << std::endl
1883  << implicit
1884  << " IF(X.GT.37.) THEN" << std::endl
1885  << " SIGMOID = 1." << std::endl
1886  << " ELSE IF(X.LT.-709.) THEN" << std::endl
1887  << " SIGMOID = 0." << std::endl
1888  << " ELSE" << std::endl
1889  << " SIGMOID = 1./(1.+EXP(-X))" << std::endl
1890  << " ENDIF" << std::endl
1891  << " END" << std::endl;
1892  sigmoid.close();
1893  TString source = filename;
1894  source += ".f";
1895  std::ofstream sourcefile(source);
1896 
1897  // Header
1898  sourcefile << " double precision function " << filename
1899  << "(x, index)" << std::endl;
1900  sourcefile << implicit;
1901  sourcefile << " double precision x(" <<
1902  fFirstLayer.GetEntriesFast() << ")" << std::endl << std::endl;
1903 
1904  // Last layer
1905  sourcefile << "C --- Last Layer" << std::endl;
1906  TNeuron *neuron;
1908  Int_t idx = 0;
1909  TString ifelseif = " if (index.eq.";
1910  while ((neuron = (TNeuron *) it->Next())) {
1911  sourcefile << ifelseif.Data() << idx++ << ") then" << std::endl
1912  << " " << filename
1913  << "=neuron" << neuron << "(x);" << std::endl;
1914  ifelseif = " else if (index.eq.";
1915  }
1916  sourcefile << " else" << std::endl
1917  << " " << filename << "=0.d0" << std::endl
1918  << " endif" << std::endl;
1919  sourcefile << " end" << std::endl;
1920 
1921  // Network
1922  sourcefile << "C --- First and Hidden layers" << std::endl;
1923  delete it;
1924  it = (TObjArrayIter *) fNetwork.MakeIterator();
1925  idx = 0;
1926  while ((neuron = (TNeuron *) it->Next())) {
1927  sourcefile << " double precision function neuron"
1928  << neuron << "(x)" << std::endl
1929  << implicit;
1930  sourcefile << " double precision x("
1931  << fFirstLayer.GetEntriesFast() << ")" << std::endl << std::endl;
1932  if (!neuron->GetPre(0)) {
1933  sourcefile << " neuron" << neuron
1934  << " = (x(" << idx+1 << ") - "
1935  << ((TNeuron *) fFirstLayer[idx])->GetNormalisation()[1]
1936  << "d0)/"
1937  << ((TNeuron *) fFirstLayer[idx])->GetNormalisation()[0]
1938  << "d0" << std::endl;
1939  idx++;
1940  } else {
1941  sourcefile << " neuron" << neuron
1942  << " = " << neuron->GetWeight() << "d0" << std::endl;
1943  TSynapse *syn;
1944  Int_t n = 0;
1945  while ((syn = neuron->GetPre(n++)))
1946  sourcefile << " neuron" << neuron
1947  << " = neuron" << neuron
1948  << " + synapse" << syn << "(x)" << std::endl;
1949  switch(neuron->GetType()) {
1950  case (TNeuron::kSigmoid):
1951  {
1952  sourcefile << " neuron" << neuron
1953  << "= (sigmoid(neuron" << neuron << ")*";
1954  break;
1955  }
1956  case (TNeuron::kLinear):
1957  {
1958  break;
1959  }
1960  case (TNeuron::kTanh):
1961  {
1962  sourcefile << " neuron" << neuron
1963  << "= (tanh(neuron" << neuron << ")*";
1964  break;
1965  }
1966  case (TNeuron::kGauss):
1967  {
1968  sourcefile << " neuron" << neuron
1969  << "= (exp(-neuron" << neuron << "*neuron"
1970  << neuron << "))*";
1971  break;
1972  }
1973  case (TNeuron::kSoftmax):
1974  {
1975  Int_t nn = 0;
1976  TNeuron* side = neuron->GetInLayer(nn++);
1977  sourcefile << " div = exp(neuron" << side << "())" << std::endl;
1978  while ((side = neuron->GetInLayer(nn++)))
1979  sourcefile << " div = div + exp(neuron" << side << "())" << std::endl;
1980  sourcefile << " neuron" << neuron ;
1981  sourcefile << "= (exp(neuron" << neuron << ") / div * ";
1982  break;
1983  }
1984  default:
1985  {
1986  sourcefile << " neuron " << neuron << "= 0.";
1987  }
1988  }
1989  sourcefile << neuron->GetNormalisation()[0] << "d0)+" ;
1990  sourcefile << neuron->GetNormalisation()[1] << "d0" << std::endl;
1991  }
1992  sourcefile << " end" << std::endl;
1993  }
1994  delete it;
1995 
1996  // Synapses
1997  sourcefile << "C --- Synapses" << std::endl;
1998  TSynapse *synapse = 0;
2000  while ((synapse = (TSynapse *) it->Next())) {
2001  sourcefile << " double precision function " << "synapse"
2002  << synapse << "(x)\n" << implicit;
2003  sourcefile << " double precision x("
2004  << fFirstLayer.GetEntriesFast() << ")" << std::endl << std::endl;
2005  sourcefile << " synapse" << synapse
2006  << "=neuron" << synapse->GetPre()
2007  << "(x)*" << synapse->GetWeight() << "d0" << std::endl;
2008  sourcefile << " end" << std::endl << std::endl;
2009  }
2010  delete it;
2011  sourcefile.close();
2012  std::cout << source << " created." << std::endl;
2013  }
2014  else if(lg == "PYTHON") {
2015  TString classname = filename;
2016  TString pyfile = filename;
2017  pyfile += ".py";
2018  std::ofstream pythonfile(pyfile);
2019  pythonfile << "from math import exp" << std::endl << std::endl;
2020  pythonfile << "from math import tanh" << std::endl << std::endl;
2021  pythonfile << "class " << classname << ":" << std::endl;
2022  pythonfile << "\tdef value(self,index";
2023  for (i = 0; i < fFirstLayer.GetEntriesFast(); i++) {
2024  pythonfile << ",in" << i;
2025  }
2026  pythonfile << "):" << std::endl;
2027  for (i = 0; i < fFirstLayer.GetEntriesFast(); i++)
2028  pythonfile << "\t\tself.input" << i << " = (in" << i << " - "
2029  << ((TNeuron *) fFirstLayer[i])->GetNormalisation()[1] << ")/"
2030  << ((TNeuron *) fFirstLayer[i])->GetNormalisation()[0] << std::endl;
2031  TNeuron *neuron;
2033  Int_t idx = 0;
2034  while ((neuron = (TNeuron *) it->Next()))
2035  pythonfile << "\t\tif index==" << idx++
2036  << ": return self.neuron" << neuron << "();" << std::endl;
2037  pythonfile << "\t\treturn 0." << std::endl;
2038  delete it;
2039  it = (TObjArrayIter *) fNetwork.MakeIterator();
2040  idx = 0;
2041  while ((neuron = (TNeuron *) it->Next())) {
2042  pythonfile << "\tdef neuron" << neuron << "(self):" << std::endl;
2043  if (!neuron->GetPre(0))
2044  pythonfile << "\t\treturn self.input" << idx++ << std::endl;
2045  else {
2046  pythonfile << "\t\tinput = " << neuron->GetWeight() << std::endl;
2047  TSynapse *syn;
2048  Int_t n = 0;
2049  while ((syn = neuron->GetPre(n++)))
2050  pythonfile << "\t\tinput = input + self.synapse"
2051  << syn << "()" << std::endl;
2052  switch(neuron->GetType()) {
2053  case (TNeuron::kSigmoid):
2054  {
2055  pythonfile << "\t\tif input<-709. : return " << neuron->GetNormalisation()[1] << std::endl;
2056  pythonfile << "\t\treturn ((1/(1+exp(-input)))*";
2057  break;
2058  }
2059  case (TNeuron::kLinear):
2060  {
2061  pythonfile << "\t\treturn (input*";
2062  break;
2063  }
2064  case (TNeuron::kTanh):
2065  {
2066  pythonfile << "\t\treturn (tanh(input)*";
2067  break;
2068  }
2069  case (TNeuron::kGauss):
2070  {
2071  pythonfile << "\t\treturn (exp(-input*input)*";
2072  break;
2073  }
2074  case (TNeuron::kSoftmax):
2075  {
2076  pythonfile << "\t\treturn (exp(input) / (";
2077  Int_t nn = 0;
2078  TNeuron* side = neuron->GetInLayer(nn++);
2079  pythonfile << "exp(self.neuron" << side << "())";
2080  while ((side = neuron->GetInLayer(nn++)))
2081  pythonfile << " + exp(self.neuron" << side << "())";
2082  pythonfile << ") * ";
2083  break;
2084  }
2085  default:
2086  {
2087  pythonfile << "\t\treturn 0.";
2088  }
2089  }
2090  pythonfile << neuron->GetNormalisation()[0] << ")+" ;
2091  pythonfile << neuron->GetNormalisation()[1] << std::endl;
2092  }
2093  }
2094  delete it;
2095  TSynapse *synapse = 0;
2097  while ((synapse = (TSynapse *) it->Next())) {
2098  pythonfile << "\tdef synapse" << synapse << "(self):" << std::endl;
2099  pythonfile << "\t\treturn (self.neuron" << synapse->GetPre()
2100  << "()*" << synapse->GetWeight() << ")" << std::endl;
2101  }
2102  delete it;
2103  pythonfile.close();
2104  std::cout << pyfile << " created." << std::endl;
2105  }
2106 }
2107 
2108 ////////////////////////////////////////////////////////////////////////////////
2109 /// Shuffle the Int_t index[n] in input.
2110 ///
2111 /// Input:
2112 /// - index: the array to shuffle
2113 /// - n: the size of the array
2114 ///
2115 /// Output:
2116 /// - index: the shuffled indexes
2117 ///
2118 /// This method is used for stochastic training
2119 
2121 {
2122  TTimeStamp ts;
2123  TRandom3 rnd(ts.GetSec());
2124  Int_t j, k;
2125  Int_t a = n - 1;
2126  for (Int_t i = 0; i < n; i++) {
2127  j = (Int_t) (rnd.Rndm() * a);
2128  k = index[j];
2129  index[j] = index[i];
2130  index[i] = k;
2131  }
2132  return;
2133 }
2134 
2135 ////////////////////////////////////////////////////////////////////////////////
2136 /// One step for the stochastic method
2137 /// buffer should contain the previous dw vector and will be updated
2138 
2140 {
2141  Int_t nEvents = fTraining->GetN();
2142  Int_t *index = new Int_t[nEvents];
2143  Int_t i,j,nentries;
2144  for (i = 0; i < nEvents; i++)
2145  index[i] = i;
2146  fEta *= fEtaDecay;
2147  Shuffle(index, nEvents);
2148  TNeuron *neuron;
2149  TSynapse *synapse;
2150  for (i = 0; i < nEvents; i++) {
2151  GetEntry(fTraining->GetEntry(index[i]));
2152  // First compute DeDw for all neurons: force calculation before
2153  // modifying the weights.
2155  for (j=0;j<nentries;j++) {
2156  neuron = (TNeuron *) fFirstLayer.UncheckedAt(j);
2157  neuron->GetDeDw();
2158  }
2159  Int_t cnt = 0;
2160  // Step for all neurons
2162  for (j=0;j<nentries;j++) {
2163  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
2164  buffer[cnt] = (-fEta) * (neuron->GetDeDw() + fDelta)
2165  + fEpsilon * buffer[cnt];
2166  neuron->SetWeight(neuron->GetWeight() + buffer[cnt++]);
2167  }
2168  // Step for all synapses
2170  for (j=0;j<nentries;j++) {
2171  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
2172  buffer[cnt] = (-fEta) * (synapse->GetDeDw() + fDelta)
2173  + fEpsilon * buffer[cnt];
2174  synapse->SetWeight(synapse->GetWeight() + buffer[cnt++]);
2175  }
2176  }
2177  delete[]index;
2178 }
2179 
2180 ////////////////////////////////////////////////////////////////////////////////
2181 /// One step for the batch (stochastic) method.
2182 /// DEDw should have been updated before calling this.
2183 
2185 {
2186  fEta *= fEtaDecay;
2187  Int_t cnt = 0;
2189  TNeuron *neuron = 0;
2190  // Step for all neurons
2191  while ((neuron = (TNeuron *) it->Next())) {
2192  buffer[cnt] = (-fEta) * (neuron->GetDEDw() + fDelta)
2193  + fEpsilon * buffer[cnt];
2194  neuron->SetWeight(neuron->GetWeight() + buffer[cnt++]);
2195  }
2196  delete it;
2198  TSynapse *synapse = 0;
2199  // Step for all synapses
2200  while ((synapse = (TSynapse *) it->Next())) {
2201  buffer[cnt] = (-fEta) * (synapse->GetDEDw() + fDelta)
2202  + fEpsilon * buffer[cnt];
2203  synapse->SetWeight(synapse->GetWeight() + buffer[cnt++]);
2204  }
2205  delete it;
2206 }
2207 
2208 ////////////////////////////////////////////////////////////////////////////////
2209 /// Sets the weights to a point along a line
2210 /// Weights are set to [origin + (dist * dir)].
2211 
2213 {
2214  Int_t idx = 0;
2215  TNeuron *neuron = 0;
2216  TSynapse *synapse = 0;
2218  while ((neuron = (TNeuron *) it->Next())) {
2219  neuron->SetWeight(origin[idx] + (dir[idx] * dist));
2220  idx++;
2221  }
2222  delete it;
2224  while ((synapse = (TSynapse *) it->Next())) {
2225  synapse->SetWeight(origin[idx] + (dir[idx] * dist));
2226  idx++;
2227  }
2228  delete it;
2229 }
2230 
2231 ////////////////////////////////////////////////////////////////////////////////
2232 /// Sets the search direction to steepest descent.
2233 
2235 {
2236  Int_t idx = 0;
2237  TNeuron *neuron = 0;
2238  TSynapse *synapse = 0;
2240  while ((neuron = (TNeuron *) it->Next()))
2241  dir[idx++] = -neuron->GetDEDw();
2242  delete it;
2244  while ((synapse = (TSynapse *) it->Next()))
2245  dir[idx++] = -synapse->GetDEDw();
2246  delete it;
2247 }
2248 
2249 ////////////////////////////////////////////////////////////////////////////////
2250 /// Search along the line defined by direction.
2251 /// buffer is not used but is updated with the new dw
2252 /// so that it can be used by a later stochastic step.
2253 /// It returns true if the line search fails.
2254 
2256 {
2257  Int_t idx = 0;
2258  Int_t j,nentries;
2259  TNeuron *neuron = 0;
2260  TSynapse *synapse = 0;
2261  // store weights before line search
2262  Double_t *origin = new Double_t[fNetwork.GetEntriesFast() +
2265  for (j=0;j<nentries;j++) {
2266  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
2267  origin[idx++] = neuron->GetWeight();
2268  }
2270  for (j=0;j<nentries;j++) {
2271  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
2272  origin[idx++] = synapse->GetWeight();
2273  }
2274  // try to find a triplet (alpha1, alpha2, alpha3) such that
2275  // Error(alpha1)>Error(alpha2)<Error(alpha3)
2276  Double_t err1 = GetError(kTraining);
2277  Double_t alpha1 = 0.;
2278  Double_t alpha2 = fLastAlpha;
2279  if (alpha2 < 0.01)
2280  alpha2 = 0.01;
2281  if (alpha2 > 2.0)
2282  alpha2 = 2.0;
2283  Double_t alpha3 = alpha2;
2284  MLP_Line(origin, direction, alpha2);
2285  Double_t err2 = GetError(kTraining);
2286  Double_t err3 = err2;
2287  Bool_t bingo = false;
2288  Int_t icount;
2289  if (err1 > err2) {
2290  for (icount = 0; icount < 100; icount++) {
2291  alpha3 *= fTau;
2292  MLP_Line(origin, direction, alpha3);
2293  err3 = GetError(kTraining);
2294  if (err3 > err2) {
2295  bingo = true;
2296  break;
2297  }
2298  alpha1 = alpha2;
2299  err1 = err2;
2300  alpha2 = alpha3;
2301  err2 = err3;
2302  }
2303  if (!bingo) {
2304  MLP_Line(origin, direction, 0.);
2305  delete[]origin;
2306  return true;
2307  }
2308  } else {
2309  for (icount = 0; icount < 100; icount++) {
2310  alpha2 /= fTau;
2311  MLP_Line(origin, direction, alpha2);
2312  err2 = GetError(kTraining);
2313  if (err1 > err2) {
2314  bingo = true;
2315  break;
2316  }
2317  alpha3 = alpha2;
2318  err3 = err2;
2319  }
2320  if (!bingo) {
2321  MLP_Line(origin, direction, 0.);
2322  delete[]origin;
2323  fLastAlpha = 0.05;
2324  return true;
2325  }
2326  }
2327  // Sets the weights to the bottom of parabola
2328  fLastAlpha = 0.5 * (alpha1 + alpha3 -
2329  (err3 - err1) / ((err3 - err2) / (alpha3 - alpha2)
2330  - (err2 - err1) / (alpha2 - alpha1)));
2331  fLastAlpha = fLastAlpha < 10000 ? fLastAlpha : 10000;
2332  MLP_Line(origin, direction, fLastAlpha);
2334  // Stores weight changes (can be used by a later stochastic step)
2335  idx = 0;
2337  for (j=0;j<nentries;j++) {
2338  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
2339  buffer[idx] = neuron->GetWeight() - origin[idx];
2340  idx++;
2341  }
2343  for (j=0;j<nentries;j++) {
2344  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
2345  buffer[idx] = synapse->GetWeight() - origin[idx];
2346  idx++;
2347  }
2348  delete[]origin;
2349  return false;
2350 }
2351 
2352 ////////////////////////////////////////////////////////////////////////////////
2353 /// Sets the search direction to conjugate gradient direction
2354 /// beta should be:
2355 ///
2356 /// \f$||g_{(t+1)}||^2 / ||g_{(t)}||^2\f$ (Fletcher-Reeves)
2357 ///
2358 /// \f$g_{(t+1)} (g_{(t+1)}-g_{(t)}) / ||g_{(t)}||^2\f$ (Ribiere-Polak)
2359 
2361 {
2362  Int_t idx = 0;
2363  Int_t j,nentries;
2364  TNeuron *neuron = 0;
2365  TSynapse *synapse = 0;
2367  for (j=0;j<nentries;j++) {
2368  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
2369  dir[idx] = -neuron->GetDEDw() + beta * dir[idx];
2370  idx++;
2371  }
2373  for (j=0;j<nentries;j++) {
2374  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
2375  dir[idx] = -synapse->GetDEDw() + beta * dir[idx];
2376  idx++;
2377  }
2378 }
2379 
2380 ////////////////////////////////////////////////////////////////////////////////
2381 /// Computes the hessian matrix using the BFGS update algorithm.
2382 /// from gamma (g_{(t+1)}-g_{(t)}) and delta (w_{(t+1)}-w_{(t)}).
2383 /// It returns true if such a direction could not be found
2384 /// (if gamma and delta are orthogonal).
2385 
2387 {
2389  if ((Double_t) gd[0][0] == 0.)
2390  return true;
2391  TMatrixD aHg(bfgsh, TMatrixD::kMult, gamma);
2394  Double_t a = 1 / (Double_t) gd[0][0];
2395  Double_t f = 1 + ((Double_t) gHg[0][0] * a);
2396  TMatrixD res( TMatrixD(delta, TMatrixD::kMult,
2397  TMatrixD(TMatrixD::kTransposed, delta)));
2398  res *= f;
2399  res -= (TMatrixD(delta, TMatrixD::kMult, tmp) +
2401  TMatrixD(TMatrixD::kTransposed, delta)));
2402  res *= a;
2403  bfgsh += res;
2404  return false;
2405 }
2406 
2407 ////////////////////////////////////////////////////////////////////////////////
2408 /// Sets the gamma \f$(g_{(t+1)}-g_{(t)})\f$ and delta \f$(w_{(t+1)}-w_{(t)})\f$ vectors
2409 /// Gamma is computed here, so ComputeDEDw cannot have been called before,
2410 /// and delta is a direct translation of buffer into a TMatrixD.
2411 
2413  Double_t * buffer)
2414 {
2416  Int_t idx = 0;
2417  Int_t j,nentries;
2418  TNeuron *neuron = 0;
2419  TSynapse *synapse = 0;
2421  for (j=0;j<nentries;j++) {
2422  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
2423  gamma[idx++][0] = -neuron->GetDEDw();
2424  }
2426  for (j=0;j<nentries;j++) {
2427  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
2428  gamma[idx++][0] = -synapse->GetDEDw();
2429  }
2430  for (Int_t i = 0; i < els; i++)
2431  delta[i].Assign(buffer[i]);
2432  //delta.SetElements(buffer,"F");
2433  ComputeDEDw();
2434  idx = 0;
2436  for (j=0;j<nentries;j++) {
2437  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
2438  gamma[idx++][0] += neuron->GetDEDw();
2439  }
2441  for (j=0;j<nentries;j++) {
2442  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
2443  gamma[idx++][0] += synapse->GetDEDw();
2444  }
2445 }
2446 
2447 ////////////////////////////////////////////////////////////////////////////////
2448 /// scalar product between gradient and direction
2449 /// = derivative along direction
2450 
2452 {
2453  Int_t idx = 0;
2454  Int_t j,nentries;
2455  Double_t output = 0;
2456  TNeuron *neuron = 0;
2457  TSynapse *synapse = 0;
2459  for (j=0;j<nentries;j++) {
2460  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
2461  output += neuron->GetDEDw() * dir[idx++];
2462  }
2464  for (j=0;j<nentries;j++) {
2465  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
2466  output += synapse->GetDEDw() * dir[idx++];
2467  }
2468  return output;
2469 }
2470 
2471 ////////////////////////////////////////////////////////////////////////////////
2472 /// Computes the direction for the BFGS algorithm as the product
2473 /// between the Hessian estimate (bfgsh) and the dir.
2474 
2476 {
2478  TMatrixD dedw(els, 1);
2479  Int_t idx = 0;
2480  Int_t j,nentries;
2481  TNeuron *neuron = 0;
2482  TSynapse *synapse = 0;
2484  for (j=0;j<nentries;j++) {
2485  neuron = (TNeuron *) fNetwork.UncheckedAt(j);
2486  dedw[idx++][0] = neuron->GetDEDw();
2487  }
2489  for (j=0;j<nentries;j++) {
2490  synapse = (TSynapse *) fSynapses.UncheckedAt(j);
2491  dedw[idx++][0] = synapse->GetDEDw();
2492  }
2493  TMatrixD direction(bfgsh, TMatrixD::kMult, dedw);
2494  for (Int_t i = 0; i < els; i++)
2495  dir[i] = -direction[i][0];
2496  //direction.GetElements(dir,"F");
2497 }
2498 
2499 ////////////////////////////////////////////////////////////////////////////////
2500 /// Draws the network structure.
2501 /// Neurons are depicted by a blue disk, and synapses by
2502 /// lines connecting neurons.
2503 /// The line width is proportional to the weight.
2504 
2506 {
2507 #define NeuronSize 2.5
2508 
2509  Int_t nLayers = fStructure.CountChar(':')+1;
2510  Float_t xStep = 1./(nLayers+1.);
2511  Int_t layer;
2512  for(layer=0; layer< nLayers-1; layer++) {
2513  Float_t nNeurons_this = 0;
2514  if(layer==0) {
2515  TString input = TString(fStructure(0, fStructure.First(':')));
2516  nNeurons_this = input.CountChar(',')+1;
2517  }
2518  else {
2519  Int_t cnt=0;
2520  TString hidden = TString(fStructure(fStructure.First(':') + 1,fStructure.Last(':') - fStructure.First(':') - 1));
2521  Int_t beg = 0;
2522  Int_t end = hidden.Index(":", beg + 1);
2523  while (end != -1) {
2524  Int_t num = atoi(TString(hidden(beg, end - beg)).Data());
2525  cnt++;
2526  beg = end + 1;
2527  end = hidden.Index(":", beg + 1);
2528  if(layer==cnt) nNeurons_this = num;
2529  }
2530  Int_t num = atoi(TString(hidden(beg, hidden.Length() - beg)).Data());
2531  cnt++;
2532  if(layer==cnt) nNeurons_this = num;
2533  }
2534  Float_t nNeurons_next = 0;
2535  if(layer==nLayers-2) {
2537  nNeurons_next = output.CountChar(',')+1;
2538  }
2539  else {
2540  Int_t cnt=0;
2541  TString hidden = TString(fStructure(fStructure.First(':') + 1,fStructure.Last(':') - fStructure.First(':') - 1));
2542  Int_t beg = 0;
2543  Int_t end = hidden.Index(":", beg + 1);
2544  while (end != -1) {
2545  Int_t num = atoi(TString(hidden(beg, end - beg)).Data());
2546  cnt++;
2547  beg = end + 1;
2548  end = hidden.Index(":", beg + 1);
2549  if(layer+1==cnt) nNeurons_next = num;
2550  }
2551  Int_t num = atoi(TString(hidden(beg, hidden.Length() - beg)).Data());
2552  cnt++;
2553  if(layer+1==cnt) nNeurons_next = num;
2554  }
2555  Float_t yStep_this = 1./(nNeurons_this+1.);
2556  Float_t yStep_next = 1./(nNeurons_next+1.);
2558  TSynapse *theSynapse = 0;
2559  Float_t maxWeight = 0;
2560  while ((theSynapse = (TSynapse *) it->Next()))
2561  maxWeight = maxWeight < theSynapse->GetWeight() ? theSynapse->GetWeight() : maxWeight;
2562  delete it;
2564  for(Int_t neuron1=0; neuron1<nNeurons_this; neuron1++) {
2565  for(Int_t neuron2=0; neuron2<nNeurons_next; neuron2++) {
2566  TLine* synapse = new TLine(xStep*(layer+1),yStep_this*(neuron1+1),xStep*(layer+2),yStep_next*(neuron2+1));
2567  synapse->Draw();
2568  theSynapse = (TSynapse *) it->Next();
2569  if (!theSynapse) continue;
2570  synapse->SetLineWidth(Int_t((theSynapse->GetWeight()/maxWeight)*10.));
2571  synapse->SetLineStyle(1);
2572  if(((TMath::Abs(theSynapse->GetWeight())/maxWeight)*10.)<0.5) synapse->SetLineStyle(2);
2573  if(((TMath::Abs(theSynapse->GetWeight())/maxWeight)*10.)<0.25) synapse->SetLineStyle(3);
2574  }
2575  }
2576  delete it;
2577  }
2578  for(layer=0; layer< nLayers; layer++) {
2579  Float_t nNeurons = 0;
2580  if(layer==0) {
2581  TString input = TString(fStructure(0, fStructure.First(':')));
2582  nNeurons = input.CountChar(',')+1;
2583  }
2584  else if(layer==nLayers-1) {
2586  nNeurons = output.CountChar(',')+1;
2587  }
2588  else {
2589  Int_t cnt=0;
2590  TString hidden = TString(fStructure(fStructure.First(':') + 1,fStructure.Last(':') - fStructure.First(':') - 1));
2591  Int_t beg = 0;
2592  Int_t end = hidden.Index(":", beg + 1);
2593  while (end != -1) {
2594  Int_t num = atoi(TString(hidden(beg, end - beg)).Data());
2595  cnt++;
2596  beg = end + 1;
2597  end = hidden.Index(":", beg + 1);
2598  if(layer==cnt) nNeurons = num;
2599  }
2600  Int_t num = atoi(TString(hidden(beg, hidden.Length() - beg)).Data());
2601  cnt++;
2602  if(layer==cnt) nNeurons = num;
2603  }
2604  Float_t yStep = 1./(nNeurons+1.);
2605  for(Int_t neuron=0; neuron<nNeurons; neuron++) {
2606  TMarker* m = new TMarker(xStep*(layer+1),yStep*(neuron+1),20);
2607  m->SetMarkerColor(4);
2608  m->SetMarkerSize(NeuronSize);
2609  m->Draw();
2610  }
2611  }
2612  const TString input = TString(fStructure(0, fStructure.First(':')));
2613  const TObjArray *inpL = input.Tokenize(" ,");
2614  const Int_t nrItems = inpL->GetLast()+1;
2615  Float_t yStep = 1./(nrItems+1);
2616  for (Int_t item = 0; item < nrItems; item++) {
2617  const TString brName = ((TObjString *)inpL->At(item))->GetString();
2618  TText* label = new TText(0.5*xStep,yStep*(item+1),brName.Data());
2619  label->Draw();
2620  }
2621  delete inpL;
2622 
2623  Int_t numOutNodes=fLastLayer.GetEntriesFast();
2624  yStep=1./(numOutNodes+1);
2625  for (Int_t outnode=0; outnode<numOutNodes; outnode++) {
2626  TNeuron* neuron=(TNeuron*)fLastLayer[outnode];
2627  if (neuron && neuron->GetName()) {
2628  TText* label = new TText(xStep*nLayers,
2629  yStep*(outnode+1),
2630  neuron->GetName());
2631  label->Draw();
2632  }
2633  }
2634 }
virtual void SetTitleOffset(Float_t offset=1)
Set distance between the axis and the axis title.
Definition: TAttAxis.cxx:294
Double_t DerivDir(Double_t *)
scalar product between gradient and direction = derivative along direction
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
This class describes a neural network.
ENeuronType GetType() const
Returns the neuron type.
Definition: TNeuron.cxx:865
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition: TAttLine.h:43
TObjArray fFirstLayer
Collection of the input neurons; subset of fNetwork.
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3258
Bool_t LineSearch(Double_t *, Double_t *)
Search along the line defined by direction.
void SteepestDir(Double_t *)
Sets the search direction to steepest descent.
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
TSynapse * GetPre(Int_t n) const
Definition: TNeuron.h:35
An array of TObjects.
Definition: TObjArray.h:37
void Randomize() const
Randomize the weights.
void BuildHiddenLayers(TString &)
Builds hidden layers.
virtual Bool_t ProcessEvents()
Process pending events (GUI, timers, sockets).
Definition: TSystem.cxx:410
void BFGSDir(TMatrixD &, Double_t *)
Computes the direction for the BFGS algorithm as the product between the Hessian estimate (bfgsh) and...
Random number generator class based on M.
Definition: TRandom3.h:27
virtual void Remove(TTreeFormula *)
Remove a formula from this manager.
void ForceExternalValue(Double_t value)
Uses the branch type to force an external value.
Definition: TNeuron.cxx:1123
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:854
Double_t fEpsilon
! Epsilon - used in stochastic minimisation - Default=0.
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:23
auto * m
Definition: textangle.C:8
TTreeFormulaManager * fManager
! TTreeFormulaManager for the weight and neurons
Double_t Log(Double_t x)
Definition: TMath.h:750
Collectable string class.
Definition: TObjString.h:28
TTree * fData
! pointer to the tree used as datasource
float Float_t
Definition: RtypesCore.h:53
const char Option_t
Definition: RtypesCore.h:62
const Ssiz_t kNPOS
Definition: RtypesCore.h:111
void SetTestDataSet(TEventList *test)
Sets the Test dataset.
void SetNewEvent() const
Inform the neuron that inputs of the network have changed, so that the buffered values have to be rec...
Definition: TNeuron.cxx:1155
virtual void Update()=0
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:423
Double_t fCurrentTreeWeight
! weight of the current tree in a chain
void SetEpsilon(Double_t eps)
Sets Epsilon - used in stochastic minimisation (look at the constructor for the complete description ...
void Export(Option_t *filename="NNfunction", Option_t *language="C++") const
Exports the NN as a function for any non-ROOT-dependant code Supported languages are: only C++ ...
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition: TMultiGraph.h:35
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1138
Regular expression class.
Definition: TRegexp.h:31
Double_t Evaluate(Int_t index, Double_t *params) const
Returns the Neural Net for a given set of input parameters #parameters must equal #input neurons...
#define R__ASSERT(e)
Definition: TError.h:96
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition: TString.h:634
Definition: test.py:1
ENeuronType
Definition: TNeuron.h:29
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
Read all branches of entry and return total number of bytes read.
Definition: TTree.cxx:5443
void AttachData()
Connects the TTree to Neurons in input and output layers.
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Load a shared library.
Definition: TSystem.cxx:1837
Basic string class.
Definition: TString.h:131
void SetEtaDecay(Double_t ed)
Sets EtaDecay - Eta *= EtaDecay at each epoch (look at the constructor for the complete description o...
Manages Markers.
Definition: TMarker.h:23
#define f(i)
Definition: RSha256.hxx:104
void ToLower()
Change string to lower-case.
Definition: TString.cxx:1125
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
time_t GetSec() const
Definition: TTimeStamp.h:135
Iterator of object array.
Definition: TObjArray.h:123
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:195
void SetEta(Double_t eta)
Sets Eta - used in stochastic minimisation (look at the constructor for the complete description of l...
Bool_t IsNaN(Double_t x)
Definition: TMath.h:882
void SetData(TTree *)
Set the data source.
Double_t GetDeDw() const
Computes the derivative of the error wrt the synapse weight.
Definition: TSynapse.cxx:90
void Shuffle(Int_t *, Int_t) const
Shuffle the Int_t index[n] in input.
TObject * At(Int_t idx) const
Definition: TObjArray.h:166
Double_t GetSumSquareError() const
Error on the output for a given event.
TObjArray fNetwork
Collection of all the neurons in the network.
virtual void Add(TTreeFormula *)
Add a new formula to the list of formulas managed The manager of the formula will be changed and the ...
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
double beta(double x, double y)
Calculates the beta function.
virtual Double_t GetWeight() const
Definition: TTree.h:523
void ConjugateGradientsDir(Double_t *, Double_t)
Sets the search direction to conjugate gradient direction beta should be:
Double_t GetBranch() const
Returns the formula value.
Definition: TNeuron.cxx:912
void DrawResult(Int_t index=0, Option_t *option="test") const
Draws the neural net output It produces an histogram with the output for the two datasets.
Used to coordinate one or more TTreeFormula objects.
TAxis * GetXaxis()
Get x axis of the graph.
bool GetBFGSH(TMatrixD &, TMatrixD &, TMatrixD &)
Computes the hessian matrix using the BFGS update algorithm.
Double_t fEta
! Eta - used in stochastic minimisation - Default=0.1
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t tau
Definition: TRolke.cxx:630
TString fextF
String containing the function name.
Double_t GetValue() const
Computes the output using the appropriate function and all the weighted inputs, or uses the branch as...
Definition: TNeuron.cxx:946
virtual Long64_t GetEntry(Int_t index) const
Return value of entry at index in the list.
Definition: TEventList.cxx:222
Bool_t IsAlnum() const
Returns true if all characters in string are alphanumeric.
Definition: TString.cxx:1746
TTreeFormula * UseBranch(TTree *, const char *)
Sets a formula that can be used to make the neuron an input.
Definition: TNeuron.cxx:876
virtual ~TMultiLayerPerceptron()
Destructor.
void SetReset(Int_t reset)
Sets number of epochs between two resets of the search direction to the steepest descent.
virtual Int_t GetN() const
Definition: TEventList.h:56
TObject * Next()
Return next object in array. Returns 0 when no more objects in array.
Definition: TObjArray.cxx:930
Bool_t EndsWith(const char *pat, ECaseCompare cmp=kExact) const
Return true if string ends with the specified string.
Definition: TString.cxx:2177
static void SetMaxima(Int_t maxop=1000, Int_t maxpar=1000, Int_t maxconst=1000)
static function to set the maximum value of 3 parameters
TNeuron * GetInLayer(Int_t n) const
Definition: TNeuron.h:37
Bool_t DumpWeights(Option_t *filename="-") const
Dumps the weights to a text file.
Base class for several text objects.
Definition: TText.h:23
Int_t fCurrentTree
! index of the current tree in a chain
TNeuron::ENeuronType fType
Type of hidden neurons.
Ssiz_t First(char c) const
Find first occurrence of a character c.
Definition: TString.cxx:499
virtual void Draw(Option_t *chopt="")
Draw this multigraph with its current attributes.
This is a simple weighted bidirectional connection between two neurons.
Definition: TSynapse.h:20
TVirtualPad is an abstract base class for the Pad and Canvas classes.
Definition: TVirtualPad.h:49
virtual void Reset(Option_t *option="")
Reset.
Definition: TH1.cxx:9707
void SetWeight(Double_t w)
Sets the neuron weight to w.
Definition: TNeuron.cxx:1146
TObjArray fSynapses
Collection of all the synapses in the network.
Bool_t IsAlpha() const
Returns true if all characters in string are alphabetic.
Definition: TString.cxx:1731
Used to pass a selection expression to the Tree drawing routine.
Definition: TTreeFormula.h:58
virtual Bool_t Notify()
This method must be overridden to handle object notification.
virtual Int_t GetTreeNumber() const
Definition: TTree.h:498
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:22
const Double_t * GetNormalisation() const
Definition: TNeuron.h:50
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
TEventList * fTraining
! EventList defining the events in the training dataset
static void GetMaxima(Int_t &maxop, Int_t &maxpar, Int_t &maxconst)
static function to get the maximum value of 3 parameters -maxop : maximum number of operations -maxpa...
Int_t GetLast() const
Return index of last object in array.
Definition: TObjArray.cxx:576
#define NeuronSize
double gamma(double x)
void BuildOneHiddenLayer(const TString &sNumNodes, Int_t &layer, Int_t &prevStart, Int_t &prevStop, Bool_t lastLayer)
Builds a hidden layer, updates the number of layers.
void AddInLayer(TNeuron *)
Tells a neuron which neurons form its layer (including itself).
Definition: TNeuron.cxx:855
void SetNormalisation(Double_t mean, Double_t RMS)
Sets the normalization variables.
Definition: TNeuron.cxx:1135
R__EXTERN TSystem * gSystem
Definition: TSystem.h:557
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2981
void Train(Int_t nEpoch, Option_t *option="text", Double_t minE=0)
Train the network.
TNeuron::ENeuronType fOutType
Type of output neurons.
virtual void Draw(Option_t *option="")
Draws the network structure.
auto * a
Definition: textangle.C:12
virtual void UnZoom()
Reset first & last bin to the full range.
Definition: TAxis.cxx:1114
Int_t GetEntriesFast() const
Definition: TObjArray.h:64
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:880
char * Form(const char *fmt,...)
Ssiz_t Length() const
Definition: TString.h:405
Double_t GetError() const
Computes the error for output neurons.
Definition: TNeuron.cxx:1061
A simple line.
Definition: TLine.h:23
TAxis * GetYaxis()
Get y axis of the graph.
void SetDEDw(Double_t in)
Sets the derivative of the total error wrt the neuron weight.
Definition: TNeuron.cxx:1166
Bool_t fTestOwner
! internal flag whether one has to delete fTest or not
Double_t GetError(Int_t event) const
Error on the output for a given event.
A TEventList object is a list of selected events (entries) in a TTree.
Definition: TEventList.h:31
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition: TList.cxx:354
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:614
Bool_t LoadWeights(Option_t *filename="")
Loads the weights from a text file conforming to the format defined by DumpWeights.
Int_t fReset
! number of epochs between two resets of the search direction to the steepest descent - Default=50 ...
const Bool_t kFALSE
Definition: RtypesCore.h:88
void SetTrainingDataSet(TEventList *train)
Sets the Training dataset.
virtual void Modified(Bool_t flag=1)=0
int Ssiz_t
Definition: RtypesCore.h:63
The Canvas class.
Definition: TCanvas.h:31
TObjArray * Tokenize(const TString &delim) const
This function is used to isolate sequential tokens in a TString.
Definition: TString.cxx:2197
TNeuron * GetPre() const
Definition: TSynapse.h:27
TIterator * MakeIterator(Bool_t dir=kIterForward) const
Returns an array iterator.
Definition: TObjArray.cxx:648
void SetWeight(Double_t w)
Sets the weight of the synapse.
Definition: TSynapse.cxx:103
TObject * UncheckedAt(Int_t i) const
Definition: TObjArray.h:90
Double_t GetWeight() const
Definition: TSynapse.h:30
This class describes an elementary neuron, which is the basic element for a Neural Network...
Definition: TNeuron.h:25
#define ClassImp(name)
Definition: Rtypes.h:365
void SetGammaDelta(TMatrixD &, TMatrixD &, Double_t *)
Sets the gamma and delta vectors Gamma is computed here, so ComputeDEDw cannot have been called bef...
void SetDEDw(Double_t in)
Sets the derivative of the total error wrt the synapse weight.
Definition: TSynapse.cxx:111
double Double_t
Definition: RtypesCore.h:55
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:330
Ssiz_t Last(char c) const
Find last occurrence of a character c.
Definition: TString.cxx:892
Double_t GetDeDw() const
Computes the derivative of the error wrt the neuron weight.
Definition: TNeuron.cxx:1082
void ComputeDEDw() const
Compute the DEDw = sum on all training events of dedw for each weight normalized by the number of eve...
Double_t fEtaDecay
! EtaDecay - Eta *= EtaDecay at each epoch - Default=1.
Double_t fDelta
! Delta - used in stochastic minimisation - Default=0.
int type
Definition: TGX11.cxx:120
unsigned long ULong_t
Definition: RtypesCore.h:51
The TTimeStamp encapsulates seconds and ns since EPOCH.
Definition: TTimeStamp.h:71
virtual void Draw(Option_t *opt)
Default Draw method for all objects.
Definition: TTree.h:411
int nentries
Definition: THbookFile.cxx:89
void SetDelta(Double_t delta)
Sets Delta - used in stochastic minimisation (look at the constructor for the complete description of...
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:619
TString fWeight
String containing the event weight.
Double_t fTau
! Tau - used in line search - Default=3.
Double_t GetTarget() const
Computes the normalized target pattern for output neurons.
Definition: TNeuron.cxx:1072
Double_t GetDEDw() const
Definition: TSynapse.h:34
void MLP_Line(Double_t *, Double_t *, Double_t)
Sets the weights to a point along a line Weights are set to [origin + (dist * dir)].
T EvalInstance(Int_t i=0, const char *stringStack[]=0)
Evaluate this treeformula.
static TClass * GetClass(const char *name, Bool_t load=kTRUE, Bool_t silent=kFALSE)
Static method returning pointer to TClass of the specified class name.
Definition: TClass.cxx:2906
TCanvas * slash()
Definition: slash.C:1
TEventList * fTest
! EventList defining the events in the test dataset
virtual void SetLineStyle(Style_t lstyle)
Set the line style.
Definition: TAttLine.h:42
void SetDecimals(Bool_t dot=kTRUE)
Sets the decimals flag By default, blank characters are stripped, and then the label is correctly ali...
Definition: TAxis.h:203
void MLP_Batch(Double_t *)
One step for the batch (stochastic) method.
virtual Long64_t GetEntries() const
Definition: TTree.h:442
Double_t GetWeight() const
Definition: TNeuron.h:48
ELearningMethod fLearningMethod
! The Learning Method
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2257
TList * GetListOfGraphs() const
Definition: TMultiGraph.h:69
Double_t * GetY() const
Definition: TGraph.h:131
TObjArray fLastLayer
Collection of the output neurons; subset of fNetwork.
Double_t GetDEDw() const
Definition: TNeuron.h:53
void MLP_Stochastic(Double_t *)
One step for the stochastic method buffer should contain the previous dw vector and will be updated...
Bool_t fTrainingOwner
! internal flag whether one has to delete fTraining or not
A TGraph is an object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
Double_t GetCrossEntropyBinary() const
Cross entropy error for sigmoid output neurons, for a given event.
TString fStructure
String containing the network structure.
#define gPad
Definition: TVirtualPad.h:286
void ExpandStructure()
Expand the structure of the first layer.
Int_t CountChar(Int_t c) const
Return number of times character c occurs in the string.
Definition: TString.cxx:476
TMultiLayerPerceptron()
Default constructor.
virtual void AddLast(TObject *obj)
Add object in the next empty slot in the array.
Definition: TObjArray.cxx:177
virtual void Reset(Option_t *option="")
Reset this histogram: contents, errors, etc.
Definition: TH2.cxx:3777
A TTree represents a columnar dataset.
Definition: TTree.h:72
#define gDirectory
Definition: TDirectory.h:223
void BuildFirstLayer(TString &)
Instantiates the neurons in input Inputs are normalised and the type is set to kOff (simple forward o...
void SetEventWeight(const char *)
Set the event weight.
void GetEntry(Int_t) const
Load an entry into the network.
virtual void Add(TGraph *graph, Option_t *chopt="")
Add a new graph to the list of graphs.
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
Double_t fLastAlpha
! internal parameter used in line search
Double_t Result(Int_t event, Int_t index=0) const
Computes the output for a given event.
static void output(int code)
Definition: gifencode.c:226
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:292
Double_t GetCrossEntropy() const
Cross entropy error for a softmax output neuron, for a given event.
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
Definition: TNamed.cxx:164
void BuildLastLayer(TString &, Int_t)
Builds the output layer Neurons are linear combinations of input, by default.
void SetLearningMethod(TMultiLayerPerceptron::ELearningMethod method)
Sets the learning method.
const Bool_t kTRUE
Definition: RtypesCore.h:87
void SetTau(Double_t tau)
Sets Tau - used in line search (look at the constructor for the complete description of learning meth...
void BuildNetwork()
Instantiates the network from the description.
const Int_t n
Definition: legend1.C:16
char name[80]
Definition: TGX11.cxx:109
TNeuron::ENeuronType GetType() const
const char * cnt
Definition: TXMLSetup.cxx:74
TString fextD
String containing the derivative name.
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:866
virtual void SetLeftMargin(Float_t leftmargin)
Set Pad left margin in fraction of the pad width.
Definition: TAttPad.cxx:110
TTreeFormula * fEventWeight
! formula representing the event weight
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
const char * Data() const
Definition: TString.h:364