gambit is hosted by Hepforge, IPPP Durham
 GAMBIT  v1.5.0-2191-ga4742ac a Global And Modular Bsm Inference Tool
composite.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
21
22 #include <cmath>
23 #include <vector>
24 #include <set>
25 #include <map>
26 #include <algorithm>
27
29 #include "gambit/ScannerBit/priors_rollcall.hpp"
31
32 namespace Gambit
33 {
34
35  // All priors are transformations which "stretch" one or more random variates
36  // sampled uniformly from the interval [0,1] (or higher dim. equivalent) into
37  // a sample from a different distribution.
38
39  // All priors will be used by pointers to the base class "BasePrior", so they
40  // must inherit from this class. Their constructors can be used to set up
41  // parameters of the transformation they perform, which should itself be
42  // actioned by the "transform" member function
43
44  // Note that before the transformation by these priors, the random number
45  // generation is totally symmetric in all parameters (this is my current
46  // assumption, may need to relax it to accommodate some fancy scanner)
47  // So the way the prior transformation is defined is what really defines which
48  // parameter in the hypercube is which physical parameter.
49
50  // However, this order has to be the order expected by the scanner wrapper of
51  // the loglikelihood function (see ScannerBit/lib/multinest.cpp for example).
52  // Parameter names are provided along with this function so that we can
53  // match them up in the prior correctly. The idea is that the constructors
54  // for the prior objects should be called in such a way as to match the
55  // required parameter order.
56
57  namespace Priors
58  {
59  std::vector<std::string> expand_dots(const std::vector<std::string> &param_names_in)
60  {
61  std::vector<std::string> param_names = param_names_in;
62
63  for (int i = 0, end = param_names.size(); i < end; i++)
64  {
65  if (param_names[i].find("...") != std::string::npos)
66  {
67  auto p_it = param_names.begin() + i;
68  std::string::size_type pos = p_it->find("...");
69  std::string prefix = p_it->substr(0, pos) + "_";
70  std::stringstream ss(p_it->substr(pos+3));
71  int N = 0;
72  if (bool(ss >> N) && N >0)
73  {
74  p_it = param_names.erase(p_it);
75  std::vector<std::string> temps;
76
77  for (int j = 0; j < N; j++)
78  {
79  std::stringstream ss;
80  ss << j;
81  temps.push_back(prefix + ss.str());
82  }
83
84  param_names.insert(p_it, temps.begin(), temps.end());
85  i += N - 1;
86  }
87  }
88  }
89
90  return param_names;
91  }
93  // Combines prior objects together, so that the Scanner can deal with just one object in a standard way.
94  // This is the class to use for setting simple 1D priors (from the library above) on individual parameters.
95  // It also allows for any combination of MD priors to be set on any combination of subspaces of
96  // the full prior.
97
98  // Constructor
99  CompositePrior::CompositePrior(const Options &model_options, const Options &prior_options)
100  {
101  std::map<std::string, std::string> sameMap;
102  std::unordered_map<std::string, std::pair<double, double>> sameMapOptions;
103  std::vector<BasePrior *> phantomPriors;
104  std::set<std::string> needSet;
105
106  // Get model parameters from the inifile
107  std::vector <std::string> modelNames = model_options.getNames();
108  std::sort(modelNames.begin(), modelNames.end());
109
110  //main loop to enter in parameter values
111  for (auto mod_it = modelNames.begin(), mod_end = modelNames.end(); mod_it != mod_end; ++mod_it)
112  {//loop over iniFile models
113  std::string &mod = *mod_it;
114  std::vector <std::string> parameterNames = model_options.getNames(mod);
115  std::sort(parameterNames.begin(), parameterNames.end());
116
117  int default_N = 0, default_i = 0;
118  bool isDefault = false;
119  std::string default_prefix;
120  std::string::size_type par_pos;
121
122  for (auto par_it = parameterNames.begin(), par_end = parameterNames.end(); par_it != par_end; par_it++)
123  {//loop over iniFile parameters
124  Options par_options;
125  std::string par_name;
126
127  par_pos = par_it->find("...");
128  if (!isDefault && par_pos != std::string::npos)
129  {
130  std::stringstream ss(par_it->substr(par_pos + 3));
131  isDefault = true;
132  if (bool(ss >> default_N) && default_N > 0)
133  {
134  isDefault = true;
135  default_prefix = par_it->substr(0, par_pos) + "_";
136  default_i = 0;
137  }
138  }
139
140  if (isDefault)
141  {
142  std::stringstream ss;
143  ss << default_i++;
144  par_name = default_prefix + ss.str();
145  }
146  else
147  {
148  par_name = *par_it;
149  }
150
151  param_names.push_back(mod + std::string("::") + par_name);
152
153  if (model_options.getNode(mod, *par_it).IsSequence() || model_options.getNode(mod, *par_it).IsScalar())
154  {
155  phantomPriors.push_back(new FixedPrior(mod + std::string("::") + par_name, model_options.getNode(mod, *par_it)));
156  }
157  else
158  {
159  par_options = Options(model_options.getNode(mod, *par_it));
160
162  {
163  std::string connectedName = par_options.getValue<std::string>("same_as");
164  std::string::size_type pos = connectedName.rfind("::");
165  if (pos == std::string::npos)
166  {
167  connectedName += std::string("::") + par_name;
168  }
169
170  sameMap[mod + std::string("::") + par_name] = connectedName;
171
172  auto opt = std::pair<double, double>(1.0, 0.0);
174  {
175  opt.first = par_options.getValue<double>("scale");
176  }
177
179  {
180  opt.second = par_options.getValue<double>("shift");
181  }
182
183  sameMapOptions[mod + std::string("::") + par_name] = opt;
184  }
186  {
187  phantomPriors.push_back(new FixedPrior(mod + std::string("::") + par_name, par_options.getNode("fixed_value")));
188  }
189  else
190  {
191  std::string joined_parname = mod + std::string("::") + par_name;
192
194  {
195  Options options = par_options.getOptions();
196  std::string priortype = par_options.getValue<std::string>("prior_type");
197
198  if(priortype == "same_as")
199  {
201  {
202  sameMap[joined_parname] = options.getValue<std::string>("same_as");
203
204  auto opt = std::pair<double, double>(1.0, 0.0);
206  {
207  opt.first = par_options.getValue<double>("scale");
208  }
209
211  {
212  opt.second = par_options.getValue<double>("shift");
213  }
214
215  sameMapOptions[joined_parname] = opt;
216  }
217  else
218  {
219  scan_err << "Same_as prior for parameter \"" << mod << "\" in model \""<< par_name << "\" has no \"same_as\" entry." << scan_end;
220  }
221  }
222  else
223  {
224  if (prior_creators.find(priortype) == prior_creators.end())
225  {
226  scan_err << "Parameter '"<< par_name <<"' of model '" << mod << "' is of type '"<<priortype
227  <<"', but no entry for this type exists in the factory function map.\n" << prior_creators.print() << scan_end;
228  }
229  else
230  {
231  BasePrior *new_prior = prior_creators.at(priortype)(std::vector<std::string>(1, joined_parname),options);
232
233  if (priortype == "fixed")
234  {
235  phantomPriors.push_back(new_prior);
236  }
237  else
238  {
239  my_subpriors.push_back( new_prior );
240  shown_param_names.push_back(joined_parname);
241  }
242  }
243  }
244  }
246  {
247  shown_param_names.push_back(joined_parname);
248
249  my_subpriors.push_back(new RangePrior1D<flatprior>(std::vector<std::string>(1, joined_parname), par_options));
250  }
251  else
252  {
253  shown_param_names.push_back(joined_parname);
254  needSet.insert(joined_parname);
255  }
256  }
257  }
258
259  if (isDefault && default_N > 0)
260  {
261  if (default_N > default_i)
262  par_it--;
263  else
264  isDefault = false;
265  }
266  }
267  }
268
269  // Get the list of priors to build from the iniFile
270  std::vector<std::string> priorNames = prior_options.getNames();
271  std::sort(priorNames.begin(), priorNames.end());
272  std::set<std::string> paramSet(shown_param_names.begin(), shown_param_names.end());
273
274  for (auto priorname_it = priorNames.begin(), priorname_end = priorNames.end(); priorname_it != priorname_end; priorname_it++)
275  {
276  std::string &priorname = *priorname_it;
278  {
279  auto params = expand_dots(prior_options.getValue<std::vector<std::string>>(priorname, "parameters"));
280
281  for (auto par_it = params.begin(), par_end = params.end(); par_it != par_end; par_it++)
282  {
283  if (paramSet.find(*par_it) == paramSet.end())
284  {
285  scan_err << "Parameter " << *par_it << " requested by " << priorname << " is either not defined by the inifile, is fixed, or is the \"same as\" another parameter." << scan_end;
286  }
287  else
288  {
289  auto find_it = needSet.find(*par_it);
290  if (find_it == needSet.end())
291  {
292  scan_err << "Parameter " << *par_it << " requested by prior '"<< priorname <<"' is reserved by a different prior." << scan_end;
293  }
294  else
295  {
296  needSet.erase(find_it);
297  }
298  }
299  }
300
301  auto options = prior_options.getOptions(priorname);
302  auto priortype = prior_options.getValue<std::string>(priorname, "prior_type");
303
304  if (prior_creators.find(priortype) == prior_creators.end())
305  {
306  scan_err << "Prior '"<< priorname <<"' is of type '"<< priortype <<"', but no entry for this type exists in the factory function map.\n" << prior_creators.print() << scan_end;
307  }
308  else
309  {
310  if (priortype == "fixed")
311  {
312  /*for (auto par_it = params.begin(), par_end = params.end(); par_it != par_end; par_it++)
313  {
314  shown_param_names.erase
315  (
316  std::find(shown_param_names.begin(), shown_param_names.end(), *par_it)
317  );
318  }*/
319
320  //my_subpriors.push_back( prior_creators.at(priortype)(params,options) );
321  phantomPriors.push_back( prior_creators.at(priortype)(params,options) );
322  }
323  else if (priortype == "same_as")
324  {
326  {
327  std::string same_name = options.getValue<std::string>("same_as");
328  for (auto par_it = params.begin(), par_end = params.end(); par_it != par_end; ++par_it)
329  {
330  /*shown_param_names.erase
331  (
332  std::find(shown_param_names.begin(), shown_param_names.end(), *par_it)
333  );*/
334  sameMap[*par_it] = same_name;
335  }
336  }
337  else
338  {
339  scan_err << "Same_as prior \"" << priorname << "\" has no \"same_as\" entry." << scan_end;
340  }
341  }
342  else
343  {
344  BasePrior* new_prior = prior_creators.at(priortype)(params,options);
345  my_subpriors.push_back( new_prior );
346
347  auto params = new_prior->getShownParameters();
348  shown_param_names.insert(shown_param_names.end(), params.begin(), params.end());
349
350  /*if (priortype == "composite")
351  {
352  auto composite_new_prior = static_cast<CompositePrior *>(new_prior);
353  auto params = composite_new_prior->getShownParameters();
354  shown_param_names.insert(shown_param_names.end(), params.begin(), params.end());
355  }
356  else
357  {
358  auto params = new_prior->getShownParameters();
359  shown_param_names.insert(shown_param_names.end(), params.begin(), params.end());
360  }*/
361  }
362  }
363  }
364  else
365  {
366  scan_err << "\"parameters\" and \"prior_type\" need to be defined for prior \"" << priorname << "\"" << scan_end;
367  }
368  }
369
370  if (needSet.size() != 0)
371  {
372  scan_err << "Priors are not defined for the following parameters: [";
373  auto it = needSet.begin();
374  scan_err << *(it++);
375  for (; it != needSet.end(); it++)
376  {
377  scan_err << ", "<< *it;
378  }
379  scan_err << "]" << scan_end;
380  }
381
382  std::map<std::string, std::string> keyMap;
383  std::string index, result;
384  unsigned int reps;
385  for (auto strMap = sameMap.begin(), strMap_end = sameMap.end(); strMap != strMap_end; ++strMap)
386  {
387  index = strMap->first;
388  result = strMap->second;
389  reps = 0;
390  while (sameMap.find(result) != sameMap.end())
391  {
392  index = result;
393  result = sameMap[index];
394
395  if (result == strMap->first)
396  {
397  scan_err << "Parameter " << strMap->first << " is \"same as\" itself." << scan_end;
398  break;
399  }
400
401  if (reps > sameMap.size())
402  {
403  scan_err << "Parameter's \"same as\"'s are loop in on each other." << scan_end;
404  break;
405  }
406  reps++;
407  }
408
409  if (keyMap.find(result) == keyMap.end())
410  keyMap[result] = strMap->first + std::string("+") + result;
411  else
412  keyMap[result] = strMap->first + std::string("+") + keyMap[result];
413  }
414
415  for (auto str_it = shown_param_names.begin(), str_end = shown_param_names.end(); str_it != str_end; str_it++)
416  {
417  auto it = keyMap.find(*str_it);
418  if (it != keyMap.end())
419  {
420  *str_it = it->second;
421  }
422  }
423
424  for (auto key_it = keyMap.begin(), key_end = keyMap.end(); key_it != key_end; key_it++)
425  {
426  if (paramSet.find(key_it->first) == paramSet.end())
427  {
428  scan_err << "same_as: " << key_it->first << " is not defined in inifile." << scan_end;
429  }
430  else
431  {
432  phantomPriors.push_back(new MultiPriors(key_it->second, sameMapOptions));
433  }
434  }
435
436  int param_size = 0;
437  for (auto subprior = my_subpriors.begin(), prior_end = my_subpriors.end(); subprior != prior_end; subprior++)
438  {
439  param_size += (*subprior)->size();
440  }
441
442  setSize(param_size);
443
444  my_subpriors.insert(my_subpriors.end(), phantomPriors.begin(), phantomPriors.end());
445  }
446
447  CompositePrior::CompositePrior(const std::vector<std::string> &params_in, const Options &options_in) : BasePrior(params_in)//, shown_param_names(params_in)
448  {
449  std::map<std::string, std::string> sameMap;
450  std::unordered_map<std::string, std::pair<double, double>> sameMapOptions;
451  std::set<std::string> needSet(params_in.begin(), params_in.end());
452  std::set<std::string> paramSet(params_in.begin(), params_in.end());
453  std::vector<BasePrior *> phantomPriors;
454
455  auto priorNames = options_in.getNames();
456  std::sort(priorNames.begin(), priorNames.end());
457
458  for (auto priorname_it = priorNames.begin(), priorname_end = priorNames.end(); priorname_it != priorname_end; ++priorname_it)
459  {
460  std::string &priorname = *priorname_it;
462  {
463  //auto params = options_in.getValue<std::vector<std::string>>(priorname, "parameters");
464
465  auto params = Gambit::Scanner::get_yaml_vector<std::string>(options_in.getNode(priorname, "parameters"));
466
467  for (auto par_it = params.begin(), par_end = params.end(); par_it != par_end; par_it++)
468  {
469  std::string &par = *par_it;
470  if (paramSet.find(par) == paramSet.end())
471  {
472  scan_err << "Parameter " << par << " requested by " << priorname
473  << " is either not defined by the inifile, is fixed, or is the \"same as\" another parameter." << scan_end;
474  }
475  else
476  {
477  auto find_it = needSet.find(par);
478  if (find_it == needSet.end())
479  {
480  scan_err << "Parameter " << par << " requested by prior '"<< priorname
481  <<"' is reserved by a different prior." << scan_end;
482  }
483  else
484  {
485  needSet.erase(find_it);
486  }
487  }
488  }
489
490  auto options = options_in.getOptions(priorname);
491  auto priortype = options_in.getValue<std::string>(priorname, "prior_type");
492
493  if (prior_creators.find(priortype) == prior_creators.end())
494  {
495  scan_err << "Prior '"<< priorname <<"' is of type '"<< priortype
496  <<"', but no entry for this type exists in the factory function map.\n" << prior_creators.print() << scan_end;
497  }
498  else
499  {
500  if (priortype == "fixed")
501  {
502  /*for (auto par_it = params.begin(), par_end = params.end(); par_it != par_end; par_it++)
503  {
504  shown_param_names.erase
505  (
506  std::find(shown_param_names.begin(), shown_param_names.end(), *par_it)
507  );
508  }*/
509
510  phantomPriors.push_back( prior_creators.at(priortype)(params,options) );
511  }
512  else if (priortype == "same_as")
513  {
515  {
516  std::string same_name = options.getValue<std::string>("same_as");
517  for (auto par_it = params.begin(), par_end = params.end(); par_it != par_end; par_it++)
518  {
519  /*shown_param_names.erase
520  (
521  std::find(shown_param_names.begin(), shown_param_names.end(), *par_it)
522  );*/
523  sameMap[*par_it] = same_name;
524  }
525  }
526  else
527  {
528  scan_err << "Same_as prior \"" << priorname << "\" has no \"same_as\" entry." << scan_end;
529  }
530  }
531  else
532  {
533  //my_subpriors.push_back( prior_creators.at(priortype)(params,options) );
534
535  BasePrior* new_prior = prior_creators.at(priortype)(params,options);
536  my_subpriors.push_back( new_prior );
537
538  auto params = new_prior->getShownParameters();
539  shown_param_names.insert(shown_param_names.end(), params.begin(), params.end());
540  }
541  }
542  }
543  else
544  {
545  scan_err << "\"parameters\" and \"prior_type\" need to be defined for prior \"" << priorname << "\"" << scan_end;
546  }
547  }
548
549  if (needSet.size() != 0)
550  {
551  scan_err << "Priors are not defined for the following parameters: [";
552  auto it = needSet.begin();
553  scan_err << *(it++);
554  for (; it != needSet.end(); it++)
555  {
556  scan_err << ", "<< *it;
557  }
558  scan_err << "]" << scan_end;
559  }
560
561  std::map<std::string, std::string> keyMap;
562  std::string index, result;
563  unsigned int reps;
564  for (auto strMap_it = sameMap.begin(), strMap_end = sameMap.end(); strMap_it != strMap_end; strMap_it++)
565  {
566  index = strMap_it->first;
567  result = strMap_it->second;
568  reps = 0;
569  while (sameMap.find(result) != sameMap.end())
570  {
571  index = result;
572  result = sameMap[index];
573
574  if (result == strMap_it->first)
575  {
576  scan_err << "Parameter " << strMap_it->first << " is \"same as\" itself." << scan_end;
577  break;
578  }
579
580  if (reps > sameMap.size())
581  {
582  scan_err << "Parameter's \"same as\"'s are loop in on each other." << scan_end;
583  break;
584  }
585  reps++;
586  }
587
588  if (keyMap.find(result) == keyMap.end())
589  keyMap[result] = strMap_it->first + std::string("+") + result;
590  else
591  keyMap[result] = strMap_it->first + std::string("+") + keyMap[result];
592  }
593
594  for (auto str_it = shown_param_names.begin(), str_end = shown_param_names.end(); str_it != str_end; str_it++)
595  {
596  auto it = keyMap.find(*str_it);
597  if (it != keyMap.end())
598  {
599  *str_it = it->second;
600  }
601  }
602
603  for (auto key_it = keyMap.begin(), key_end = keyMap.end(); key_it != key_end; key_it++)
604  {
605  if (paramSet.find(key_it->first) == paramSet.end())
606  {
607  scan_err << "same_as: " << key_it->first << " is not defined in inifile." << scan_end;
608  }
609  else
610  {
611  phantomPriors.push_back(new MultiPriors(key_it->second, sameMapOptions));
612  }
613  }
614
615  int param_size = 0;
616  for (auto subprior = my_subpriors.begin(), subprior_end = my_subpriors.end(); subprior != subprior_end; subprior++)
617  {
618  param_size += (*subprior)->size();
619  }
620
621  setSize(param_size);
622
623  my_subpriors.insert(my_subpriors.end(), phantomPriors.begin(), phantomPriors.end());
624  }
625  } // end namespace Priors
626 } // end namespace Gambit
627
Abstract base class for priors.
Definition: base_prior.hpp:40
std::vector< std::string > shown_param_names
Definition: composite.hpp:53
void setSize(const unsigned int size)
Definition: base_prior.hpp:79
std::vector< std::string > param_names
Definition: base_prior.hpp:46
const std::vector< str > getNames(const args &... keys) const
Retrieve values from key-value pairs in options node.
reg_elem< create_prior_function > prior_creators
Definition: priors.hpp:43
CompositePrior(const Options &model_options, const Options &prior_options)
Special "build-a-prior" classi.
Definition: composite.cpp:99
Declarations for the YAML options class.
const Options getOptions(const args &... keys) const
Recursive options retrieval.
bool hasKey(const args &... keys) const
Getters for key/value pairs (which is all the options node should contain)
TYPE getValue(const args &... keys) const
#define scan_err
Defined to macros to output errors in the form: scan_err << "error" << scan_end; scan_warn << "warnin...
Utility Functions for the Gambit Scanner.
std::vector< BasePrior * > my_subpriors
Definition: composite.hpp:52
A parameter that is fixed to a different parameter.
Template class for 1d priors which need only a "range" option in their constructor.
Definition: flat_log.hpp:151
virtual std::vector< std::string > getShownParameters() const
Definition: base_prior.hpp:75
A fixed parameter.
#define scan_end
std::vector< std::string > expand_dots(const std::vector< std::string > &param_names_in)
Definition: composite.cpp:59
TODO: see if we can use this one:
Definition: Analysis.hpp:33
A small wrapper object for &#39;options&#39; nodes.
YAML::Node getNode(const args &... keys) const
Retrieve raw YAML node.