gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
MC_convergence.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
21 
22 #include <omp.h>
27 
28 // #define COLLIDERBIT_DEBUG
29 
30 namespace Gambit
31 {
32  namespace ColliderBit
33  {
34 
36  std::map<const MC_convergence_checker* const, bool> MC_convergence_checker::convergence_map;
37 
39  MC_convergence_checker::MC_convergence_checker() : n_threads(omp_get_max_threads()), converged(false)
40  {
41  n_signals = new std::vector<int>[n_threads];
42  convergence_map[this] = false;
43  }
44 
47  {
48  delete[] n_signals;
49  }
50 
53  {
54  clear();
55  set_settings(settings);
56  }
57 
60  {
61  if (omp_get_thread_num() > 0) utils_error().raise(LOCAL_INFO, "Cannot call this function from inside an OpenMP block.");
62  _settings = &settings;
63  }
64 
67  {
68  if (omp_get_thread_num() > 0) utils_error().raise(LOCAL_INFO, "Cannot call this function from inside an OpenMP block.");
69  converged = false;
70  convergence_map[this] = false;
71  for (int i = 0; i != n_threads; ++i)
72  {
73  n_signals[i].clear();
74  }
75  }
76 
77 
80  {
81  // Abort if the analysis container tracked by this object is already fully converged
82  if (converged) return;
83 
84  // Work out the thread number.
85  int my_thread = omp_get_thread_num();
86 
87  // Loop over all the analyses and populate their current signal predictions
88  n_signals[my_thread].clear();
89  for (auto& analysis_pointer_pair : ac.get_current_analyses_map())
90  {
91  // Loop over all the signal regions in this analysis
92  for (auto& sr : analysis_pointer_pair.second->get_results())
93  {
94  // Update the number of accepted events in this signal region
95  n_signals[my_thread].push_back(sr.n_sig_MC);
96  }
97  }
98  }
99 
100 
103  {
104  if (not converged)
105  {
106 
107  int SR_index = -1;
108  // Loop over all analyses
109  bool analysis_converged;
110  bool all_analyses_converged = true; // Will be set to false if any analysis is not converged
111  for (auto& analysis_pointer_pair : ac.get_current_analyses_map())
112  {
113 
114  analysis_converged = false;
115 
116  // Loop over all the signal regions in this analysis
117  bool SR_converged;
118  bool all_SR_converged = true; // Will be set to false if any SR is not converged
119  for (auto& sr : analysis_pointer_pair.second->get_results())
120  {
121  SR_converged = false;
122  SR_index += 1;
123 
124  // Sum signal count across threads
125  int total_counts = 0;
126  for (int j = 0; j != n_threads; j++)
127  {
128  // Tally up the counts across all threads
129  total_counts += n_signals[j][SR_index];
130  }
131 
132  double fractional_stat_uncert = (total_counts == 0 ? 1.0 : 1.0/sqrt(total_counts));
133  double absolute_stat_uncert = total_counts * fractional_stat_uncert;
134  SR_converged = (_settings->stop_at_sys and total_counts > 0 and absolute_stat_uncert <= sr.n_sig_MC_sys) or
135  (fractional_stat_uncert <= _settings->target_stat);
136 
137  if (not SR_converged) all_SR_converged = false;
138 
139  #ifdef COLLIDERBIT_DEBUG
140  cerr << endl;
141  cerr << "DEBUG: SIGNAL REGION " << SR_index << " of " << n_signals[0].size() << endl;
142  cerr << "DEBUG: SR label: " << sr.sr_label << " in analysis " << analysis_pointer_pair.first << endl;
143  cerr << "DEBUG: absolute_stat_uncert vs sys: " << absolute_stat_uncert << " vs " << sr.n_sig_MC_sys << endl;
144  cerr << "DEBUG: fractional_stat_uncert vs target: " << fractional_stat_uncert << " vs " << _settings->target_stat << endl;
145  cerr << "DEBUG: Is this SR done? " << SR_converged << endl;
146  #endif
147 
148  if (SR_converged)
149  {
150  // Shortcut
152  {
153  converged = true;
154  convergence_map[this] = true;
155  return true;
156  }
157 
159  {
160  analysis_converged = true;
161  break; // break signal region loop
162  }
163  }
164  else // SR not converged
165  {
166  // Shortcut
168  {
169  return false;
170  }
171  }
172  } // End loop over SRs
173 
174  if (_settings->all_SR_must_converge) analysis_converged = all_SR_converged;
175 
176  #ifdef COLLIDERBIT_DEBUG
177  cerr << endl;
178  cerr << "DEBUG: Done looping over SRs for analysis " << analysis_pointer_pair.first << endl;
179  cerr << "DEBUG: analysis_converged = " << analysis_converged << endl;
180  #endif
181 
182  if (not analysis_converged) all_analyses_converged = false;
183 
184  // Shortcut
185  if (analysis_converged and not _settings->all_analyses_must_converge)
186  {
187  converged = true;
188  convergence_map[this] = true;
189  return true;
190  }
191  else if (not analysis_converged and _settings->all_analyses_must_converge)
192  {
193  return false;
194  }
195 
196  } // End loop over analyses
197 
198  #ifdef COLLIDERBIT_DEBUG
199  cerr << endl;
200  cerr << "DEBUG: Done looping over analyses in this container" << endl;
201  cerr << "DEBUG: Current variable values:" << endl;
202  cerr << "DEBUG: analysis_converged = " << analysis_converged << endl;
203  cerr << "DEBUG: all-analysis_converged = " << all_analyses_converged << endl;
204  #endif
205 
206  if (not all_analyses_converged) return false;
207  converged = true;
208  convergence_map[this] = true;
209  } // end: if (not converged
210 
211  // Now check if all instances of this class have also set their entry in the convergence map to true,
212  // implying that all analyses in all containers have reached convergence.
214  {
215  for (auto& it : convergence_map)
216  {
217  if (not it.second) return false;
218  }
219  return true;
220  }
221  return true;
222  }
223 
224  }
225 }
const std::map< str, Analysis * > & get_current_analyses_map() const
Get analyses map for the current collider.
EXPORT_SYMBOLS error & utils_error()
Utility errors.
void clear()
Clear all convergence data (for all threads)
void set_settings(const convergence_settings &)
Provide a pointer to the convergence settings.
#define LOCAL_INFO
Definition: local_info.hpp:34
int n_threads
Total number of threads that the checker is configured to deal with.
const convergence_settings * _settings
A pointer to the convergence settings to use.
void init(const convergence_settings &)
Initialise (or re-initialise) the object.
static std::map< const MC_convergence_checker *const, bool > convergence_map
A map containing pointers to all instances of this class.
Class for holding ColliderBit analyses.
std::vector< int > * n_signals
Pointer to an array holding the signal counts on each thread.
ColliderBit Monte Carlo convergence routines.
bool achieved(const AnalysisContainer &ac)
Check if convergence has been achieved across threads, and across all instances of this class...
Type for holding Monte Carlo convergence settings.
Exception objects required for standalone compilation.
void update(const AnalysisContainer &)
Update the convergence data. This is the only routine meant to be called in parallel.
A class for managing collections of Analysis instances.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
Class for ColliderBit analyses.
bool converged
Flag indicating if everything tracked by this instance is converged.