gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
likelihood_container.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
30 
35 
36 //#define CORE_DEBUG
37 
38 namespace Gambit
39 {
40 
41  // Methods for Likelihood_Container class.
42 
44  Likelihood_Container::Likelihood_Container(const std::map<str, primary_model_functor *> &functorMap,
45  DRes::DependencyResolver &dependencyResolver, IniParser::IniFile &iniFile,
46  const str &purpose, Printers::BaseBasePrinter& printer)
47  : dependencyResolver (dependencyResolver),
48  printer (printer),
49  functorMap (functorMap),
50  min_valid_lnlike (iniFile.getValueOrDef<double>(0.9*std::numeric_limits<double>::lowest(), "likelihood", "model_invalid_for_lnlike_below")),
51  alt_min_valid_lnlike (iniFile.getValueOrDef<double>(0.5*min_valid_lnlike, "likelihood", "model_invalid_for_lnlike_below_alt")),
52  active_min_valid_lnlike (min_valid_lnlike), // can be switched to the alternate value by the scanner
53  print_invalid_points (iniFile.getValueOrDef<bool>(true, "likelihood", "print_invalid_points")),
54  intralooptime_label ("Runtime(ms) intraloop"),
55  interlooptime_label ("Runtime(ms) interloop"),
56  totallooptime_label ("Runtime(ms) totalloop"),
57  /* Note, likelihood container should be constructed after dependency
58  resolution, so that new printer IDs can be safely acquired without
59  risk of collision with graph vertex IDs */
60  intraloopID(Printers::get_main_param_id(intralooptime_label)),
61  interloopID(Printers::get_main_param_id(interlooptime_label)),
62  totalloopID(Printers::get_main_param_id(totallooptime_label)),
63  #ifdef CORE_DEBUG
64  debug (true)
65  #else
66  debug (iniFile.getValueOrDef<bool>(false, "debug") or iniFile.getValueOrDef<bool>(false, "likelihood", "debug"))
67  #endif
68  {
69  // Set the list of valid return types of functions that can be used for 'purpose' by this container class.
70  const std::vector<str> allowed_types_for_purpose = initVector<str>("double", "std::vector<double>", "float", "std::vector<float>");
71  // Find subset of vertices that match requested purpose
72  auto all_vertices = dependencyResolver.getObsLikeOrder();
73  for (auto it = all_vertices.begin(); it != all_vertices.end(); ++it)
74  {
75  if (dependencyResolver.getIniEntry(*it)->purpose == purpose)
76  {
77  return_types[*it] = dependencyResolver.checkTypeMatch(*it, purpose, allowed_types_for_purpose);
78  target_vertices.push_back(std::move(*it));
79  }
80  else
81  {
82  aux_vertices.push_back(std::move(*it));
83  }
84  }
85  }
86 
88  void Likelihood_Container::setParameters (const std::unordered_map<std::string, double> &parameterMap)
89  {
90  // Set up a stream containing the parameter values, for diagnostic output
91  std::ostringstream parstream;
92 
93  // Iterate over the primary_model_parameters functors of all the models being scanned.
94  for (auto act_it = functorMap.begin(), act_end = functorMap.end(); act_it != act_end; act_it++)
95  {
96  parstream << " " << act_it->first << ":" << endl;
97  // Get the names of the parameters for this model.
98  auto paramkeys = act_it->second->getcontentsPtr()->getKeys();
99  // Iterate over the parameters, setting their values in the primary_model_parameters functors from the parameterMap.
100  for (auto par_it = paramkeys.begin(), par_end = paramkeys.end(); par_it != par_end; par_it++)
101  {
102  str key = act_it->first + "::" + *par_it;
103  auto tmp_it = parameterMap.find(key);
104  if(tmp_it == parameterMap.end())
105  {
106  std::ostringstream err;
107  err << "Error! Failed to set parameter '"<<key<<"' following prior transformation! The parameter could not be found in the map returned by the prior. This probably means that the prior you are using contains a bug." << std::endl;
108  err << "The parameters and values that *were* returned by the prior were:" <<std::endl;
109  if(parameterMap.size()==0){ err << "None! Size of map was zero." << std::endl; }
110  else {
111  for (auto par_jt = parameterMap.begin(); par_jt != parameterMap.end(); ++par_jt)
112  {
113  err << par_jt->first << "=" << par_jt->second << std::endl;
114  }
115  }
116  core_error().raise(LOCAL_INFO,err.str());
117  }
118  parstream << " " << *par_it << ": " << tmp_it->second << endl;
119  act_it->second->getcontentsPtr()->setValue(*par_it, tmp_it->second);
120  }
121  }
122 
123  // Notify all exceptions of the values of the parameters for this point.
124  exception::set_parameters("\n\nYAML-ready parameter values at failed point:\n"+parstream.str());
125 
126  // Print out the MPI rank and values of the parameters for this point if in debug mode.
127  if (debug)
128  {
129  #ifdef WITH_MPI
130  GMPI::Comm COMM_WORLD;
131  std::cout << "MPI process rank: "<< COMM_WORLD.Get_rank() << std::endl;
132  #endif
133  cout << parstream.str();
134  logger() << LogTags::core << "\nBeginning computations for parameter point:\n" << parstream.str() << EOM;
135  }
136  // Print the parameter point to the logs, even if not in debug mode
137  //logger() << LogTags::core << "\nBeginning computations for parameter point:\n" << parstream.str() << EOM;
138 
139 
140  }
141 
143  double Likelihood_Container::main(std::unordered_map<std::string, double> &in)
144  {
145  logger() << LogTags::core << LogTags::debug << "Entered Likelihood_Container::main" << EOM;
146 
147  double lnlike = 0;
148  bool point_invalidated = false;
149 
150  // Check for signals from the scanner to switch to an alternate minimum log likelihood value. TODO: could let scanner plugin set the actual value?
151  static bool switch_done(false); // Disable this check once the switch occurs
152  if(not switch_done)
153  {
154  if(check_for_switch_to_alternate_min_LogL())
155  {
156  active_min_valid_lnlike = alt_min_valid_lnlike; // starts off equal to min_valid_lnlike
157  logger() << "Switched to using alt_min_valid_lnlike ("<<alt_min_valid_lnlike<<") instead of original value ("<<min_valid_lnlike<<")" << EOM;
158  switch_done = true;
159  }
160  }
161 
162  // Check for signals to abort run
163  if(signaldata().check_if_shutdown_begun())
164  {
165  tell_scanner_early_shutdown_in_progress(); // e.g. sets 'quit' flag in Diver
166  logger() << "Informed scanner that early shutdown is in progress and it should secure all its output files if possible." << EOM;
167  }
168 
169  // Decide if we need to skip the likelihood calculation due to shutdown procedure
170  if(signaldata().shutdown_begun() and not scanner_can_quit())
171  {
172  // If the scanner does not have a built-in mechanism for halting the scan early, then we will assume
173  // responsiblity for the process and attempt to shut the scan down from our side.
175  lnlike = alt_min_valid_lnlike; // Always use this larger value to avoid scanner deadlocks (e.g. MultiNest refuses to progress without a likelihood above its minimum threshold)
176  point_invalidated = true; // Will prevent this likelihood value from being flagged as 'valid' by the printer
177  logger() << "Shutdown in progess! The scanner is not flagged as being able to shut itself down, so are managing the shutdown from the likelihood container side. Returning min_valid_lnlike to ScannerBit instead of computing likelihood." << EOM;
178  }
179  else // Do the normal likelihood calculation
180  {
181  // If the shutdown has been triggered but the quit flag is present, then we let the likelihood evaluation proceed as normal.
182 
183  bool compute_aux = true;
184 
185  // Set the values of the parameter point in the PrimaryParameters functor, and log them to cout and/or the logs if desired.
186  setParameters(in);
187 
188  // Logger debug output; things labelled 'LogTags::debug' only get logged if the logger::debug or master debug flags are true, not if only 'likelihood::debug' is true.
189  logger() << LogTags::core << LogTags::debug << "Number of target vertices to calculate: " << target_vertices.size() << endl
190  << "Number of auxiliary vertices to calculate: " << aux_vertices.size() << EOM;
191 
192  // Begin timing of total likelihood evaluation
193  std::chrono::time_point<std::chrono::system_clock> startL = std::chrono::system_clock::now();
194 
195  // Compute time since the previous likelihood evaluation ended
196  std::chrono::duration<double> interloop_time = startL - previous_endL;
197 
198  // First work through the target functors, i.e. the ones contributing to the likelihood.
199  for (auto it = target_vertices.begin(), end = target_vertices.end(); it != end; ++it)
200  {
201  // Log the likelihood being tried.
202  str likelihood_tag = "ikelihood contribution from " + dependencyResolver.get_functor(*it)->origin()
203  + "::" + dependencyResolver.get_functor(*it)->name();
204  if (debug) logger() << LogTags::core << "Calculating l" << likelihood_tag << "." << EOM;
205 
206  try
207  {
208  // Set up debug output streams.
209  std::ostringstream debug_to_cout;
210  if (debug) debug_to_cout << " L" << likelihood_tag << ": ";
211 
212  // Calculate the likelihood component.
214 
215  // Switch depending on whether the functor returns floats or doubles and a single likelihood or a vector of them.
216  str rtype = return_types[*it];
217  if (rtype == "double")
218  {
219  double result = dependencyResolver.getObsLike<double>(*it);
220  if (debug) debug_to_cout << result;
221  lnlike += result;
222  }
223  else if (rtype == "std::vector<double>")
224  {
225  std::vector<double> result = dependencyResolver.getObsLike<std::vector<double> >(*it);
226  for (auto jt = result.begin(); jt != result.end(); ++jt)
227  {
228  if (debug) debug_to_cout << *jt << " ";
229  lnlike += *jt;
230  }
231  }
232  else if (rtype == "float")
233  {
234  float result = dependencyResolver.getObsLike<float>(*it);
235  if (debug) debug_to_cout << result;
236  lnlike += result;
237  }
238  else if (rtype == "std::vector<float>")
239  {
240  std::vector<float> result = dependencyResolver.getObsLike<std::vector<float> >(*it);
241  for (auto jt = result.begin(); jt != result.end(); ++jt)
242  {
243  if (debug) debug_to_cout << *jt << " ";
244  lnlike += *jt;
245  }
246  }
247  else core_error().raise(LOCAL_INFO, "Unexpected target functor type.");
248 
249  // Print debug info
250  if (debug) cout << debug_to_cout.str() << endl;
251 
252  // Don't just roll over if it's a NaN, kill the scan and force the developer to fix it.
253  if (Utils::isnan(lnlike))
254  {
255  core_error().raise(LOCAL_INFO, "L" + likelihood_tag + " is NaN!");
256  }
257 
258  // If we've dropped below the likelihood corresponding to effective zero already, skip the rest of the vertices.
260 
261  // Log completion of this likelihood.
262  if (debug) logger() << LogTags::core << "Computed l" << likelihood_tag << "." << EOM;
263  }
264 
265  // Catch points that are invalid, either due to low like or pathology. Skip the rest of the vertices if a point is invalid.
266  catch(invalid_point_exception& e)
267  {
268  logger() << LogTags::core << "Point invalidated by " << e.thrower()->origin() << "::" << e.thrower()->name() << ": " << e.message() << EOM;
270  lnlike = active_min_valid_lnlike;
271  compute_aux = false;
272  point_invalidated = true;
273  // If print_ivalid_points is false disable the printer
275  printer.disable();
276  if (debug) cout << "Point invalid." << endl;
277  break;
278  }
279  }
280 
281 
282  // If none of the likelihood calculations have invalidated the point, calculate the additional auxiliary observables.
283  if (compute_aux)
284  {
285  if (debug) logger() << LogTags::core << "Completed likelihoods. Calculating additional observables." << EOM;
286 
287  for (auto it = aux_vertices.begin(), end = aux_vertices.end(); it != end; ++it)
288  {
289  // Log the observables being tried.
290  str aux_tag = "dditional observable from " + dependencyResolver.get_functor(*it)->origin()
291  + "::" + dependencyResolver.get_functor(*it)->name();
292  if (debug) logger() << LogTags::core << "Calculating a" << aux_tag << "." << EOM;
293 
294  try
295  {
297  if (debug) logger() << LogTags::core << "Computed a" << aux_tag << "." << EOM;
298  }
300  {
301  logger() << LogTags::core << "Additional observable invalidated by " << e.thrower()->origin()
302  << "::" << e.thrower()->name() << ": " << e.message() << EOM;
303  }
304  }
305  }
306 
307  // If the point is invalid and print_invalid_points = false disable the printer, otherwise print vertices
308  if(point_invalidated and !print_invalid_points)
309  printer.disable();
310  else
311  {
312  for (auto it = target_vertices.begin(), end = target_vertices.end(); it != end; ++it)
313  dependencyResolver.printObsLike(*it,getPtID());
314  for (auto it = aux_vertices.begin(), end = aux_vertices.end(); it != end; ++it)
315  dependencyResolver.printObsLike(*it,getPtID());
316  }
317 
318  // End timing of total likelihood evaluation
319  std::chrono::time_point<std::chrono::system_clock> endL = std::chrono::system_clock::now();
320 
321  // Compute time since the previous likelihood evaluation ended
322  // I.e. computing time of this likelihood, plus overhead from previous inter-loop time.
323  std::chrono::duration<double> true_total_loop_time = endL - previous_endL;
324 
325  // Update stored timing information for use in next loop
326  previous_startL = startL;
327  previous_endL = endL;
328 
329  std::chrono::duration<double> runtimeL = endL - startL;
330  typedef std::chrono::milliseconds ms;
331 
332  // Print timing data
334  {
335  int rank = printer.getRank();
336  // Convert time counts to doubles (had weird problem with long long ints on some systems)
337  double d_runtime = std::chrono::duration_cast<ms>(runtimeL).count();
338  double d_interloop = std::chrono::duration_cast<ms>(interloop_time).count();
339  double d_total = std::chrono::duration_cast<ms>(true_total_loop_time).count();
340  printer.print(d_runtime, intralooptime_label, intraloopID, rank, getPtID());
341  printer.print(d_interloop, interlooptime_label, interloopID, rank, getPtID());
342  printer.print(d_total, totallooptime_label, totalloopID, rank, getPtID());
343  }
344 
345  }
346 
347  if (debug) cout << "Total log-likelihood: " << lnlike << endl << endl;
348  logger() << "Total lnL: " << lnlike << EOM;
350 
351  // Disable the printer so that it doesn't try to output the min_valid_lnlike as a valid likelihood value. ScannerBit will re-enable it when needed again.
352  // Disable only for the next print call
353  if(point_invalidated) printer.disable(1);
354 
355  logger() << LogTags::core << LogTags::debug << "Returning control to ScannerBit" << EOM;
356 
357  return lnlike;
358  }
359 
360 
361 }
362 
void print(T const &in, const std::string &label, const int vertexID, const uint rank, const ulong pointID)
DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry DecayTable::Entry double
Printers::BaseBasePrinter & printer
Bound printer object.
std::chrono::time_point< std::chrono::system_clock > previous_startL
Global record of time that last likelihood evaluation began, for computing true total iteration time...
DRes::DependencyResolver & dependencyResolver
Bound dependency resolver object.
const IniParser::ObservableType * getIniEntry(VertexID)
bool printTiming()
Getter for print_timing flag (used by LikelihoodContainer)
std::vector< DRes::VertexID > aux_vertices
Graph vertices corresponding to additional functors not in ObsLike part of yaml file.
Helper functions for dealing with POSIX signals.
#define LOCAL_INFO
Definition: local_info.hpp:34
STL namespace.
std::map< std::string, double > parameterMap
error & core_error()
Core errors.
EXPORT_SYMBOLS int get_main_param_id(const std::string &)
Returns unique positive parameter id; just a thin wrapper for get_param_id.
str origin() const
Getter for the wrapped function&#39;s origin (module or backend name)
Definition: functors.cpp:121
functor * get_functor(VertexID)
Get the functor corresponding to a single VertexID.
void attempt_soft_shutdown()
Perform soft shutdown if processes can be synchronised.
Gambit invalid point exception class.
Definition: exceptions.hpp:226
std::string message()
Retrieve the message that this exception was raised with.
Definition: exceptions.cpp:374
double main(std::unordered_map< std::string, double > &in)
Evaluate total likelihood function.
std::vector< DRes::VertexID > target_vertices
Graph vertices corresponding to functors in the ObsLike section of yaml file.
std::vector< VertexID > getObsLikeOrder()
Retrieve the order in which target vertices are to be evaluated.
functor * thrower()
Retrieve pointer to the functor that threw the invalid point exception.
Definition: exceptions.cpp:410
GAMBIT signal handling functions.
std::chrono::milliseconds ms
std::chrono::time_point< std::chrono::system_clock > previous_endL
Global record of time that last likelihood evaluation ended, for computing intra-iteration overhead t...
str name() const
Getter for the wrapped function&#39;s name.
Definition: functors.cpp:115
double active_min_valid_lnlike
Active value for the minimum log likelihood (one of the above two values, whichever is currently in-u...
Main dependency resolver.
std::map< DRes::VertexID, str > return_types
Map of return types of target functors.
EXPORT_SYMBOLS SignalData & signaldata()
Retrieve global instance of signal handler options struct.
int getRank()
Retrieve/Set MPI rank (setting is useful for e.g. the postprocessor to re-print points from other ran...
const Logging::endofmessage EOM
Explicit const instance of the end of message struct in Gambit namespace.
Definition: logger.hpp:99
void disable(int n=-1)
"Turn off" printer; i.e.
Logging::LogMaster & logger()
Function to retrieve a reference to the Gambit global log object.
Definition: logger.cpp:95
Printers::BaseBasePrinter printer
Type of the printer objects.
double alt_min_valid_lnlike
Alternate value for the minimum log likelihood (scanner can trigger a switch to this in special circu...
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
Likelihood container declarations.
void calcObsLike(VertexID)
Calculate a single target vertex.
void printObsLike(VertexID, const int)
Print a single target vertex.
double min_valid_lnlike
Primary value of the log likelihood at which a point is considered so unlikely that it can be ruled o...
const int intraloopID
printer IDs for timing data
A simple C++ wrapper for the MPI C bindings.
std::map< str, primary_model_functor * > functorMap
Map of parameter names to values.
const str intralooptime_label
Print labels for timing data.
Likelihood_Container(const std::map< str, primary_model_functor *> &functorMap, DRes::DependencyResolver &dependencyResolver, IniParser::IniFile &iniFile, const str &purpose, Printers::BaseBasePrinter &printer)
Constructor.
void invalidatePointAt(VertexID, bool)
Main inifile class.
Definition: yaml_parser.hpp:89
str checkTypeMatch(VertexID, const str &, const std::vector< str > &)
Ensure that the type of a given vertex is equivalent to at least one of a provided list...
void setParameters(const std::unordered_map< std::string, double > &)
Do the prior transformation and populate the parameter map.
TODO: see if we can use this one:
Definition: Analysis.hpp:33
bool print_invalid_points
Switch to print or not print invalid points to the output file.
TYPE getObsLike(VertexID vertex)
Return the result of a functor.
bool debug
Run in likelihood debug mode?