gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
great.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
17 
18 #include <vector>
19 #include <limits>
20 #include <cstdio>
21 
27 
28 #include "TGreatModel.h"
29 #include "TGreatManager.h"
30 #include "TGreatMCMCAlgorithmCovariance.h"
31 
32 
33 // Set up the global scan data container
34 namespace Gambit
35 {
36  namespace GreAT
37  {
39  }
40 }
41 
42 
46 
47 scanner_plugin(great, version(1, 0, 0))
48 {
49  // Access standard Gambit things
50  using namespace Gambit;
51  using namespace Gambit::Scanner;
52 
53  // const static PriorTransform prior;
54 
55  // Error thrown if the following entries are not present in the inifile
56  reqd_inifile_entries(); // None at the momement. Might need some later.
57 
58  // Tell CMake to search for the GreAT library.
59  reqd_libraries("great");
60  reqd_headers("fparser.hh", "TGreatModel.h", "TGreatManager.h", "TGreatMCMCAlgorithmCovariance.h", "ROOT");
61 
62  // Code to execute when the plugin is loaded.
64  {
65  const int MPIrank = get_printer().get_stream()->getRank();
66  if (MPIrank == 0) std::cout << "\033[1;31mLoading GreAT plugin for ScannerBit.\033[0m" << std::endl;
67  // Retrieve the external likelihood calculator
68  GreAT::data.likelihood_function = get_purpose(get_inifile_value<std::string>("like"));
69  // Retrieve the external printer
70  GreAT::data.printer = &(get_printer());
71  }
72 
73  int plugin_main(void)
74  {
75 
76  // Run parameters
77  const int nPar = get_dimension(); // Dimensionality of the parameter space
78  const int nTrialLists = get_inifile_value<int> ("nTrialLists", 10); // Number of trial lists (e.g. Markov chains)
79  const int nTrials = get_inifile_value<int> ("nTrials", 20000); // Number of trials (e.g. Markov steps)
80  const int MPIrank = get_printer().get_stream()->getRank(); // MPI rank of this process
81  const bool resume_mode = get_printer().resume_mode(); // Resuming run or not
82  const str outpath = Gambit::Utils::ensure_path_exists(get_inifile_value<std::string>("default_output_path")+"GreAT-native/");
83  GreAT::data.min_logLike = get_inifile_value<double>("likelihood: model_invalid_for_lnlike_below");
84 
85  // Set up output and MultiRun log filenames
86  std::ostringstream ss1, ss2, ss3;
87  ss1 << outpath << "MCMC_" << MPIrank << ".root";
88  ss2 << outpath << "MultiRun.txt";
89  ss3 << outpath << "MultiRun.txt.lock";
90  std::string outputfilename = ss1.str();
91  std::string multifilename = ss2.str();
92  std::string lockfilename = ss3.str();
93 
94  if (resume_mode and MPIrank == 0)
95  {
96  // Clear GreAT lock file
97  remove(lockfilename.c_str());
98  }
99  else
100  {
101  // Wipe entire previous GreAT output if not resuming
103  }
104  #ifdef WITH_MPI
105  GMPI::Comm mpi;
106  mpi.Barrier();
107  #endif
108 
109  // Creating GreAT Model, i.e. parameter space and function to be minimised.
110  // This is new'd because a TGreATManager later takes ownership of it and deletes it internally.
111  TGreatModel* MyModel = new TGreatModel();
112 
113  // Setting up the hypercube parameter space
114  // MyModel->AddParameter(std::string name, std::string unit, double start_value, double min_value, double max_value);
115  std::string x = "";
116  for(int i = 0; i < nPar; i++)
117  {
118  x = "x" + std::to_string(i+1);
119  MyModel->AddParameter(x, "", 0.5, 0., 1.);
120  }
121 
122  // Setting up the logarithmic likelihoodfunction
123  MyModel->SetLogLikelihoodFunction(GreAT::LogLikelihoodFunction);
124 
125  // Setting up the GreAT Manager
126  if (MPIrank == 0) std::cout << "\033[1;31mCreating GreAT Manager\033[0m" << std::endl;
127  // TGreatManager<typename T> MyManager(TGreatModel*);
128  // Using here a multivariate Gaussian distribution (TGreatMCMCAlgorithmCovariance)
129  TGreatManager<TGreatMCMCAlgorithmCovariance> MyManager(MyModel);
130  // Tell the algorithm to use former points to update its prior
131  MyManager.GetAlgorithm()->SetUpdateStatus(true);
132  // Set the output path, file name, and name for the TTree
133  MyManager.SetOutputFileName(outputfilename);
134  MyManager.SetTreeName("mcmc");
135  // Set number of trials (steps) and triallists (chains)
136  MyManager.SetNTrialLists(nTrialLists);
137  MyManager.SetNTrials(nTrials);
138  // Run GreAT
139  if (MPIrank == 0) std::cout << "\033[1;31mRunning GreAT...\033[0m" << std::endl;
140  //MyManager.ActivateMultiRun(multifilename.c_str());
141  MyManager.Run();
142 
143  // Analyse
144  // 1) Fetch the ROOT file
145  TFile *file;
146  file = TFile::Open(MyManager.GetOutputFileName().c_str());
147  TTree *mcmc = (TTree *) file->Get("mcmc");
148 
149  // 2) Define the estimator
150  // TGreatEstimator<typename T>(TTree*)
151  TGreatEstimator<TGreatMCMCAlgorithmCovariance> estimator(mcmc);
152  // Show the scan statistics
153  estimator.ShowStatistics();
154 
155  // Setup auxilliary stream. It is only needed by the master process
156  if(MPIrank == 0)
157  {
158  Gambit::Options ind_samples_options;// = get_inifile_node("aux_printer_ind_samples_options");
159 
160  // Options to desynchronise print streams from the main Gambit iterations. This allows for random access writing, or writing of global scan data.
161  ind_samples_options.setValue("synchronised", false);
162 
163  if (MPIrank == 0) std::cout << "\033[1;31mWriting points...\033[0m" << std::endl;
164  // Initialise auxiliary print stream
165  GreAT::data.printer->new_stream("ind_samples", ind_samples_options);
166 
167  Scanner::printer* ind_samples_printer(GreAT::data.printer->get_stream("ind_samples"));
168  static const int MPIrank = GreAT::data.likelihood_function->getRank();
169 
170  TGreatMCMCSample *prev_sample = estimator.GetFirstIndSample();
171  unsigned int multiplicity = 0;
172 
173  for(TGreatMCMCSample *sample = prev_sample; sample != 0; sample = estimator.GetNextIndSample())
174  {
175  // count samples to get their posterior weight and save them
176  if(prev_sample->GetID() == sample->GetID())
177  ++multiplicity;
178  else
179  {
180  ind_samples_printer->print(multiplicity, "multiplicity", MPIrank, prev_sample->GetID());
181  ind_samples_printer->print(prev_sample->GetID(), "Point ID", MPIrank, prev_sample->GetID());
182  prev_sample = sample;
183  multiplicity = 1;
184  }
185  }
186  // Save the last point
187  ind_samples_printer->print(multiplicity, "multiplicity", MPIrank, prev_sample->GetID());
188  ind_samples_printer->print(prev_sample->GetID(), "Point ID", MPIrank, prev_sample->GetID());
189 
190  // Finish.
191  std::cout << "\033[1;31mGreAT finished successfully!\033[0m" << std::endl;
192  }
193 
194  return 0;
195  }
196 }
197 
198 namespace Gambit
199 {
200  namespace GreAT
201  {
202  double LogLikelihoodFunction(TGreatPoint& point)
203  {
204  // check if point is within the unit cube
205  std::vector<double> parameter_vector = point.GetPoint();
206 
207  bool outside = false;
208  for (unsigned int i=0;i<parameter_vector.size();i++)
209  {
210  if (parameter_vector[i]<0 || parameter_vector[i]>1)
211  {
212  outside = true;
213  }
214  }
215 
216  if (outside)
217  {
218  // at least one dimension is outside the unit cube so return -1e100 for LogLike
219  return GreAT::data.min_logLike;
220  }
221  else
222  {
223  point.SetID(GreAT::data.likelihood_function->getNextPtID()); // Need to use the *next* PtID because PtID will not move on until the likelihood function is called.
224  return GreAT::data.likelihood_function(parameter_vector);
225  }
226  }
227  }
228 }
greatScanData data
Definition: great.cpp:38
double LogLikelihoodFunction(TGreatPoint &point)
Function to be minimised by GreAT.
Definition: great.cpp:202
General small utility classes, typedefs, etc.
scanner_plugin(great, version(1, 0, 0))
================================================= Interface to ScannerBit
Definition: great.cpp:47
Declarations for the YAML options class.
virtual void new_stream(const std::string &, const Options &)=0
Create auxiliary printer object.
#define plugin_constructor
Runs when the plugin is loaded.
#define reqd_inifile_entries(...)
Tells ScannerBit that these tags are required.
#define plugin_main(...)
Declaration of the main function which will be ran by the interface.
EXPORT_SYMBOLS int remove_all_files_in(const str &dirname, bool error_if_absent=true)
Delete all files in a directory (does not act recursively)
ScannerBit interface to GreAT 1.0.0.
Variadic utilty functions.
std::string str
Shorthand for a standard string.
Definition: Analysis.hpp:35
EXPORT_SYMBOLS const str & ensure_path_exists(const str &)
Ensure that a path exists (and then return the path, for chaining purposes)
A simple C++ wrapper for the MPI C bindings.
Scanner::printer_interface * printer
Definition: great.hpp:35
#define reqd_libraries(...)
Tells ScannerBit that these libraries are requested.
#define reqd_headers(...)
Tells ScannerBit that these include files must exist.
Scanner::like_ptr likelihood_function
Definition: great.hpp:34
void setValue(const KEYTYPE &key, const VALTYPE &val)
Basic setter, for adding extra options.
Structure for passing likelihood and printer data through GreAT to the objective function.
Definition: great.hpp:32
TODO: see if we can use this one:
Definition: Analysis.hpp:33
A small wrapper object for &#39;options&#39; nodes.