gambit is hosted by Hepforge, IPPP Durham
GAMBIT  v1.5.0-252-gf9a3f78
a Global And Modular Bsm Inference Tool
random_tools.hpp
Go to the documentation of this file.
1 #ifndef RANDOM_TOOLS_HPP
2 #define RANDOM_TOOLS_HPP
3 #include <iostream>
4 #include <cstdio>
5 
7 
8 template <typename T>
9 T *matrix(const int xN)
10 {
11  T *temp = new T[xN];
12 
13  return temp;
14 }
15 
16 template <typename T>
17 T **matrix(const int xN, const int yN)
18 {
19  T **temp = new T*[xN];
20 
21  for (int i = 0; i < xN; i++)
22  temp[i] = new T[yN];
23 
24  return temp;
25 }
26 
27 template <typename T>
28 T ***matrix(const int xN, const int yN, const int zN)
29 {
30  T ***temp = new T**[xN];
31 
32  for (int i = 0; i < xN; i++)
33  {
34  temp[i] = new T*[yN];
35  for (int j = 0; j < yN; j++)
36  temp[i][j] = new T[zN];
37  }
38  return temp;
39 }
40 
41 template <typename T>
42 T *matrix(const int xN, T in)
43 {
44  T *temp = new T[xN];
45 
46  for (int i = 0; i < xN; i++)
47  temp[i] = in;
48 
49  return temp;
50 }
51 
52 template <typename T>
53 T **matrix(const int xN, const int yN, T in)
54 {
55  T **temp = new T*[xN];
56 
57  for (int i = 0; i < xN; i++)
58  {
59  temp[i] = new T[yN];
60  for (int j = 0; j < yN; j++)
61  temp[i][j] = in;
62  }
63 
64  return temp;
65 }
66 
67 template <typename T>
68 T ***matrix(const int xN, const int yN, const int zN, T in)
69 {
70  T ***temp = new T**[xN];
71 
72  for (int i = 0; i < xN; i++)
73  {
74  temp[i] = new T*[yN];
75  for (int j = 0; j < yN; j++)
76  {
77  temp[i][j] = new T[zN];
78  for (int k = 0; k < zN; k++)
79  temp[i][j][k] = in;
80  }
81  }
82  return temp;
83 }
84 
85 template <typename T>
86  T *matrix(const int xN, T *in)
87 {
88  T *temp = new T[xN];
89  for (int i = 0; i < xN; i++)
90  temp[i] = in[i];
91 
92  return temp;
93 }
94 
95 template <typename T>
96  T **matrix(const int xN, const int yN, T **in)
97 {
98  T **temp = new T*[xN];
99 
100  for (int i = 0; i < xN; i++)
101  {
102  temp[i] = new T[yN];
103  for (int j = 0; j < yN; j++)
104  temp[i][j] = in[i][j];
105  }
106 
107  return temp;
108 }
109 
110 template <typename T>
111  T ***matrix(const int xN, const int yN, const int zN, T ***in)
112 {
113  T ***temp = new T**[xN];
114 
115  for (int i = 0; i < xN; i++)
116  {
117  temp[i] = new T*[yN];
118  for (int j = 0; j < yN; j++)
119  {
120  temp[i][j] = new T[zN];
121  for (int k = 0; k < zN; k++)
122  temp[i][j][k] = in[i][j][k];
123  }
124  }
125  return temp;
126 }
127 
128 
129 template <typename T>
130 void del(T *temp)
131 {
132  delete[] temp;
133  temp = NULL;
134 }
135 
136 template <typename T>
137 void del(T **temp, int xN)
138 {
139  for (int i = 0; i < xN; i++)
140  delete[] temp[i];
141 
142  delete[] temp;
143  temp = NULL;
144 }
145 
146 template <typename T>
147 void del(T ***temp, int xN, int yN)
148 {
149  for (int i = 0; i < xN; i++)
150  {
151  for (int j = 0; j < yN; j++)
152  delete[] temp[i][j];
153  delete[] temp[i];
154  }
155 
156  delete[] temp;
157  temp = NULL;
158 }
159 
160 template<class T>
161 inline const T SQ(const T a) {return a*a;}
162 
163 template<class T>
164 inline const T SQR(const T a) {return a*a;}
165 
166 class Cholesky
167 {
168  private:
169  double **el;
170 
171  protected:
172  int num;
173 
174 
175  public:
176  Cholesky(const int nin) : num(nin)
177  {
178  el = matrix <double> (num, num);
179  }
180 
181  Cholesky(double **a, const int nin) : num(nin)
182  {
183  el = matrix <double> (num, num);
184  double sum = 0;
185  int i, j, k;
186 
187  for (i = 0; i < num; i++)
188  for (j = 0; j < num; j++)
189  el[i][j] = a[i][j];
190 
191  for (i = 0; i < num; i++)
192  {
193  for (j = i; j < num; j++)
194  {
195  for(sum = el[i][j], k = i - 1; k >= 0; k--)
196  sum -= el[i][k]*el[j][k];
197  if(i ==j)
198  {
199  if(sum <= 0.0)
200  {
201  std::cout << "Cholesky failed " << sum << std::endl;
202  getchar();
203  }
204  el[i][i] = sqrt(sum);
205  }
206  else
207  el[j][i] = sum/el[i][i];
208  }
209  }
210  for (i = 0; i < num; i++)
211  for (j = 0; j < i; j++)
212  el[j][i] = 0.0;
213  }
214 
215  bool EnterMatM(double **a, const int min)
216  {
217  double sum = 0;
218  int i, j, k;
219  //ofstream out1("cov.dat");
220  for (i = 0; i < num; i++)
221  {
222  for (j = 0; j < num; j++)
223  {
224  el[i][j] = a[i][j];
225  //out1 << el[i][j] << " " << endl;
226  }
227  //out1 << endl;
228  }
229 
230  if (num >= min)
231  {
232  k = min -1;
233  int l = num - min + 1;
234  for (i = 0; i < l; i++)
235  {
236  for (j = k+i; j < num; j++)
237  el[i][j] = el[j][i] = 0.0;
238  }
239  }
240 
241  for (i = 0; i < num; i++)
242  {
243  for (j = i; j < num; j++)
244  {
245  for(sum = el[i][j], k = i - 1; k >= 0; k--)
246  sum -= el[i][k]*el[j][k];
247  if(i ==j)
248  {
249  if(sum <= 0.0)
250  return true;
251  el[i][i] = sqrt(sum);
252  }
253  else
254  el[j][i] = sum/el[i][i];
255  }
256  }
257  for (i = 0; i < num; i++)
258  for (j = 0; j < i; j++)
259  el[j][i] = 0.0;
260 
261 // ofstream out2("cov2.dat");
262 // double **ct = matrix <double> (n, n, 0.0);
263 // for (i = 0; i < num; i++)
264 // for (j = 0; j < n; j++)
265 // for (k = 0; k < num; k++)
266 // ct[i][j] += el[i][k]*el[j][k];
267 // for (i = 0; i < num; i++)
268 // {
269 // for (j = 0; j < num; j++)
270 // out2 << ct[i][j] << " ";
271 // out2 << endl;
272 // }
273 // cout << "finished" << endl;
274 // getchar();
275 
276  return false;
277  }
278 
279  bool EnterMat(double **a)
280  {
281  double sum = 0;
282  int i, j, k;
283 
284  for (i = 0; i < num; i++)
285  {
286  for (j = 0; j < num; j++)
287  {
288  el[i][j] = a[i][j];
289  }
290  }
291 
292  for (i = 0; i < num; i++)
293  {
294  for (j = i; j < num; j++)
295  {
296  for(sum = el[i][j], k = i - 1; k >= 0; k--)
297  sum -= el[i][k]*el[j][k];
298  if(i ==j)
299  {
300  if(sum <= 0.0)
301  return false;
302  el[i][i] = sqrt(sum);
303  }
304  else
305  el[j][i] = sum/el[i][i];
306  }
307  }
308  for (i = 0; i < num; i++)
309  for (j = 0; j < i; j++)
310  el[j][i] = 0.0;
311 
312  return true;
313  }
314 
315  void EnterMat(double **a, int nin)
316  {
317  del <double> (el, num);
318  num = nin;
319  el = matrix <double> (num, num);
320  EnterMat(a);
321  }
322 
323  void ElMult (double *y, double *b)
324  {
325  int i, j;
326  for(i = 0; i < num; i++)
327  {
328  b[i] = 0.0;
329  for (j = 0; j <= i; j++)
330  b[i] += el[i][j]*y[j];
331  }
332  }
333 
334  void ElMult (double *y)
335  {
336  int i, j;
337  double *b = new double[num];
338  for(i = 0; i < num; i++)
339  {
340  b[i] = 0.0;
341  for (j = 0; j <= i; j++)
342  {
343  b[i] += el[i][j]*y[j];
344  }
345  }
346  for (i = 0; i < num; i++)
347  {
348  y[i] = b[i];
349  }
350  delete[] b;
351  }
352 
353  void Solve (double *b, double *x)
354  {
355  int i, k;
356  double sum;
357 
358  for (i = 0; i < num; i++)
359  {
360  for (sum = b[i], k=i-1; k >=0; k--)
361  sum -= el[i][k]*x[k];
362  x[i]=sum/el[i][i];
363  }
364 
365  for (i = num-1; i >=0; i--)
366  {
367  for (sum = x[i], k=i+1; k<num; k++)
368  sum -= el[k][i]*x[k];
369  x[i]=sum/el[i][i];
370  }
371  }
372 
373  double Square(double *y, double *y0)
374  {
375  int i, j;
376  double sum;
377  double *x = new double[num];
378 
379  for (i = 0; i < num; i++)
380  {
381  for (sum = (y[i]-y0[i]), j=0; j < i; j++)
382  sum -= el[i][j]*x[j];
383  x[i]=sum/el[i][i];
384  }
385 
386  sum = 0.0;
387  for (i = 0; i < num; i++)
388  sum += x[i]*x[i];
389 
390  delete[] x;
391 
392  return sum;
393  }
394 
395  double Square(double *y, double *y0, int *map)
396  {
397  int i, j;
398  double sum;
399  double *x = new double[num];
400 
401  for (i = 0; i < num; i++)
402  {
403  for (sum = (y[map[i]]-y0[i]), j=0; j < i; j++)
404  sum -= el[i][j]*x[j];
405  x[i]=sum/el[i][i];
406  }
407 
408  sum = 0.0;
409  for (i = 0; i < num; i++)
410  sum += x[i]*x[i];
411 
412  delete[] x;
413 
414  return sum;
415  }
416 
417  void Inverse(double **ainv)
418  {
419  double sum;
420 
421  for (int i = 0; i < num; i++)
422  for (int j = 0; j <= i; j++)
423  {
424  sum = (i == j ? 1.0 : 0.0);
425  for (int k = i-1; k >= j; k--)
426  sum -= el[i][k]*ainv[j][k];
427  ainv[j][i] = sum/el[i][i];
428  }
429 
430  for (int i = num - 1; i >= 0; i--)
431  for (int j = 0; j <= i; j++)
432  {
433  sum = (i < j ? 0.0 : ainv[j][i]);
434  for (int k = i + 1; k < num; k++)
435  sum -= el[k][i]*ainv[j][k];
436  ainv[i][j] = ainv[j][i] = sum/el[i][i];
437  }
438  }
439 
440  double DetSqrt()
441  {
442  double temp = 1.0;
443  for (int i = 0; i < num; i++)
444  temp *= el[i][i];
445  return temp;
446  }
448  {
449  del <double> (el, num);
450  }
451 };
452 
453 class Ran_old
454 {
455  private:
456  unsigned long long int u, v, w;
457 
458  public:
459  Ran_old(unsigned long long int j) : v(4101842887655102017LL), w(1)
460  {
461  u = j ^ v; int64();
462  v = u; int64(); int64();
463  w = v; int64(); int64();
464  }
465  inline unsigned long long int int64()
466  {
467  u = u * 2862933555777941757LL + 7046029254386353087LL;
468  v ^= v >> 17; v ^= v << 31; v ^= v >> 8;
469  w = 4294957665U*(w & 0xffffffff) + (w >> 32);
470  unsigned long long int x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;
471  return (x + v) ^ w;
472  }
473  inline double Doub(){return 5.42101086242752217E-20 * int64();}
474  inline unsigned int int32(){return (unsigned int)int64();}
475 };
476 
477 class Ran
478 {
479 public:
480  Ran(unsigned long long int){}
481  inline double Doub(){return Gambit::Random::draw();}
482 };
483 
484 class ExponDev : public Ran
485 {
486  private:
487  double beta;
488  public:
489  ExponDev(double din, unsigned long long ii) : Ran(ii), beta(din){}
490  double Dev()
491  {
492  double u;
493  do
494  {
495  u = Doub();
496  }
497  while(u == 0);
498 
499  return -log(u)/beta;
500  }
501 };
502 
503 class NormalDev : public Ran
504 {
505  private:
506  double mu, sig;
507 
508  public:
509  NormalDev(double mmu, double ssig, unsigned long long i) : Ran(i), mu(mmu), sig(ssig){}
510  double Dev()
511  {
512  double u, v, x, y, q;
513  do
514  {
515  u = Doub();
516  v = 1.7156*(Doub() - 0.5);
517  x = u - 0.449871;
518  y = fabs(v) + 0.386595;
519  q = x*x + y*(0.19600*y-0.25472*x);
520  }
521  while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u));
522 
523  return mu + sig*v/u;
524  }
525 };
526 
527 class BasicDevs : public Ran
528 {
529  private:
530 
531  public:
532  BasicDevs(unsigned long long i) : Ran(i) {}
533 
534  double Dev()
535  {
536  double u, v, x, y, q;
537  do
538  {
539  u = Doub();
540  v = 1.7156*(Doub() - 0.5);
541  x = u - 0.449871;
542  y = fabs(v) + 0.386595;
543  q = x*x + y*(0.19600*y-0.25472*x);
544  }
545  while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u));
546 
547  return v/u;
548  }
549 
550  double ExpDev()
551  {
552  double u;
553  do
554  {
555  u = Doub();
556  }
557  while(u == 0);
558 
559  return -log(u);
560  }
561 };
562 
563 class MultiNormalDev : public Ran, public Cholesky
564 {
565  private:
566  int mm;
567  double f;
568  double *spt;
569 
570  public:
571  MultiNormalDev (double **vvar, double fin, unsigned long long int j, int nin) : Ran(j), Cholesky(vvar, nin), mm(nin), f(fin)
572  {
573  spt = matrix <double> (mm);
574  }
575  void Dev(double *pt, double *mean)
576  {
577  double u, v, x, y, q;
578  for (int i = 0; i < mm; i++)
579  {
580  do
581  {
582  u = Doub();
583  v = 1.7156*(Doub() - 0.5);
584  x = u - 0.449871;
585  y = fabs(v) + 0.386595;
586  q = x*x + y*(0.19600*y-0.25472*x);
587  }
588  while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u));
589  spt[i] = v/u;
590  }
591  ElMult(spt, pt);
592  for (int i = 0; i < mm; i++)
593  pt[i] = f*pt[i] + mean[i];
594  }
595  void Dev(double **cvar, double *pt, double *mean)
596  {
597  double u, v, x, y, q;
598  for (int i = 0; i < mm; i++)
599  {
600  do
601  {
602  u = Doub();
603  v = 1.7156*(Doub() - 0.5);
604  x = u - 0.449871;
605  y = fabs(v) + 0.386595;
606  q = x*x + y*(0.19600*y-0.25472*x);
607  }
608  while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u));
609  spt[i] = v/u;
610  }
611  EnterMat(cvar);
612  ElMult(spt, pt);
613  for (int i = 0; i < mm; i++)
614  pt[i] = f*pt[i] + mean[i];
615  }
617  {
618  del <double> (spt);
619  }
620 };
621 
622 class AdvanceDevs : public BasicDevs, public Cholesky
623 {
624  private:
625  double fac;
626 
627  public:
628  AdvanceDevs(int nin, double din, unsigned long long iin) : BasicDevs(iin), Cholesky(nin), fac(din)
629  {
630  }
631 
632  AdvanceDevs(double **vvar, const int nin, double din, unsigned long long iin) : BasicDevs(iin), Cholesky(vvar, nin), fac(din)
633  {
634  }
635 
636  void MultiDev(double *pin, double *p0)
637  {
638  int i;
639  double dist = 0.0;
640 
641  double vec[num];
642  double norm = 0.0;
643  for (i = 0; i < num; i++)
644  {
645  vec[i] = Dev();
646  norm += vec[i]*vec[i];
647  }
648  norm = sqrt(norm);
649 
650  if(Doub() < 0.33)
651  dist = ExpDev();
652  else
653  {
654  double temp;
655  for (i = 0; i < 2; i++)
656  {
657  temp = Dev();
658  dist += temp*temp;
659  }
660  dist = sqrt(dist/2.0);
661  }
662 
663  ElMult(vec, pin);
664  for (i = 0; i < num; i++)
665  {
666  pin[i] = fac*dist*pin[i]/norm + p0[i];
667  }
668 
669  return;
670  }
671 
672  void MultiDev(double **cvar, double *pin, double *p0)
673  {
674  int i;
675  double dist = 0.0;
676 
677  double vec[num];
678  double norm = 0.0;
679  for (i = 0; i < num; i++)
680  {
681  vec[i] = Dev();
682  norm += vec[i]*vec[i];
683  }
684  norm = sqrt(norm);
685 
686  if(Doub() < 0.33)
687  dist = ExpDev();
688  else
689  {
690  double temp;
691  for (i = 0; i < 2; i++)
692  {
693  temp = Dev();
694  dist += temp*temp;
695  }
696  dist = sqrt(dist/2.0);
697  }
698 
699  EnterMat(cvar);
700  ElMult(vec, pin);
701  for (i = 0; i < num; i++)
702  {
703  pin[i] = fac*dist*pin[i]/norm + p0[i];
704  }
705 
706  return;
707  }
708 
709  void EllipseDev(double *pin, double *p0, double fin)
710  {
711  int i;
712  double dist = 0.0;
713 
714  double vec[num];
715  double norm = 0.0;
716  for (i = 0; i < num; i++)
717  {
718  vec[i] = Dev();
719  norm += vec[i]*vec[i];
720  }
721  norm = sqrt(norm);
722 
723  dist = pow(Doub(), 1.0/num);
724 
725  ElMult(vec, pin);
726  for (i = 0; i < num; i++)
727  {
728  pin[i] = fin*dist*pin[i] + p0[i];
729  }
730 
731  return;
732  }
733 
734  void EllipseDev(double **cvar, double *pin, double *p0, double fin)
735  {
736  int i;
737  double dist = 0.0;
738 
739  double vec[num];
740  double norm = 0.0;
741  for (i = 0; i < num; i++)
742  {
743  vec[i] = Dev();
744  norm += vec[i]*vec[i];
745  }
746  norm = sqrt(norm);
747 
748  dist = pow(Doub(), 1.0/num);
749 
750  EnterMat(cvar);
751  ElMult(vec, pin);
752  for (i = 0; i < num; i++)
753  {
754  pin[i] = fin*dist*pin[i] + p0[i];
755  }
756 
757  return;
758  }
759 };
760 
761 class RandomPlane : public AdvanceDevs
762 {
763 private:
764  double **rotVec;
765  double **currentVec;
766  double **endVec;
767  double **endEndVec;
768  double **sEndVec;
769  int proj, extra;
770  double alim, alimt;
771  double sqrtAlim, powAlim, ratioAlim, ratioAlimt;
772 
773 public:
774  RandomPlane(const int projin, const int nin, const double din, const double alim, const double alimt, unsigned long long iin) : AdvanceDevs(nin, din, iin), proj(projin), alim(alim), alimt(alimt)
775  {
776  rotVec = matrix <double> (nin, nin);
777  RandRot();
778  if (proj > num) proj = num;
779  extra = num % proj;
780  endVec = currentVec + num - extra;
781  endEndVec = currentVec + num;
782  sEndVec = currentVec + num - proj;
783  sqrtAlim = sqrt(alim);
784  powAlim = pow(alim, 0.5 - proj);
785  double vol1 = 2.0*(sqrt(alim) - 1.0);
786  double vol2 = (1.0 - pow(alim, 0.5 - proj))/(proj - 0.5);
787  ratioAlim = vol1/(vol1 + vol2);
788  ratioAlimt = (alimt - 1.0)/(2.0*alimt + proj - 2.0);
789  }
790 
791  double WalkDev()
792  {
793  if (Doub() <= ratioAlim)
794  return SQ(1.0 + Doub()*(sqrtAlim - 1.0));
795  else
796  return pow(powAlim + Doub()*(1.0 - powAlim), 1.0/(proj - 0.5));
797  }
798 
799  double TransDev()
800  {
801  if (Doub() < ratioAlimt)
802  {
803  return -pow(Doub(), 1.0/(alimt + proj - 1.0));
804  }
805  else
806  {
807  return -pow(Doub(), 1.0/(1.0 - alimt));
808  }
809  }
810 
811  double KWalkDev()
812  {
813  if (Doub() < 0.5)
814  {
815  if (Doub() <= ratioAlim)
816  {
817  return SQ(1.0 + Doub()*(sqrtAlim - 1.0));
818  }
819  else
820  {
821  return pow(powAlim + Doub()*(1.0 - powAlim), 1.0/(proj - 0.5));
822  }
823  }
824  else
825  {
826  if (Doub() < ratioAlimt)
827  {
828  return -pow(Doub(), 1.0/(alimt + proj - 1.0));
829  }
830  else
831  {
832  return -pow(Doub(), 1.0/(1.0 - alimt));
833  }
834  }
835  }
836 
837  bool KWalkDev(double &Z)
838  {
839  if (Doub() < 0.5)
840  {
841  Z = SQ(1.0/sqrtAlim + Doub()*(sqrtAlim - 1.0/sqrtAlim));
842 
843  return pow(Z, proj - 1) > Doub();
844  }
845  else
846  {
847  if (Doub() < (alimt - 1.0)/(2.0*alimt))
848  {
849  Z = -pow(Doub(), 1.0/(alimt + 1.0));
850  }
851  else
852  {
853  Z = -pow(Doub(), 1.0/(1.0 - alimt));
854  }
855 
856  return pow(Z, proj - 2) > Doub();
857  }
858  }
859 
860  double WalkDev(double *ptrOut, double *ptr, double *ptr0)
861  {
862  double Z;
863 
864  Z = SQ(1.0/sqrtAlim + Doub()*(sqrtAlim - 1.0/sqrtAlim));
865 
866  Mult(ptrOut, ptr, ptr0, Z);
867 
868  return (proj - 1.0)*log(Z);
869  }
870 
871  double TransDev(double *ptrOut, double *ptr, double *ptr0)
872  {
873  double Z;
874 
875  if (Doub() < (alimt - 1.0)/(2.0*alimt))
876  {
877  Z = -pow(Doub(), 1.0/(alimt + 1.0));
878  }
879  else
880  {
881  Z = -pow(Doub(), 1.0/(1.0 - alimt));
882  }
883 
884  Mult(ptrOut, ptr, ptr0, Z);
885 
886  return (proj - 1.0)*log(-Z);
887  }
888 
889  void Mult(double *ptrOut, double *ptr, double *ptr0, const double Z)
890  {
891  double z[proj];
892  for (int i = 0; i < proj; i++)
893  {
894  z[i] = 0.0;
895  for (int j = 0; j < num; j++)
896  {
897  z[i] += (Z-1.0)*(ptr[j] - ptr0[j])*currentVec[i][j];
898  }
899  }
900 
901  for (int j = 0; j < num; j++)
902  {
903  ptrOut[j] = ptr[j];
904  for (int i = 0; i < proj; i++)
905  {
906  ptrOut[j] += z[i]*currentVec[i][j];
907  }
908  }
909 
910  currentVec++;// += proj;
911  if (currentVec == sEndVec)
912  {
913  for (double **vec = currentVec, **vec2 = rotVec; vec != endEndVec; vec++, vec2++)
914  {
915  double *temp = *vec2;
916  *vec2 = *vec;
917  *vec = temp;
918  }
919  RandRot(extra);
920  }
921  }
922 
923  void Mult2(double *ptrOut, double *ptr, double *ptr0, const double Z)
924  {
925  RandRot(0, proj);
926 
927  double z[proj];
928  for (int i = 0; i < proj; i++)
929  {
930  z[i] = 0.0;
931  for (int j = 0; j < num; j++)
932  {
933  z[i] += (Z-1.0)*(ptr[j] - ptr0[j])*currentVec[i][j];
934  }
935  }
936 
937  for (int j = 0; j < num; j++)
938  {
939  ptrOut[j] = ptr[j];
940  for (int i = 0; i < proj; i++)
941  {
942  ptrOut[j] += z[i]*currentVec[i][j];
943  }
944  }
945  }
946 
947  void RandRot(const int start = 0)
948  {
949  double temp;
950  double vec[num];
951  int i, j, k;
952 
953  for (i = start; i < num; i++)
954  {
955  for (j = 0; j < num; j++)
956  vec[j] = Dev();
957  for (j = 0; j < i; j++)
958  {
959  temp = 0.0;
960  for (k = 0; k < num; k++)
961  temp += vec[k]*rotVec[j][k];
962  for (k = 0; k < num; k++)
963  vec[k] -= temp*rotVec[j][k];
964  }
965  temp = 0.0;
966  for (j = 0; j < num; j++)
967  temp += vec[j]*vec[j];
968  temp = sqrt(temp);
969  for (j = 0; j < num; j++)
970  rotVec[i][j] = vec[j]/temp;
971  }
972 
973  currentVec = rotVec;
974  }
975 
976  void RandRot(const int start, const int end)
977  {
978  double temp;
979  double vec[num];
980  int i, j, k;
981 
982  for (i = start; i < end; i++)
983  {
984  for (j = 0; j < num; j++)
985  vec[j] = Dev();
986  for (j = 0; j < i; j++)
987  {
988  temp = 0.0;
989  for (k = 0; k < num; k++)
990  temp += vec[k]*rotVec[j][k];
991  for (k = 0; k < num; k++)
992  vec[k] -= temp*rotVec[j][k];
993  }
994  temp = 0.0;
995  for (j = 0; j < num; j++)
996  temp += vec[j]*vec[j];
997  temp = sqrt(temp);
998  for (j = 0; j < num; j++)
999  rotVec[i][j] = vec[j]/temp;
1000  }
1001 
1002  currentVec = rotVec;
1003  }
1004 
1005  int Dim() const {return proj;}
1006 
1008  {
1009  del <double> (rotVec, num);
1010  }
1011 };
1012 
1013 class RandomBasis : public BasicDevs
1014 {
1015  private:
1016  double **rotVec;
1017 
1018  protected:
1019  int num;
1020  double **currentVec;
1021  double **endVec;
1022 
1023  public:
1024  RandomBasis(int nin, unsigned long long iin) : BasicDevs(iin), num(nin)
1025  {
1026  rotVec = matrix <double> (nin, nin);
1027  RandRot();
1028  endVec = currentVec + num;
1029  }
1030 
1031  void ChangeDim(const int nin)
1032  {
1033  del <double> (rotVec, num);
1034  num = nin;
1035  rotVec = matrix <double> (nin, nin);
1036  RandRot();
1037  endVec = currentVec + num;
1038  }
1039 
1040  void RandRot()
1041  {
1042  double temp;
1043  double vec[num];
1044  int i, j, k;
1045 
1046  for (i = 0; i < num; i++)
1047  {
1048  for (j = 0; j < num; j++)
1049  vec[j] = Dev();
1050  for (j = 0; j < i; j++)
1051  {
1052  temp = 0.0;
1053  for (k = 0; k < num; k++)
1054  temp += vec[k]*rotVec[j][k];
1055  for (k = 0; k < num; k++)
1056  vec[k] -= temp*rotVec[j][k];
1057  }
1058  temp = 0.0;
1059  for (j = 0; j < num; j++)
1060  temp += vec[j]*vec[j];
1061  temp = sqrt(temp);
1062  for (j = 0; j < num; j++)
1063  rotVec[i][j] = vec[j]/temp;
1064  }
1065 
1066  currentVec = rotVec;
1067  }
1068 
1069  double RanMult(double **cin)
1070  {
1071  int i, j;
1072  double temp = 0.0;
1073 
1074  for (i = 0; i < num; i++)
1075  {
1076  for (j = 0; j < num; j++)
1077  {
1078  temp += (*currentVec)[i]*cin[i][j]*(*currentVec)[j];
1079  }
1080  }
1081  return temp;
1082  }
1083 
1084  void RanMult(const double in, double *out)
1085  {
1086  int i;
1087  double *ptr = *currentVec;
1088 
1089  for (i = 0; i < num; i++)
1090  {
1091  *out++ = in*(*ptr++);
1092 
1093  }
1094  return;
1095  }
1096 
1097  void RanMult(double *in, const double w, double *out)
1098  {
1099  int i;
1100  double *ptr = *currentVec;
1101 
1102  for (i = 0; i < num; i++)
1103  {
1104  *out++ = *in++ + w*(*ptr++);
1105  }
1106  return;
1107  }
1108 
1109  double Mag(double *a, double *a0)
1110  {
1111  double temp = 0.0;
1112  int i;
1113  double *ptr = *currentVec;
1114 
1115  for (i = 0; i < num; i++)
1116  temp += (*a++-*a0++)*(*ptr++);
1117 
1118  return temp;
1119  }
1120 
1121  void Adjust(double *a, const double lim, const int iin)
1122  {
1123  double temp = (lim - a[iin])/(*currentVec)[iin];
1124  double *ptr = *currentVec;
1125 
1126  for (int i = 0; i < num; i++)
1127  a[i] += temp*(*ptr++);
1128  return;
1129  }
1130 
1131  virtual void operator ++ (int)
1132  {
1133  if (++currentVec == endVec)
1134  {
1135  RandRot();
1136  }
1137  }
1138 
1139  virtual ~RandomBasis()
1140  {
1141  del <double> (rotVec, num);
1142  }
1143 };
1144 
1146 {
1147  private:
1148  int num;
1149 
1150  public:
1151  TransformRandomBasis(double **vvar, int nin, unsigned long long iin) : RandomBasis(nin, iin), Cholesky(vvar, nin), num(nin)
1152  {
1153  double **ptr = currentVec;
1154  for (int i = 0; i < num; i++)
1155  {
1156  ElMult(*ptr++);
1157  }
1158  }
1159 
1160  void operator ++ (int)
1161  {
1162  if (++currentVec == endVec)
1163  {
1164  RandRot();
1165  double **ptr = currentVec;
1166  for (int i = 0; i < num; i++)
1167  {
1168  ElMult(*ptr++);
1169  }
1170  }
1171  }
1172 };
1173 
1174 class MultiNormDev : public RandomBasis, public Cholesky
1175 {
1176  private:
1177  double fac;
1178  int num;
1179 
1180  public:
1181  MultiNormDev(int nin, double din, unsigned long long iin) : RandomBasis(nin, iin), Cholesky(nin), fac(din), num(nin)
1182  {
1183  }
1184 
1185  MultiNormDev(double **vvar, const int nin, double din, unsigned long long iin) : RandomBasis(nin, iin), Cholesky(vvar, nin), fac(din), num(nin)
1186  {
1187  }
1188 
1189  void MultiDev(double *pin, double *p0)
1190  {
1191  int i;
1192  double dist = 0.0;
1193 
1194  if(Doub() < 0.33)
1195  dist = ExpDev();
1196  else
1197  {
1198  double temp;
1199  for (i = 0; i < 2; i++)
1200  {
1201  temp = Dev();
1202  dist += temp*temp;
1203  }
1204  dist = sqrt(dist/2.0);
1205  }
1206 
1207  ElMult(*currentVec, pin);
1208  for (i = 0; i < num; i++)
1209  {
1210  pin[i] = fac*dist*pin[i] + p0[i];
1211  }
1212  (*this)++;
1213 
1214  return;
1215  }
1216 
1217  void MultiDev(double **cvar, double *pin, double *p0)
1218  {
1219  int i;
1220  double dist = 0.0;
1221 
1222  if(Doub() < 0.33)
1223  dist = ExpDev();
1224  else
1225  {
1226  double temp;
1227  for (i = 0; i < 2; i++)
1228  {
1229  temp = Dev();
1230  dist += temp*temp;
1231  }
1232  dist = sqrt(dist/2.0);
1233  }
1234 
1235  EnterMat(cvar);
1236  ElMult(*currentVec, pin);
1237  for (i = 0; i < num; i++)
1238  {
1239  pin[i] = fac*dist*pin[i] + p0[i];
1240  }
1241  (*this)++;
1242 
1243  return;
1244  }
1245 
1246  void MultiDevGauss(double **cvar, double *pin, double *p0)
1247  {
1248  int i;
1249  double dist = 0.0;
1250 
1251  double temp;
1252  for (i = 0; i < 2; i++)
1253  {
1254  temp = Dev();
1255  dist += temp*temp;
1256  }
1257  dist = sqrt(dist/2.0);
1258 
1259 
1260  EnterMat(cvar);
1261  ElMult(*currentVec, pin);
1262  for (i = 0; i < num; i++)
1263  {
1264  pin[i] = fac*dist*pin[i] + p0[i];
1265  }
1266  (*this)++;
1267 
1268  return;
1269  }
1270 
1271  void EllipseDev(double *pin, double *p0, double fin)
1272  {
1273  int i;
1274  double dist = 0.0;
1275 
1276  dist = pow(Doub(), 1.0/num);
1277 
1278  ElMult(*currentVec, pin);
1279  for (i = 0; i < num; i++)
1280  {
1281  pin[i] = fin*dist*pin[i] + p0[i];
1282  }
1283  (*this)++;
1284 
1285  return;
1286  }
1287 
1288  void EllipseDev(double **cvar, double *pin, double *p0, double fin)
1289  {
1290  int i;
1291  double dist = 0.0;
1292 
1293  dist = pow(Doub(), 1.0/num);
1294 
1295  EnterMat(cvar);
1296  ElMult(*currentVec, pin);
1297  for (i = 0; i < num; i++)
1298  {
1299  pin[i] = fin*dist*pin[i] + p0[i];
1300  }
1301  (*this)++;
1302 
1303  return;
1304  }
1305 };
1306 
1307 #endif
void del(T *temp)
unsigned long long int int64()
Cholesky(const int nin)
Ran(unsigned long long int)
unsigned int int32()
double ** rotVec
Ran_old(unsigned long long int j)
void EllipseDev(double *pin, double *p0, double fin)
void Adjust(double *a, const double lim, const int iin)
START_MODEL beta
Definition: Axions.hpp:35
double Dev()
unsigned long long int w
double TransDev(double *ptrOut, double *ptr, double *ptr0)
RandomBasis(int nin, unsigned long long iin)
NormalDev(double mmu, double ssig, unsigned long long i)
void Dev(double **cvar, double *pt, double *mean)
const T SQ(const T a)
double RanMult(double **cin)
void EllipseDev(double **cvar, double *pin, double *p0, double fin)
bool EnterMat(double **a)
double ** currentVec
MultiNormDev(int nin, double din, unsigned long long iin)
double ** endVec
void Mult2(double *ptrOut, double *ptr, double *ptr0, const double Z)
double TransDev()
double ** sEndVec
void ElMult(double *y)
T * matrix(const int xN)
Definition: random_tools.hpp:9
BasicDevs(unsigned long long i)
AdvanceDevs(int nin, double din, unsigned long long iin)
RandomPlane(const int projin, const int nin, const double din, const double alim, const double alimt, unsigned long long iin)
double WalkDev(double *ptrOut, double *ptr, double *ptr0)
MultiNormalDev(double **vvar, double fin, unsigned long long int j, int nin)
void Inverse(double **ainv)
START_MODEL b
Definition: demo.hpp:235
void Dev(double *pt, double *mean)
double ** currentVec
void RandRot(const int start, const int end)
void MultiDev(double *pin, double *p0)
double ** rotVec
int Dim() const
ExponDev(double din, unsigned long long ii)
double Dev()
double ** endVec
double WalkDev()
Cholesky(double **a, const int nin)
double beta
static double draw()
Draw a single uniform random deviate from the interval (0,1) using the chosen RNG engine...
bool KWalkDev(double &Z)
void ChangeDim(const int nin)
double Square(double *y, double *y0, int *map)
const double mu
Definition: SM_Z.hpp:42
void MultiDevGauss(double **cvar, double *pin, double *p0)
dictionary fin
double Doub()
double ** endEndVec
TransformRandomBasis(double **vvar, int nin, unsigned long long iin)
void MultiDev(double **cvar, double *pin, double *p0)
void MultiDev(double *pin, double *p0)
double KWalkDev()
A threadsafe interface to the STL random number generators.
void Solve(double *b, double *x)
void Mult(double *ptrOut, double *ptr, double *ptr0, const double Z)
void EllipseDev(double *pin, double *p0, double fin)
MultiNormDev(double **vvar, const int nin, double din, unsigned long long iin)
double Doub()
double ** el
AdvanceDevs(double **vvar, const int nin, double din, unsigned long long iin)
void ElMult(double *y, double *b)
void RandRot(const int start=0)
void EllipseDev(double **cvar, double *pin, double *p0, double fin)
virtual ~RandomBasis()
double pow(const double &a)
Outputs a^i.
double Mag(double *a, double *a0)
double Square(double *y, double *y0)
void EnterMat(double **a, int nin)
void RanMult(double *in, const double w, double *out)
std::vector< T > vec(std::vector< T > vector)
Definition: daFunk.hpp:142
double Dev()
double ExpDev()
bool EnterMatM(double **a, const int min)
double DetSqrt()
const T SQR(const T a)
void MultiDev(double **cvar, double *pin, double *p0)
void RanMult(const double in, double *out)