29 #include <Base/Common/W32Compat.hh>
38 #include <gsl/gsl_version.h>
39 #include <gsl/gsl_rng.h>
40 #include <gsl/gsl_randist.h>
70 : Seed_(0), Idum_(0), LongIdum_(0)
76 #ifdef RANDOM_SEED_ZERO_FIXED
83 : Seed_(0), Idum_(0), LongIdum_(0)
88 #ifdef RANDOM_SEED_ZERO_FIXED
89 static bool first =
true;
94 #ifdef RANDOM_SEED_ZERO_FIXED
105 srand((
long int)
Seed_);
115 #ifdef RANDOM_SEED_ZERO_FIXED
127 Seed_=(
unsigned int)(uc);
141 int min = (min_arg<max_arg)?(min_arg):(max_arg);
142 int max = (min_arg<max_arg)?(max_arg):(min_arg);
149 double nr = (double)r / (
double)(RAND_MAX);
150 int val = min + (int)(nr * (
double)(max-min));
152 BIASASSERT(val>=min && val<=max);
159 unsigned min = (min_arg<max_arg)?(min_arg):(max_arg);
160 unsigned max = (min_arg<max_arg)?(max_arg):(min_arg);
166 unsigned val = (unsigned)(dval);
167 if (val<min || val>max) val = (val % (max+1-min)) + min;
168 BIASASSERT(val>=min && val<=max);
179 for (
int i=0; i<num; i++)
187 static bool initialized =
false;
188 static gsl_rng * r = NULL;
190 cout <<
"initializing gsl_rng ! \n";
193 const gsl_rng_type * T = gsl_rng_mt19937;
195 r = gsl_rng_alloc (T);
200 return gsl_ran_gaussian(r, sigma) + mean;
206 BIASERR(
"GetNormalDistributedGSL(): No GSL support in BIAS, returning "
207 <<
"own implementation");
214 const int num, std::vector<double>& res)
217 for (
int i=0; i<num; i++)
229 vector<Vector2<double> > vran;
232 if (res==0) ran = vran[1];
244 double eigen_val1, eigen_val2;
252 R[0][0] = eigen_vec1[0];
253 R[1][0] = eigen_vec1[1];
254 R[0][1] = eigen_vec2[0];
255 R[1][1] = eigen_vec2[1];
258 vector<double> ran1, ran2;
268 ran2.resize(num, 0.0);
271 BIASERR(
"not a covariance matrix ("<<eigen_res<<
") "<<cov<<endl);
283 vector<Vector2<double> >::iterator it;
284 for (
int i=(
int)ran.size()-1; i>=0; i--){
285 ran[i].Set(ran1[i], ran2[i]);
306 BIASERR(
"quasi random sequences can be only generated for spaces with "
307 "dimension less than 7. Requested dimension is "<<req_dim);
310 static int dimension = -1;
312 if (dimension != req_dim || req_dim<=0){
322 if (min.
Size()!=max.
Size() || req_dim!=(int)min.
Size()){
323 BIASERR(
"argument vector size mismatch: "<<min.
Size()<<
" "<<max.
Size()
327 for (
int d=0; d<req_dim; d++){
329 BIASERR(
"max["<<d<<
"] smaller than min["<<d<<
"] ("<<max[d]
330 <<
" < "<<min[d]<<
")");
341 for (
int i=0; i<dimension; i++)
342 res[i] = r[i] * scale[i] + min[i];
354 for (
int d=0; d<req_dim; d++){
356 BIASERR(
"max["<<d<<
"] smaller than min["<<d<<
"] ("<<max[d]
357 <<
" < "<<min[d]<<
")");
361 if (min.
Size()!=max.
Size() || req_dim!=(int)min.
Size()){
362 BIASERR(
"argument vector size mismatch: "<<min.
Size()<<
" "<<max.
Size()
372 const int num_runs = req_dim/6 + (req_dim%6==0)?(0):(1);
375 res.resize(num, zero);
376 for (
int k=0; k<num_runs; k++){
381 dimension = (k==(num_runs-1))?(req_dim%6):(6);
382 for (
unsigned i=0; i<num; i++){
384 for (
int d=0; d<dimension; d++){
385 res[i][k6+d] = r[d] * scale[k6+d] + min[k6+d];
388 INFNANCHECK(res[i][k6+d]);
399 BIASERR(
"quasi random sequences can be only generated for spaces with "
400 "dimension less than 7. Requested dimension is "<<res.
Size());
404 BIASERR(
"invalid dimension: "<<dim);
409 for (
int i=0; i<dim; i++){
416 for (
int i=0; i<dim; i++){
429 BIASERR(
"invalid dimension: "<<dim);
434 for (
int i=0; i<dim; i++){
440 if (mres!=0)
return mres;
441 BIASASSERT(res.size()==(unsigned)num);
442 for (
unsigned i=0; i<num; i++){
443 for (
int d=0; d<dim; d++){
458 if (src<-1.0 || src>1.0)
return -1;
459 const double abs_src = fabs(src);
465 static const double sqrt_2 =
466 1.4142135623730950488016887242096980785696718753769480731766797379907325;
468 dst = sqrt(-2.0 * log(abs_src)) / sqrt_2;
469 if (src<0.0) dst = -dst;
478 double& dst1,
double& dst2)
const
480 if (src1<-1.0 || src1>1.0)
return -2;
481 if (src2<-1.0 || src2>1.0)
return -2;
482 double rsq = src1*src1 + src2*src2;
483 if (rsq>=1.0 || rsq==0.0)
return -1;
484 double fac = sqrt(-2.0 * log(rsq) / rsq);
494 #define IM 2147483647
499 #define NDIV (1+(IM-1)/NTAB)
500 #define EPS DBL_EPSILON
501 #define RNMX (1.0-EPS)
508 static long iv[NTAB];
511 if (*idum <= 0 || !iy){
512 if (-(*idum) < 1) *idum = 1;
513 else *idum = -(*idum);
514 for (j=NTAB+7; j>=0; j--){
516 *idum=IA*(*idum-k*IQ)-IR*k;
517 if (*idum <0) *idum += IM;
518 if (j<NTAB) iv[j] = *idum;
523 *idum=IA*(*idum-k*IQ)-IR*k;
524 if (*idum < 0) *idum += IM;
528 if ((temp=
float(AM*iy)) > RNMX)
return RNMX;
556 static long ix1,ix2,ix3;
563 if (*idum < 0 || iff == 0) {
565 ix1=(IC1-(*idum)) % M1;
566 ix1=(IA1*ix1+IC1) % M1;
568 ix1=(IA1*ix1+IC1) % M1;
570 for (j=1;j<=97;j++) {
571 ix1=(IA1*ix1+IC1) % M1;
572 ix2=(IA2*ix2+IC2) % M2;
573 r[j]=(ix1+ix2*RM2)*RM1;
577 ix1=(IA1*ix1+IC1) % M1;
578 ix2=(IA2*ix2+IC2) % M2;
579 ix3=(IA3*ix3+IC3) % M3;
580 j=(int)rint(1 + ((97.0*ix3)/M3));
581 if (j > 97 || j < 1) BIASERR(
"RAN1: This cannot happen. "<<j);
583 r[j]=(ix1+ix2*RM2)*RM1;
606 static long iy=0, ir[98];
610 BIASASSERT(idum!=NULL);
611 if ((*idum) < 0 || iff == 0)
614 if (((*idum)=(IC-(*idum)) % M) < 0) *idum = -(*idum);
615 for (j=1;j<=97;j++) {
616 *idum=(IA*(*idum)+IC) % M;
619 (*idum) = (IA*(*idum)+IC) % M;
622 j = (int)floor(1 + 97.0*iy / M);
625 if (j > 97 || j < 1) j = (j % 97) + 1;
627 *idum=(IA*(*idum)+IC) % M;
629 return (
double) iy/M;
645 v1=2.0*
ran1_(idum)-1.0;
646 v2=2.0*
ran1_(idum)-1.0;
649 fac=sqrt(-2.0*log(r)/r);
668 v1=2.0*
ran2_(idum)-1.0;
669 v2=2.0*
ran2_(idum)-1.0;
672 fac=sqrt(-2.0*log(r)/r);
692 v1=2.0*
ran1_(idum)-1.0;
693 v2=2.0*
ran1_(idum)-1.0;
696 fac=sqrt(-2.0*log(r)/r);
707 double &mean,
double &covariance)
const {
709 for (
unsigned int i=0; i<samples.size(); i++) {
713 if (samples.empty())
return;
714 mean /= samples.size();
715 for (
unsigned int i=0; i<samples.size(); i++) {
716 covariance += (mean - samples[i])*(mean - samples[i]);
718 covariance /= samples.size();
726 BIASASSERT(!samples.empty());
727 const int dim = samples[0].Size();
729 for (
unsigned int i=0; i<samples.size(); i++) {
732 mean /= (double)samples.
size();
736 for (
unsigned int i=0; i<samples.size(); i++) {
741 covariance *= 1.0/samples.
size();
750 unsigned long i,im,ipp;
752 static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1];
753 static unsigned long mdeg[MAXDIM+1]={0,1,2,3,3,4,4};
754 static unsigned long ip[MAXDIM+1]={0,0,1,1,2,1,4};
755 static unsigned long iv[MAXDIM*MAXBIT+1]={
756 0,1,1,1,1,1,1,3,1,3,3,1,1,5,7,7,3,3,5,15,11,5,15,13,9};
760 for (k=1;k<=MAXDIM;k++) ix[k]=0;
762 if (iv[1] != 1)
return;
763 fac=1.0/(1L << MAXBIT);
764 for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k];
765 for (k=1;k<=MAXDIM;k++) {
766 for (j=1; j<=(int)mdeg[k]; j++)
767 iu[j][k] <<= (MAXBIT-j);
769 for (j=mdeg[k]+1;j<=MAXBIT;j++) {
773 for (l=mdeg[k]-1;l>=1;l--) {
774 if (ipp & 1) i ^= iu[j-l][k];
782 for (j=1;j<=MAXBIT;j++) {
783 if (!(im & 1))
break;
787 BIASERR(
"MAXBIT too small in sobseq");
791 for (k=1;k<=min(*n,MAXDIM);k++) {
double GetNormalDistributedGSL(const double mean, const double sigma)
interface to GSL for comparison purposes
Random()
initializes vars/seed with deterministic values, NOT calling reset
double gasdevf_(long *idum)
int GetQuasiUniformDistributed(const int dim, const Vector< double > &min, const Vector< double > &max, Vector< double > &res)
On subsequent calls returns a Sobol' sequence of quasi random numbers 'uniformly distributed' in a N-...
Matrix< T > & newsize(Subscript M, Subscript N)
unsigned int Size() const
length of the vector
double GetUniformDistributed(const double min, const double max)
on succesive calls return uniform distributed random variable between min and max ...
int GetUniformDistributedInt(const int min, const int max)
get uniform distributed random variable including min/max
double gasdev3_(int *idum)
Vector< T > & newsize(Subscript N)
void SetZero()
Sets all values to zero.
double gasdev_(int *idum)
produces normal distributed random variable with unit variance and zero mean
double ran1_(int *idum)
produces uniform random variable between 0 and 1 exclusive
int EigenvalueDecomposition(T &value1, Vector2< T > &vector1, T &value2, Vector2< T > &vector2) const
Eigenvalue decomposition.
void Reset()
calls srand() with a seed generated from system call time, also initializes some other variables ...
int Uniform2Normal_(const double src, double &dst) const
transfroms a uniform disributed random variable into a normal distributed random variable ...
double GetNormalDistributed(const double mean, const double sigma)
on succesive calls return normal distributed random variable with mean and standard deviation sigma ...
void GetMeanAndCovariance(const std::vector< double > &samples, double &mean, double &covariance) const
check-function to compute mean and cov from samples
void sobseq_(int *n, double x[])
Generates a sobol sequence of n dimensional 'quasi' random vectors on succesive calls.
int GetQuasiNormalDistributed(const int dim, Vector< double > &res)
On subsequent calls returns a Sobol' sequence of quasi random numbers 'normal distributed' with mean ...