|
GAMBIT
v1.5.0-2191-ga4742ac
a Global And Modular Bsm Inference Tool
|
Go to the documentation of this file. 1 #ifndef RANDOM_TOOLS_HPP 2 #define RANDOM_TOOLS_HPP 17 T ** matrix( const int xN, const int yN) 19 T **temp = new T*[xN]; 21 for ( int i = 0; i < xN; i++) 28 T *** matrix( const int xN, const int yN, const int zN) 30 T ***temp = new T**[xN]; 32 for ( int i = 0; i < xN; i++) 35 for ( int j = 0; j < yN; j++) 36 temp[i][j] = new T[zN]; 46 for ( int i = 0; i < xN; i++) 53 T ** matrix( const int xN, const int yN, T in) 55 T **temp = new T*[xN]; 57 for ( int i = 0; i < xN; i++) 60 for ( int j = 0; j < yN; j++) 68 T *** matrix( const int xN, const int yN, const int zN, T in) 70 T ***temp = new T**[xN]; 72 for ( int i = 0; i < xN; i++) 75 for ( int j = 0; j < yN; j++) 77 temp[i][j] = new T[zN]; 78 for ( int k = 0; k < zN; k++) 89 for ( int i = 0; i < xN; i++) 96 T ** matrix( const int xN, const int yN, T **in) 98 T **temp = new T*[xN]; 100 for ( int i = 0; i < xN; i++) 103 for ( int j = 0; j < yN; j++) 104 temp[i][j] = in[i][j]; 110 template < typename T> 111 T *** matrix( const int xN, const int yN, const int zN, T ***in) 113 T ***temp = new T**[xN]; 115 for ( int i = 0; i < xN; i++) 117 temp[i] = new T*[yN]; 118 for ( int j = 0; j < yN; j++) 120 temp[i][j] = new T[zN]; 121 for ( int k = 0; k < zN; k++) 122 temp[i][j][k] = in[i][j][k]; 129 template < typename T> 136 template < typename T> 137 void del(T **temp, int xN) 139 for ( int i = 0; i < xN; i++) 146 template < typename T> 147 void del(T ***temp, int xN, int yN) 149 for ( int i = 0; i < xN; i++) 151 for ( int j = 0; j < yN; j++) 161 inline const T SQ( const T a) { return a*a;} 164 inline const T SQR( const T a) { return a*a;} 178 el = matrix <double> ( num, num); 183 el = matrix <double> ( num, num); 187 for (i = 0; i < num; i++) 188 for (j = 0; j < num; j++) 191 for (i = 0; i < num; i++) 193 for (j = i; j < num; j++) 195 for(sum = el[i][j], k = i - 1; k >= 0; k--) 196 sum -= el[i][k]*el[j][k]; 201 std::cout << "Cholesky failed " << sum << std::endl; 204 el[i][i] = sqrt(sum); 207 el[j][i] = sum/el[i][i]; 210 for (i = 0; i < num; i++) 211 for (j = 0; j < i; j++) 220 for (i = 0; i < num; i++) 222 for (j = 0; j < num; j++) 233 int l = num - min + 1; 234 for (i = 0; i < l; i++) 236 for (j = k+i; j < num; j++) 237 el[i][j] = el[j][i] = 0.0; 241 for (i = 0; i < num; i++) 243 for (j = i; j < num; j++) 245 for(sum = el[i][j], k = i - 1; k >= 0; k--) 246 sum -= el[i][k]*el[j][k]; 251 el[i][i] = sqrt(sum); 254 el[j][i] = sum/el[i][i]; 257 for (i = 0; i < num; i++) 258 for (j = 0; j < i; j++) 284 for (i = 0; i < num; i++) 286 for (j = 0; j < num; j++) 292 for (i = 0; i < num; i++) 294 for (j = i; j < num; j++) 296 for(sum = el[i][j], k = i - 1; k >= 0; k--) 297 sum -= el[i][k]*el[j][k]; 302 el[i][i] = sqrt(sum); 305 el[j][i] = sum/el[i][i]; 308 for (i = 0; i < num; i++) 309 for (j = 0; j < i; j++) 317 del <double> ( el, num); 319 el = matrix <double> ( num, num); 326 for(i = 0; i < num; i++) 329 for (j = 0; j <= i; j++) 330 b[i] += el[i][j]*y[j]; 337 double * b = new double[ num]; 338 for(i = 0; i < num; i++) 341 for (j = 0; j <= i; j++) 343 b[i] += el[i][j]*y[j]; 346 for (i = 0; i < num; i++) 358 for (i = 0; i < num; i++) 360 for (sum = b[i], k=i-1; k >=0; k--) 361 sum -= el[i][k]*x[k]; 365 for (i = num-1; i >=0; i--) 367 for (sum = x[i], k=i+1; k< num; k++) 368 sum -= el[k][i]*x[k]; 377 double *x = new double[ num]; 379 for (i = 0; i < num; i++) 381 for (sum = (y[i]-y0[i]), j=0; j < i; j++) 382 sum -= el[i][j]*x[j]; 387 for (i = 0; i < num; i++) 395 double Square( double *y, double *y0, int *map) 399 double *x = new double[ num]; 401 for (i = 0; i < num; i++) 403 for (sum = (y[map[i]]-y0[i]), j=0; j < i; j++) 404 sum -= el[i][j]*x[j]; 409 for (i = 0; i < num; i++) 421 for ( int i = 0; i < num; i++) 422 for ( int j = 0; j <= i; j++) 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]; 430 for ( int i = num - 1; i >= 0; i--) 431 for ( int j = 0; j <= i; j++) 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]; 443 for ( int i = 0; i < num; i++) 449 del <double> ( el, num); 456 unsigned long long int u, v, w; 459 Ran_old( unsigned long long int j) : v(4101842887655102017LL), w(1) 462 v = u; int64(); int64(); 463 w = v; int64(); int64(); 465 inline unsigned long long int int64() 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; 473 inline double Doub(){ return 5.42101086242752217E-20 * int64();} 474 inline unsigned int int32(){ return ( unsigned int)int64();} 480 Ran( unsigned long long int){} 489 ExponDev( double din, unsigned long long ii) : Ran(ii), beta(din){} 509 NormalDev( double mmu, double ssig, unsigned long long i) : Ran(i), mu(mmu), sig(ssig){} 512 double u, v, x, y, q; 516 v = 1.7156*(Doub() - 0.5); 518 y = fabs(v) + 0.386595; 519 q = x*x + y*(0.19600*y-0.25472*x); 521 while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u)); 536 double u, v, x, y, q; 540 v = 1.7156*(Doub() - 0.5); 542 y = fabs(v) + 0.386595; 543 q = x*x + y*(0.19600*y-0.25472*x); 545 while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u)); 573 spt = matrix <double> (mm); 575 void Dev( double *pt, double *mean) 577 double u, v, x, y, q; 578 for ( int i = 0; i < mm; i++) 583 v = 1.7156*(Doub() - 0.5); 585 y = fabs(v) + 0.386595; 586 q = x*x + y*(0.19600*y-0.25472*x); 588 while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u)); 592 for ( int i = 0; i < mm; i++) 593 pt[i] = f*pt[i] + mean[i]; 595 void Dev( double **cvar, double *pt, double *mean) 597 double u, v, x, y, q; 598 for ( int i = 0; i < mm; i++) 603 v = 1.7156*(Doub() - 0.5); 605 y = fabs(v) + 0.386595; 606 q = x*x + y*(0.19600*y-0.25472*x); 608 while(q > 0.27597 && (q > 0.27846 || v*v > -4.0*log(u)*u*u)); 613 for ( int i = 0; i < mm; i++) 614 pt[i] = f*pt[i] + mean[i]; 643 for (i = 0; i < num; i++) 646 norm += vec[i]*vec[i]; 655 for (i = 0; i < 2; i++) 660 dist = sqrt(dist/2.0); 664 for (i = 0; i < num; i++) 666 pin[i] = fac*dist*pin[i]/norm + p0[i]; 672 void MultiDev( double **cvar, double *pin, double *p0) 679 for (i = 0; i < num; i++) 682 norm += vec[i]*vec[i]; 691 for (i = 0; i < 2; i++) 696 dist = sqrt(dist/2.0); 701 for (i = 0; i < num; i++) 703 pin[i] = fac*dist*pin[i]/norm + p0[i]; 716 for (i = 0; i < num; i++) 719 norm += vec[i]*vec[i]; 723 dist = pow(Doub(), 1.0/num); 726 for (i = 0; i < num; i++) 728 pin[i] = fin*dist*pin[i] + p0[i]; 741 for (i = 0; i < num; i++) 744 norm += vec[i]*vec[i]; 748 dist = pow(Doub(), 1.0/num); 752 for (i = 0; i < num; i++) 754 pin[i] = fin*dist*pin[i] + p0[i]; 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) 776 rotVec = matrix <double> (nin, nin); 778 if (proj > num) proj = num; 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); 793 if (Doub() <= ratioAlim) 794 return SQ(1.0 + Doub()*(sqrtAlim - 1.0)); 796 return pow(powAlim + Doub()*(1.0 - powAlim), 1.0/(proj - 0.5)); 801 if (Doub() < ratioAlimt) 803 return - pow(Doub(), 1.0/(alimt + proj - 1.0)); 807 return - pow(Doub(), 1.0/(1.0 - alimt)); 815 if (Doub() <= ratioAlim) 817 return SQ(1.0 + Doub()*(sqrtAlim - 1.0)); 821 return pow(powAlim + Doub()*(1.0 - powAlim), 1.0/(proj - 0.5)); 826 if (Doub() < ratioAlimt) 828 return - pow(Doub(), 1.0/(alimt + proj - 1.0)); 832 return - pow(Doub(), 1.0/(1.0 - alimt)); 841 Z = SQ(1.0/sqrtAlim + Doub()*(sqrtAlim - 1.0/sqrtAlim)); 843 return pow(Z, proj - 1) > Doub(); 847 if (Doub() < (alimt - 1.0)/(2.0*alimt)) 849 Z = - pow(Doub(), 1.0/(alimt + 1.0)); 853 Z = - pow(Doub(), 1.0/(1.0 - alimt)); 856 return pow(Z, proj - 2) > Doub(); 860 double WalkDev( double *ptrOut, double *ptr, double *ptr0) 864 Z = SQ(1.0/sqrtAlim + Doub()*(sqrtAlim - 1.0/sqrtAlim)); 866 Mult(ptrOut, ptr, ptr0, Z); 868 return (proj - 1.0)*log(Z); 871 double TransDev( double *ptrOut, double *ptr, double *ptr0) 875 if (Doub() < (alimt - 1.0)/(2.0*alimt)) 877 Z = - pow(Doub(), 1.0/(alimt + 1.0)); 881 Z = - pow(Doub(), 1.0/(1.0 - alimt)); 884 Mult(ptrOut, ptr, ptr0, Z); 886 return (proj - 1.0)*log(-Z); 889 void Mult( double *ptrOut, double *ptr, double *ptr0, const double Z) 892 for ( int i = 0; i < proj; i++) 895 for ( int j = 0; j < num; j++) 897 z[i] += (Z-1.0)*(ptr[j] - ptr0[j])*currentVec[i][j]; 901 for ( int j = 0; j < num; j++) 904 for ( int i = 0; i < proj; i++) 906 ptrOut[j] += z[i]*currentVec[i][j]; 911 if (currentVec == sEndVec) 913 for ( double ** vec = currentVec, **vec2 = rotVec; vec != endEndVec; vec++, vec2++) 915 double *temp = *vec2; 923 void Mult2( double *ptrOut, double *ptr, double *ptr0, const double Z) 928 for ( int i = 0; i < proj; i++) 931 for ( int j = 0; j < num; j++) 933 z[i] += (Z-1.0)*(ptr[j] - ptr0[j])*currentVec[i][j]; 937 for ( int j = 0; j < num; j++) 940 for ( int i = 0; i < proj; i++) 942 ptrOut[j] += z[i]*currentVec[i][j]; 953 for (i = start; i < num; i++) 955 for (j = 0; j < num; j++) 957 for (j = 0; j < i; j++) 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]; 966 for (j = 0; j < num; j++) 967 temp += vec[j]*vec[j]; 969 for (j = 0; j < num; j++) 970 rotVec[i][j] = vec[j]/temp; 982 for (i = start; i < end; i++) 984 for (j = 0; j < num; j++) 986 for (j = 0; j < i; j++) 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]; 995 for (j = 0; j < num; j++) 996 temp += vec[j]*vec[j]; 998 for (j = 0; j < num; j++) 999 rotVec[i][j] = vec[j]/temp; 1002 currentVec = rotVec; 1005 int Dim() const { return proj;} 1009 del <double> (rotVec, num); 1026 rotVec = matrix <double> (nin, nin); 1028 endVec = currentVec + num; 1033 del <double> (rotVec, num); 1035 rotVec = matrix <double> (nin, nin); 1037 endVec = currentVec + num; 1046 for (i = 0; i < num; i++) 1048 for (j = 0; j < num; j++) 1050 for (j = 0; j < i; j++) 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]; 1059 for (j = 0; j < num; j++) 1060 temp += vec[j]*vec[j]; 1062 for (j = 0; j < num; j++) 1063 rotVec[i][j] = vec[j]/temp; 1066 currentVec = rotVec; 1074 for (i = 0; i < num; i++) 1076 for (j = 0; j < num; j++) 1078 temp += (*currentVec)[i]*cin[i][j]*(*currentVec)[j]; 1087 double *ptr = *currentVec; 1089 for (i = 0; i < num; i++) 1091 *out++ = in*(*ptr++); 1097 void RanMult( double *in, const double w, double *out) 1100 double *ptr = *currentVec; 1102 for (i = 0; i < num; i++) 1104 *out++ = *in++ + w*(*ptr++); 1109 double Mag( double *a, double *a0) 1113 double *ptr = *currentVec; 1115 for (i = 0; i < num; i++) 1116 temp += (*a++-*a0++)*(*ptr++); 1121 void Adjust( double *a, const double lim, const int iin) 1123 double temp = (lim - a[iin])/(*currentVec)[iin]; 1124 double *ptr = *currentVec; 1126 for ( int i = 0; i < num; i++) 1127 a[i] += temp*(*ptr++); 1131 virtual void operator ++ ( int) 1133 if (++currentVec == endVec) 1141 del <double> (rotVec, num); 1153 double **ptr = currentVec; 1154 for ( int i = 0; i < num; i++) 1160 void operator ++ ( int) 1162 if (++currentVec == endVec) 1165 double **ptr = currentVec; 1166 for ( int i = 0; i < num; i++) 1199 for (i = 0; i < 2; i++) 1204 dist = sqrt(dist/2.0); 1207 ElMult(*currentVec, pin); 1208 for (i = 0; i < num; i++) 1210 pin[i] = fac*dist*pin[i] + p0[i]; 1227 for (i = 0; i < 2; i++) 1232 dist = sqrt(dist/2.0); 1236 ElMult(*currentVec, pin); 1237 for (i = 0; i < num; i++) 1239 pin[i] = fac*dist*pin[i] + p0[i]; 1252 for (i = 0; i < 2; i++) 1257 dist = sqrt(dist/2.0); 1261 ElMult(*currentVec, pin); 1262 for (i = 0; i < num; i++) 1264 pin[i] = fac*dist*pin[i] + p0[i]; 1276 dist = pow(Doub(), 1.0/num); 1278 ElMult(*currentVec, pin); 1279 for (i = 0; i < num; i++) 1281 pin[i] = fin*dist*pin[i] + p0[i]; 1293 dist = pow(Doub(), 1.0/num); 1296 ElMult(*currentVec, pin); 1297 for (i = 0; i < num; i++) 1299 pin[i] = fin*dist*pin[i] + p0[i];
unsigned long long int int64()
Ran(unsigned long long int)
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)
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)
double RanMult(double **cin)
void EllipseDev(double **cvar, double *pin, double *p0, double fin)
bool EnterMat(double **a)
MultiNormDev(int nin, double din, unsigned long long iin)
void Mult2(double *ptrOut, double *ptr, double *ptr0, const double Z)
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)
void Dev(double *pt, double *mean)
void RandRot(const int start, const int end)
void MultiDev(double *pin, double *p0)
ExponDev(double din, unsigned long long ii)
Cholesky(double **a, const int nin)
static double draw() Draw a single uniform random deviate from the interval (0,1) using the chosen RNG engine...
void ChangeDim(const int nin)
double Square(double *y, double *y0, int *map)
void MultiDevGauss(double **cvar, double *pin, double *p0)
void MultiDev(double **cvar, double *pin, double *p0)
void MultiDev(double *pin, double *p0)
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)
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)
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)
bool EnterMatM(double **a, const int min)
void MultiDev(double **cvar, double *pin, double *p0)
void RanMult(const double in, double *out)
|