gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
diver.cpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
16 
17 #include <vector>
18 #include <limits>
19 #include <fstream>
20 
26 
30 
31 scanner_plugin(diver, version(1, 0, 4))
32 {
33 
34  // Access Diver stuff and standard Gambit things
35  using namespace Gambit;
36  using namespace Gambit::Diver_1_0_4;
37 
38  // Error thrown if the following entries are not present in the inifile
40 
41  // Tell cmake to search for the diver library.
42  reqd_libraries("diver");
43 
44  // Set up the scan data container
46 
47  // Code to execute when the plugin is loaded.
49  {
50  // Retrieve the external likelihood calculator
51  data.likelihood_function = get_purpose(get_inifile_value<std::string>("like"));
52  if (data.likelihood_function->getRank() == 0) cout << "Loading Diver differential evolution plugin for ScannerBit." << std::endl;
53  // Retrieve the external printer
54  data.printer = &(get_printer());
55  // Do not allow GAMBIT's own likelihood calculator to directly shut down the scan.
56  // Diver will assume responsibility for this process, triggered externally by
57  // the 'plugin_info.early_shutdown_in_progress()' function.
58  data.likelihood_function->disable_external_shutdown();
59  }
60 
61  int plugin_main (void)
62  {
63  // Path to save Diver samples, resume files, etc
64  str defpath = get_inifile_value<str>("default_output_path");
65  str root = Utils::ensure_path_exists(get_inifile_value<str>("path",defpath+"Diver/native"));
66 
67  // Ask the printer if this is a resumed run or not, and check that the necessary files exist if so.
68  bool resume = get_printer().resume_mode();
69  if (resume)
70  {
71  bool good = true;
72  static const std::vector<str> names = initVector<str>(root+".rparam", root+".devo", root+".raw");
73  for (auto it = names.begin(); it != names.end(); ++it)
74  {
75  std::ifstream file(*it);
76  good = good and file.good() and (file.peek() != std::ifstream::traits_type::eof());
77  file.close();
78  }
79  if (not good)
80  {
81  std::ostringstream warning;
82  warning << "Cannot resume previous Diver run because one or all of" << endl;
83  for (auto it = names.begin(); it != names.end(); ++it) warning << " " << *it << endl;
84  warning << "is missing or empty. This is probably because your last run didn't " << endl
85  << "complete even one generation. Diver will start from scratch, " << endl
86  << "as if you had specified -r.";
87  if (data.likelihood_function->getRank() == 0) cout << "WARNING: " << warning.str() << endl;
88  scan_warn << warning.str() << scan_end;
89  resume = false;
90  }
91  }
92 
93  // Retrieve the global option specifying the minimum interesting likelihood.
94  double gl0 = get_inifile_value<double>("likelihood: model_invalid_for_lnlike_below");
95  // Retrieve the global option specifying the likelihood offset to use
96  double offset = get_inifile_value<double>("likelihood: lnlike_offset", 1e-4*gl0);
97  // Make sure the likleihood functor knows to apply the offset internally in ScannerBit
98  data.likelihood_function->setPurposeOffset(offset);
99  // Offset the minimum interesting likelihood by the offset, and flip it to match diver sign convention.
100  gl0 = -1.0 * (gl0 + offset);
101 
102  // Other Diver run parameters
103  int nPar = get_dimension(); // Dimensionality of the parameter space
104  int nDerived = 0; // Number of derived quantities to output (GAMBIT printers handle these).
105  int nDiscrete = get_inifile_value<int> ("nDiscrete", 0); // Number of parameters that are to be treated as discrete
106  bool partitionDiscrete = get_inifile_value<bool> ("partitionDiscrete", false); // Split the population evenly amongst discrete parameters and evolve separately
107  int maxciv = get_inifile_value<int> ("maxciv", 1); // Maximum number of civilisations
108  int maxgen = get_inifile_value<int> ("maxgen", 5000); // Maximum number of generations per civilisation
109  int NP = get_inifile_value<int> ("NP"); // Population size (individuals per generation)
110  double Cr = get_inifile_value<double>("Cr", 0.9); // Crossover factor
111  double lambda = get_inifile_value<double>("lambda", 0.0); // Mixing factor between best and rand/current
112  bool current = get_inifile_value<bool> ("current", false); // Use current vector for mutation
113  bool expon = get_inifile_value<bool> ("expon", false); // Use exponential crossover
114  int bndry = get_inifile_value<int> ("bndry", 3); // Boundary constraint: 1=brick wall, 2=random re-initialization, 3=reflection
115  bool jDE = get_inifile_value<bool> ("jDE", true); // Use self-adaptive choices for rand/1/bin parameters as per Brest et al 2006
116  bool lambdajDE = get_inifile_value<bool> ("lambdajDE", true); // Use self-adaptive rand-to-best/1/bin parameters; based on Brest et al 2006
117  double convthresh = get_inifile_value<double>("convthresh", 1.e-3); // Threshold for gen-level convergence: smoothed fractional improvement in the mean population value
118  int convsteps = get_inifile_value<int> ("convsteps", 10); // Number of steps to smooth over when checking convergence
119  bool removeDuplicates = get_inifile_value<bool> ("removeDuplicates", true); // Weed out duplicate vectors within a single generation
120  bool doBayesian = false; // Calculate approximate log evidence and posterior weightings
121  double maxNodePop = 1.9; // Population at which node is partitioned in binary space partitioning for posterior
122  double Ztolerance = 0.01; // Input tolerance in log-evidence
123  int savecount = get_inifile_value<int> ("savecount", 1); // Save progress every savecount generations
124  bool native_output = get_inifile_value<bool> ("full_native_output", true); // Output .raw file (Diver native sample output format)
125  int init_pop_strategy = get_inifile_value<int> ("init_population_strategy", 2);// Initialisation strategy: 0=one shot, 1=n-shot, 2=n-shot with error if no valid vectors found.
126  int max_ini_attempts = get_inifile_value<int> ("max_initialisation_attempts", 10000); // Maximum number of times to try to find a valid vector for each slot in the initial population.
127  double max_acceptable_value= get_inifile_value<double>("max_acceptable_value",0.9999*gl0); // Maximum function value to accept for the initial generation if init_population_strategy > 0.
128  int verbose = get_inifile_value<int> ("verbosity", 0); // Output verbosity: 0=only error messages, 1=basic info, 2=civ-level info, 3+=population info
129  int seed = get_inifile_value<int> ("seed", -1); // Base seed for random number generation; non-positive means seed from the system clock
130  double (*prior)(const double[], const int, void*&) = NULL; // Pointer to prior function, only used if doBayesian = true.
131  void* context = &data; // Pointer to GAMBIT likelihood function and printers, passed through to objective function.
132 
133  // Copy the contents of root to a char array.
134  char path[root.length()+1];
135  strcpy(path, root.c_str());
136 
137  // Unit cube boundaries
138  double lowerbounds[nPar]; // Lower boundaries of parameter space to scan
139  double upperbounds[nPar]; // Upper boundaries of parameter space to scan
140  for (int i = 0; i < nPar; i++)
141  {
142  lowerbounds[i] = 0.0;
143  upperbounds[i] = 1.0;
144  }
145 
146  // Scale factors
147  std::vector<double> Fvec = get_inifile_value<std::vector<double> >("F", initVector<double>(0.7));
148  int nF = Fvec.size(); // Size of the array indicating scale factors
149  double F[nF]; // Scale factor(s).
150  std::copy(Fvec.begin(), Fvec.end(), F);
151 
152  // Discrete parameters
153  int discrete[nDiscrete]; // Indices of discrete parameters, Fortran style, i.e. starting at 1!!
154  for (int i = 0; i < nDiscrete; i++)
155  {
156  discrete[i] = 0; //TODO Needs to be set automatically somehow? Not yet sure how to deal with discrete parameters in GAMBIT.
157  }
158 
159  // Run Diver
160  if (data.likelihood_function->getRank() == 0) cout << "Starting Diver run..." << std::endl;
161  cdiver(&objective, nPar, lowerbounds, upperbounds, path, nDerived, nDiscrete,
162  discrete, partitionDiscrete, maxciv, maxgen, NP, nF, F, Cr, lambda, current,
163  expon, bndry, jDE, lambdajDE, convthresh, convsteps, removeDuplicates, doBayesian,
164  prior, maxNodePop, Ztolerance, savecount, resume, native_output, init_pop_strategy,
165  max_ini_attempts, max_acceptable_value, seed, context, verbose);
166  if (data.likelihood_function->getRank() == 0) cout << "Diver run finished!" << std::endl;
167  return 0;
168 
169  }
170 
171 }
172 
176 
177 namespace Gambit
178 {
179 
180  namespace Diver_1_0_4
181  {
182 
183  //Function to be minimized. Corresponds to -ln(Likelihood). Redirects to the target of context pointer.
184  double objective(double params[], const int param_dim, int &fcall, bool &quit, const bool validvector, void*& context)
185  {
186  // Return the worst possible likelihood if the point is outside the prior box.
187  if (not validvector) return std::numeric_limits<double>::max();
188 
189  // Put the parameters into a C++ vector
190  std::vector<double> param_vec(params, params + param_dim);
191 
192  // Retrieve the likelihood function from the context pointer and call it
193  diverScanData* data = static_cast<diverScanData*>(context);
194  double lnlike = data->likelihood_function(param_vec);
195 
196  // Increment the number of function calls, tell Diver to continue and return the likelihood
197  fcall += 1;
198 
199  // Check whether the calling code wants us to shut down early
201 
202  return -lnlike;
203 
204  }
205 
206  }
207 
208 }
209 
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
greatScanData data
Definition: great.cpp:38
double objective(double params[], const int param_dim, int &fcall, bool &quit, const bool validvector, void *&context)
Function to be minimised by Diver.
Definition: diver.cpp:184
#define scan_warn
General small utility classes, typedefs, etc.
ScannerBit interface to Diver 1.0.2.
Scanner::printer_interface * printer
Definition: diver.hpp:40
Declarations for the YAML options class.
GAMBIT warning class.
Definition: exceptions.hpp:165
General small utility functions.
#define plugin_constructor
Runs when the plugin is loaded.
Scanner::like_ptr likelihood_function
Definition: diver.hpp:39
#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.
scanner_plugin(diver, version(1, 0, 0))
================================================= Interface to ScannerBit
Definition: diver.cpp:31
const bool verbose
Definition: logging.cpp:52
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)
DS5_MSPCTM DS_INTDOF int
double lambda(double x, double y, double z)
Definition: MSSM_H.hpp:38
Structure for passing likelihood and printer data through Diver to the objective function.
Definition: diver.hpp:37
#define scan_end
#define reqd_libraries(...)
Tells ScannerBit that these libraries are requested.
void cdiver(double(*)(double[], const int, int &, bool &, const bool, void *&), int, const double[], const double[], const char[], int, int, const int[], bool, const int, const int, int, int, const double[], double, double, bool, bool, int, bool, bool, double, int, bool, bool, double(*)(const double[], const int, void *&), double, double, int, bool, bool, int, int, double, void *&, int)
EXPORT_SYMBOLS pluginInfo plugin_info
Access Functor for plugin info.
TODO: see if we can use this one:
Definition: Analysis.hpp:33