Belle II Software  release-05-01-25
HopfieldNetwork.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010-2017 Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Martin Heck, Nils Braun *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <tracking/trackFindingVXD/trackSetEvaluator/HopfieldNetwork.h>
12 
13 #include <framework/logging/Logger.h>
14 
15 #include <Eigen/Dense>
16 
17 #include <numeric>
18 
19 using namespace Belle2;
20 
22  std::vector<OverlapResolverNodeInfo>& overlapResolverNodeInfos, unsigned short nIterations)
23 {
24  //Start value for neurons if they are compatible.
25  //Each compatible connection activates a node with this value.
26  //As the sum of all the activations shall be less than one, we divide the
27  //activiation by the total number of Nodes.
28  //Incompatible Nodes get later minus one, which counteracts all activation,
29  //if the incompatible Node is active.
30  if (overlapResolverNodeInfos.size() < 2) {
31  B2DEBUG(20, "No reason to doHopfield with less than 2 nodes!");
32  return 0;
33  }
34 
35  const float compatibilityValue = (1.0 - m_omega) / static_cast<float>(overlapResolverNodeInfos.size() - 1);
36 
37  const size_t overlapSize = overlapResolverNodeInfos.size();
38 
39  //Weight matrix; knows compatibility between each possible pair of Nodes
40  Eigen::MatrixXd W(overlapSize, overlapSize);
41  //A): Set all elements to compatible:
42  W.fill(compatibilityValue);
43 
44  //B): Inform weight matrix elements of incompatible neurons:
45  for (const auto& aTC : overlapResolverNodeInfos) {
46  for (unsigned int overlapIndex : aTC.overlaps) {
47  W(aTC.trackIndex, overlapIndex) = -1.0;
48  }
49  }
50 
51 
52  // Neuron values
53  Eigen::VectorXd x(overlapSize);
54  // randomize neuron values for first iteration:
55  for (unsigned int i = 0; i < overlapSize; i++) {
56  x(i) = gRandom->Uniform(1.0); // WARNING: original does Un(0;0.1) not Un(0;1)!
57  }
58 
59  //Store for results from last round:
60  Eigen::VectorXd xOld(overlapSize);
61 
62  //Order of execution for neuron values:
63  std::vector<unsigned short> sequenceVector(overlapSize);
64  //iota fills the vector with 0, 1, 2, ... , (size-1)
65  std::iota(sequenceVector.begin(), sequenceVector.end(), 0);
66 
67  //The following block will be evaluated to empty, if LOG_NO_B2DEBUG is defined:
68  B2DEBUG(100, "sequenceVector with length " << sequenceVector.size());
69  B2DEBUG(100, "Entries are from begin to end:");
70  for (auto && entry : sequenceVector) {
71  B2DEBUG(100, std::to_string(entry) + ", ");
72  }
73 
74  //Store all values of c for protocolling:
75  std::vector<float> cValues(nIterations);
76  //Store for maximum change of weights between iterations.
77  float c = 1.0;
78 
79  //Iterate until change in weights is small:
80  unsigned iIterations = 0;
81 
82  float T = m_T;
83 
84 
85  while (c > m_cmax) {
86  std::shuffle(sequenceVector.begin(), sequenceVector.end(), TRandomWrapper());
87 
88  xOld = x;
89 
90  for (unsigned int i : sequenceVector) {
91  float aTempVal = W.row(i).dot(x);
92  float act = aTempVal + m_omega * overlapResolverNodeInfos[i].qualityIndicator;
93  x(i) = 0.5 * (1. + tanh(act / T));
94  }
95 
96  T = 0.5 * (T + m_Tmin);
97 
98  //Determine maximum change in weights:
99  c = (x - xOld).cwiseAbs().maxCoeff();
100  B2DEBUG(10, "c value is " << c << " at iteration " << iIterations);
101  cValues.at(iIterations) = c;
102 
103  if (iIterations + 1 == nIterations) {
104  B2INFO("Hopfield reached maximum number of iterations without convergence. cValues are:");
105  for (auto && entry : cValues) {
106  B2INFO(std::to_string(entry));
107  }
108  break;
109  }
110  iIterations++;
111  }
112 
113  //Copy Node values into the activity state of the OverlapResolverNodeInfo objects:
114  for (unsigned int i = 0; i < overlapSize; i++) {
115  overlapResolverNodeInfos[i].activityState = x(i);
116  }
117 
118  return iIterations;
119 }
Belle2::HopfieldNetwork::m_cmax
float m_cmax
maximal change of weights between iterations
Definition: HopfieldNetwork.h:76
Belle2::HopfieldNetwork::m_T
float m_T
start temperature of annealing
Definition: HopfieldNetwork.h:74
Belle2::HopfieldNetwork::TRandomWrapper
Wrap TRandom to be useable as a uniform random number generator with std algorithms like std::shuffle...
Definition: HopfieldNetwork.h:80
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::HopfieldNetwork::doHopfield
unsigned short doHopfield(std::vector< OverlapResolverNodeInfo > &overlapResolverNodeInfos, unsigned short nIterations=20)
Performance of the actual algorithm.
Definition: HopfieldNetwork.cc:21
Belle2::HopfieldNetwork::m_omega
float m_omega
tuning parameter of the hopfield network
Definition: HopfieldNetwork.h:73
Belle2::HopfieldNetwork::m_Tmin
float m_Tmin
minimal temperature allowed
Definition: HopfieldNetwork.h:75