gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
twalk.hpp
Go to the documentation of this file.
1 // GAMBIT: Global and Modular BSM Inference Tool
2 // *********************************************
16 #ifndef TWALK_HPP
17 #define TWALK_HPP
18 
19 #include <iostream>
20 #include <fstream>
21 #include <sstream>
22 #include <vector>
23 
24 #include "scanner_plugin.hpp"
25 #include "random_tools.hpp"
26 
27 namespace Gambit
28 {
29 
30  namespace Scanner
31  {
32 
33  inline bool notUnit(const std::vector<double> &in)
34  {
35  for (auto it = in.begin(), end = in.end(); it != end; ++it)
36  {
37  if(*it < 0.0 || *it > 1.0 || *it == -(*it))
38  {
39  return true;
40  }
41  }
42 
43  return false;
44  }
45 
46  template <typename T>
47  inline typename T::iterator::pointer c_ptr(T &it){return &(*it.begin());}
48 
49  inline std::vector<std::vector<double>> calcCov(const std::vector<std::vector<double>> &pts)
50  {
51  static size_t dim = pts[0].size();
52  static size_t N = pts.size();
53 
54  std::vector<std::vector<double>> covar(dim, std::vector<double>(dim, 0.0));
55  std::vector<double> avg(dim, 0.0);
56  size_t i, j;
57 
58  for (auto it = pts.begin(), pt_end = pts.end(); it != pt_end; ++it)
59  {
60  for (i = 0; i < dim; i++)
61  {
62  avg[i] += (*it)[i];
63  for (j = 0; j < dim; j++)
64  {
65  covar[i][j] += (*it)[i]*(*it)[j];
66  }
67  }
68  }
69 
70  for (i = 0; i < dim; i++)
71  {
72  for (j = 0; j < dim; j++)
73  {
74  covar[i][j] = covar[i][j]/N - avg[i]*avg[j]/N/N;
75  }
76  }
77 
78  return covar;
79  }
80 
81  inline std::vector<std::vector<double>> calcIndent (const std::vector<std::vector<double>> &pts)
82  {
83  static size_t dim = pts[0].size();
84 
85  std::vector<std::vector<double>> covar(dim, std::vector<double>(dim, 0.0));
86  std::vector<double> hi(dim, 1.0), low(dim, 0.0);
87  size_t i;
88 
89  for (auto it = pts.begin(), end = pts.end(); it != end; ++it)
90  {
91  for (i = 0; i < dim; i++)
92  {
93  if (hi[i] < (*it)[i])
94  {
95  hi[i] = (*it)[i];
96  }
97 
98  if (low[i] > (*it)[i])
99  {
100  low[i] = (*it)[i];
101  }
102  }
103  }
104 
105  for (i = 0; i < dim; i++)
106  {
107  covar[i][i] = (hi[i] - low[i])*(hi[i] - low[i])/12.0;
108  }
109 
110  return covar;
111  }
112 
113  class RanNumGen
114  {
115  private:
116  std::vector<RandomPlane *> gDev;
117  double div;
118 
119  public:
120  RanNumGen(int proj, int ma, double din, double alim, double alimt, double div, unsigned long seed = 0) : div(div)
121  {
122  gDev.push_back(new RandomPlane(proj, ma, din, alim, alimt, seed++));
123  }
124 
125  double Doub() {return gDev[0]->Doub();}
126 
127  double ExpDev() {return gDev[0]->ExpDev();}
128 
129  #ifdef WITH_MPI
130  double Dev(std::vector<double> &aNext, std::vector<std::vector<double>> &a0, int t, int tt, int freePts, std::vector<int> &tints)
131  #else
132  double Dev(std::vector<double> &aNext, std::vector<std::vector<double>> &a0, int t, int tt)
133  #endif
134  {
135  double ran = gDev[0]->Doub();
136  double logZ;
137  double b0 = div/2.0;
138  double b1 = div;
139  double b2 = (1.0 + div)/2.0;
140 
141  if (ran < b0)
142  {
143  logZ = gDev[0]->WalkDev(&aNext[0], &a0[t][0], &a0[tt][0]);
144  }
145  else if (ran < b1)
146  {
147  logZ = gDev[0]->TransDev(&aNext[0], &a0[t][0], &a0[tt][0]);
148  }
149  else if (ran < b2)
150  {
151  std::vector<std::vector<double>> temp;
152  #ifdef WITH_MPI
153  for (int i = 0, end = freePts; i < end; i++)
154  {
155  temp.push_back(std::vector<double>(a0[tints[i]]));
156  }
157  #else
158  for (int i = 0, end = a0.size(); i < end; i++)
159  {
160  if (i != t)
161  temp.push_back(std::vector<double>(a0[i]));
162  }
163  #endif
164  //if (!gDev[0]->EnterMat(calcCov(temp)))
165  if (true)
166  {
167  #ifdef WITH_MPI
168  int ttt;
169  while(tints[ttt = freePts*gDev[0]->Doub()] == tt);
170 
171  //double ct = gDev[0]->Max(&a0[tt][0], &a0[tints[ttt]][0]);
172  #else
173  int ttt = (a0.size()-2)*gDev[0]->Doub();
174  if (t < tt)
175  {
176  if (ttt >= t) ttt++;
177  if (ttt >= tt) ttt++;
178  }
179  else
180  {
181  if (ttt >= tt) ttt++;
182  if (ttt >= t) ttt++;
183  }
184  //double ct = gDev[0]->Max(&a0[tt][0], &a0[ttt][0]);
185  //std::cout << ct << std::endl;
186  #endif
187  //int ma = a0[0].size();
188  //std::vector<std::vector<double>> cov (ma, std::vector<double>(ma, 0.0));
189  //for (int i = 0; i < ma; i++)
190  // cov[i][i] = ct*ct;
191  //gDev[0]->EnterMat(cov);
192  gDev[0]->HopBlow(&aNext[0], &a0[t][0], &a0[tt][0], &a0[ttt][0]);
193  //gDev[0]->EnterMat(calcIndent(temp));
194  }
195 
196  //gDev[0]->MultiDev(&aNext[0], &a0[t][0]);
197  logZ = 0.0;
198  }
199  else
200  {
201  std::vector<std::vector<double>> temp;
202  #ifdef WITH_MPI
203  for (int i = 0, end = freePts; i < end; i++)
204  {
205  temp.push_back(std::vector<double>(a0[tints[i]]));
206  }
207  #else
208  for (int i = 0, end = a0.size(); i < end; i++)
209  {
210  if (i != t)
211  temp.push_back(std::vector<double>(a0[i]));
212  }
213  #endif
214  //if (!gDev[0]->EnterMat(calcCov(temp)))
215  if (true)
216  {
217  #ifdef WITH_MPI
218  int ttt;
219  while(tints[ttt = freePts*gDev[0]->Doub()] == tt);
220 
221  //double ct = gDev[0]->Max(&a0[tt][0], &a0[tints[ttt]][0]);
222  //std::cout << ct << " " << t << " " << tints[ttt] << std::endl;getchar();
223  #else
224  int ttt = (a0.size()-2)*gDev[0]->Doub();
225  if (t < tt)
226  {
227  if (ttt >= t) ttt++;
228  if (ttt >= tt) ttt++;
229  }
230  else
231  {
232  if (ttt >= tt) ttt++;
233  if (ttt >= t) ttt++;
234  }
235  //double ct = gDev[0]->Max(&a0[tt][0], &a0[ttt][0]);
236  #endif
237  //int ma = a0[0].size();
238  //std::vector<std::vector<double>> cov (ma, std::vector<double>(ma, 0.0));
239  //for (int i = 0; i < ma; i++)
240  // cov[i][i] = ct*ct;
241  //gDev[0]->EnterMat(cov);
242  gDev[0]->HopBlow(&aNext[0], &a0[tt][0], &a0[tt][0], &a0[ttt][0]);
243  //gDev[0]->EnterMat(calcIndent(temp));
244  }
245 
246  //gDev[0]->MultiDev(&aNext[0], &a0[tt][0]);
247  logZ = 0.0;
248  }
249 
250  return logZ;
251  }
252 
254  {
255  for (auto it = gDev.begin(), end = gDev.end(); it != end; ++it)
256  {
257  delete *it;
258  }
259  }
260  };
261 
262  // struct twalk_resume_info
263  // {
264  // bool resume;
265  // std::string file_path;
266  // int N;
267  //
268  // template <typename... T>
269  // inline void operator(int count, T&&... params)
270  // {
271  // if (resume)
272  // {
273  // ifstream in(file_path.c_str(), ios::binary);
274  // }
275  // }
276  // }
277 
278 
279 
280  void TWalk(Gambit::Scanner::like_ptr LogLike,
282  Gambit::Scanner::resume_params_func set_resume_params,
283  const int &dimension,
284  const double &div,
285  const int &proj,
286  const double &din,
287  const double &alim,
288  const double &alimt,
289  const long long &rand,
290  const double &sqrtR,
291  const int &NChains,
292  const bool &hyper_grid,
293  const int &burn_in,
294  const int &save_freq,
295  const double &hrs_max);
296 
297  #endif
298 
299  /*
300  double ran = gDev.Doub();
301  if (ran < b0)
302  {
303  logZ = gDev.WalkDev(c_ptr(aNext), c_ptr(a0[t]), c_ptr(a0[tt]));
304  }
305  else if (ran < b1)
306  {
307  logZ = gDev.TransDev(c_ptr(aNext), c_ptr(a0[t]), c_ptr(a0[tt]));
308  }
309  else if (ran < b2)
310  {
311  std::vector<std::vector<double>> temp = a0;
312  #ifdef WITH_MPI
313  for (int i = 0, end = NThreads - numtasks; i < end; i++)
314  {
315  temp.push_back(a0[tints[i]]);
316  }
317  #else
318  for (int i = 0, end = a0.size(); i < end; i++)
319  {
320  if (i != tt)
321  temp.push_back(a0[i]);
322  }
323  #endif
324  if (!gDev.EnterMat(calcCov(temp)))
325  {
326  gDev.EnterMat(calcIndent(temp));
327  }
328 
329  gDev.MultiDev(c_ptr(aNext), c_ptr(a0[t]));
330  logZ = 0.0;
331  }
332  else
333  {
334  std::vector<std::vector<double>> temp;
335  #ifdef WITH_MPI
336  for (int i = 0, end = NThreads - numtasks; i < end; i++)
337  {
338  temp.push_back(a0[tints[i]]);
339  }
340  #else
341  for (int i = 0, end = a0.size(); i < end; i++)
342  {
343  if (i != tt)
344  temp.push_back(a0[i]);
345  }
346  #endif
347  if (!gDev.EnterMat(calcCov(temp)))
348  {
349  gDev.EnterMat(calcIndent(temp));
350  }
351 
352  gDev.MultiDev(c_ptr(aNext), c_ptr(a0[tt]));
353  logZ = 0.0;
354  }*/
355  }
356 }
void TWalk(Gambit::Scanner::like_ptr LogLike, Gambit::Scanner::printer_interface &printer, Gambit::Scanner::resume_params_func set_resume_params, const int &dimension, const double &div, const int &proj, const double &din, const double &alim, const double &alimt, const long long &rand, const double &sqrtR, const int &NChains, const bool &hyper_grid, const int &burn_in, const int &save_freq, const double &hrs_max)
Definition: twalk.cpp:90
std::vector< std::vector< double > > calcCov(const std::vector< std::vector< double >> &pts)
Definition: twalk.hpp:49
double Dev(std::vector< double > &aNext, std::vector< std::vector< double >> &a0, int t, int tt)
Definition: twalk.hpp:132
class to interface with the plugin manager resume functions.
Definition: plugin_defs.hpp:52
bool notUnit(const std::vector< double > &in)
Definition: twalk.hpp:33
std::vector< RandomPlane * > gDev
Definition: twalk.hpp:116
std::vector< std::vector< double > > calcIndent(const std::vector< std::vector< double >> &pts)
Definition: twalk.hpp:81
RanNumGen(int proj, int ma, double din, double alim, double alimt, double div, unsigned long seed=0)
Definition: twalk.hpp:120
T::iterator::pointer c_ptr(T &it)
Definition: twalk.hpp:47
Manager class for creating printer objects.
declaration for scanner module
TODO: see if we can use this one:
Definition: Analysis.hpp:33
likelihood container for scanner plugins.