Random.cpp

/* 
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
  Multimediale Systeme der Informationsverarbeitung
  Institut fuer Informatik
  Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/


#include "Random.hh"

#ifdef WIN32
#include <Base/Common/W32Compat.hh>
#endif //_WIN32

#include <cmath>
#include <cstdlib>
#include <stdlib.h>
#include <float.h>

#ifdef BIAS_HAVE_GSL
#include <gsl/gsl_version.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#endif
using namespace BIAS;
using namespace std;

//#define RANDOM_SEED_ZERO_FIXED 

// #ifdef RANDOM_SEED_ZERO_FIXED
// #include <fstream>
// namespace BIAS { 
//   int rand(void){
//     static std::ofstream os("random.txt");
//     if (!os.good()){ BIASERR("error writing to file "); BIASABORT; }
//     static int num = 0;
//     int ran = ::rand();
//     cout << setw(5) << num << setw(15) << ran << endl << flush;
//     os << setw(5) << num << setw(15) << ran << endl << flush;
//     num++;
//     return ran;
//   }
// }
//#endif

Random::
~Random()
{}


Random::
Random()
  : Seed_(0), Idum_(0), LongIdum_(0)
{
  // to avoid nondeterministics by only instantiating a new object
  // never call srand() here
  // evers 2005-06-30
  //  Reset();
#ifdef RANDOM_SEED_ZERO_FIXED 
  Reset(0);
#endif
}

Random::
Random(const long int seed)
  : Seed_(0), Idum_(0), LongIdum_(0)
{
  Reset(seed);
}

#ifdef RANDOM_SEED_ZERO_FIXED
static bool first = true;
#endif

void Random::Reset(const long int seed)
{
#ifdef RANDOM_SEED_ZERO_FIXED 
  if (!first){
    return;
  }
  first = false;
  Seed_ = 0;
#else
  Seed_ = seed;
#endif
  //cout << "srand("<<Seed_<<")\n";
  // initialize random source
  srand((long int)Seed_);
  // initialize seed for ran1_
  Idum_=rand();
  LongIdum_ = (long)Idum_;
  //cerr << "Idum_: " << Idum_ << "\tLongIdum_: "<<LongIdum_<<endl;
}

void Random::
Reset()
{
#ifdef RANDOM_SEED_ZERO_FIXED 
  if (!first){
    return;
  }
  first = false;
  Seed_ = 0;
#else
    //long int uc;
  time_t uc;

  // get the seed (seconds since 1.1.1970 0:00)
  time(&uc);
  Seed_=(unsigned int)(uc);
#endif 
  //cout << "srand("<<Seed_<<")\n";
  // initialize random source
  srand( Seed_ );
  // initialize seed for ran1_
  Idum_=rand();
  LongIdum_ = (long)Idum_;
  //cerr << "Idum_: " << Idum_ << "\tLongIdum_: "<<LongIdum_<<endl;
}

int Random::
GetUniformDistributedInt(const int min_arg, const int max_arg)
{
  int min = (min_arg<max_arg)?(min_arg):(max_arg);
  int max = (min_arg<max_arg)?(max_arg):(min_arg);
  long int r;
#ifdef WIN32
  r = rand();
#else 
  r = random();
#endif
  double nr = (double)r / (double)(RAND_MAX);
  int val = min + (int)(nr * (double)(max-min));
  //int old_val = int(min+(double(rand())/(RAND_MAX)*(max-min+1)));
  BIASASSERT(val>=min && val<=max);
  return val;
}

unsigned int Random::
GetUniformDistributedInt(const unsigned min_arg, const unsigned max_arg)
{
  unsigned min = (min_arg<max_arg)?(min_arg):(max_arg);
  unsigned max = (min_arg<max_arg)?(max_arg):(min_arg);
//   double bin = (double)(1+max-min);
//   double res = GetUniformDistributed(0.0, 1.0) * bin;
//   //cerr << setprecision(14) << res<<" "<<(unsigned int)floor(res)+min << endl;
//   return (unsigned int)floor(res)+min;
  double dval = GetUniformDistributed((double)min, (double)max+1.0);
  unsigned val = (unsigned)(dval);
  if (val<min || val>max) val = (val % (max+1-min)) + min;
  BIASASSERT(val>=min && val<=max);
  return val;
}


void Random::
GetUniformDistributed(const double min, const double max, const int num, 
                      vector<double>& res)
{
  res.resize(num);
  double size=max-min;
  for (int i=0; i<num; i++)
    res[i]=min + ran2_(&LongIdum_) * size;
}

double Random::
GetNormalDistributedGSL(const double mean, const double sigma)
{
#ifdef BIAS_HAVE_GSL
    static bool initialized = false;
    static gsl_rng * r = NULL;
    if (!initialized){
      cout << "initializing gsl_rng ! \n";
      gsl_rng_env_setup();
      //      const gsl_rng_type *T = gsl_rng_taus2;
      const gsl_rng_type * T = gsl_rng_mt19937; 
      BIASASSERT(T);
      r = gsl_rng_alloc (T);
      BIASASSERT(r);
      initialized = true;
    }     

    return gsl_ran_gaussian(r, sigma) + mean;
  
    // this is missing
    // gsl_rng_free (r);
    //}
#else
  BIASERR("GetNormalDistributedGSL(): No GSL support in BIAS, returning "
          <<"own implementation");
  return GetNormalDistributed(mean, sigma);
#endif
}

void Random::
GetNormalDistributed(const double mean, const double sigma, 
                     const int num, std::vector<double>& res)
{
  res.resize(num);
  for (int i=0; i<num; i++)
    res[i]=mean + gasdevf_(&LongIdum_) * sigma;
}



int Random::
GetNormalDistributed(const Vector2<double>& mean, const Matrix2x2<double>& cov,
                     Vector2<double>& ran)
{
  int res=0;

  vector<Vector2<double> > vran;
  res = GetNormalDistributed(mean, cov, 1, vran);

  if (res==0) ran = vran[1];

  return res;
}

int Random::
GetNormalDistributed(const Vector2<double>& mean, 
                     const Matrix2x2<double>& cov,
                     const int num, vector<Vector2<double> >& ran)
{
  int res=0;
  Vector2<double> eigen_vec1, eigen_vec2, tmp_ran;
  double eigen_val1, eigen_val2;
  int eigen_res = 
    cov.EigenvalueDecomposition(eigen_val1,eigen_vec1, eigen_val2, eigen_vec2);
  
  //cout << "eignevalues of covariance "<<eigen_val1<<", "<<eigen_val2<<endl;
  //cout << "eignevectors of covariance "<<eigen_vec1<<", "<<eigen_vec2<<endl;
  
  Matrix2x2<double> R;
  R[0][0] = eigen_vec1[0];
  R[1][0] = eigen_vec1[1];
  R[0][1] = eigen_vec2[0];
  R[1][1] = eigen_vec2[1];
  //cout << "R: "<<R<<endl;

  vector<double> ran1, ran2;
  switch (eigen_res){
  case 2:
    //cout << "2D distribution" << endl;
    GetNormalDistributed(0.0, sqrt(eigen_val1), num, ran1);
    GetNormalDistributed(0.0, sqrt(eigen_val2), num, ran2);
    break;
  case 1:
    //cout << "pseudo 2D distribution = 1D dsitribution" << endl;
    GetNormalDistributed(0.0, sqrt(eigen_val1), num, ran1);
    ran2.resize(num, 0.0);
    break;
  case 0:
    BIASERR("not a covariance matrix ("<<eigen_res<<") "<<cov<<endl);
    //BIASASSERT(false);
    BIASABORT;
    //res = -1;
    break;
  default:
    res = eigen_res;
    break;
  }

  if (res>=0){
    ran.resize(num);
    vector<Vector2<double> >::iterator it;
    for (int i=(int)ran.size()-1; i>=0; i--){
      ran[i].Set(ran1[i], ran2[i]);
      //cout <<"---\nraw: "<< ran[i]<<endl;
      ran[i] = R * ran[i];
      //cout <<"rot: "<< ran[i]<<endl;
      ran[i] += mean;
      //cout <<"sht: "<< ran[i]<<endl;
      
    }
  }

  return res;
}


int Random::
GetQuasiUniformDistributed(const int req_dim,
                           const Vector<double>& min, 
                           const Vector<double>& max,
                           Vector<double>& res)
{
  if (req_dim>6){
    BIASERR("quasi random sequences can be only generated for spaces with "
            "dimension less than 7. Requested dimension is "<<req_dim);
    return -2;
  }
  static int dimension = -1;
  static double r[6];
  if (dimension != req_dim || req_dim<=0){
    dimension = -1;
    // initialize
    // cout << "reinitializing quasi random number sequence generator\n";
    sobseq_(&dimension, &r[0]);
    dimension = req_dim;
  }
  if (req_dim<=0){ // just initialization is requested
    return 0; 
  }
  if (min.Size()!=max.Size() || req_dim!=(int)min.Size()){
    BIASERR("argument vector size mismatch: "<<min.Size()<<" "<<max.Size()
            <<"  "<<req_dim);
    return -1;
  }
  for (int d=0; d<req_dim; d++){
    if (max[d]<min[d]) {
      BIASERR("max["<<d<<"] smaller than min["<<d<<"] ("<<max[d]
              <<" < "<<min[d]<<")");
      return -4;
    }
  }

  sobseq_(&dimension, r);
  if (res.Size() != min.Size())
    res.newsize(dimension);
  
  const Vector<double> scale = max - min;

  for (int i=0; i<dimension; i++)
    res[i] = r[i] * scale[i] + min[i];
  return 0;
}


int Random::
GetQuasiUniformDistributed(const int req_dim,
                           const Vector<double>& min, 
                           const Vector<double>& max,
                           const unsigned num, 
                           std::vector<Vector<double> >& res)
{
  for (int d=0; d<req_dim; d++){
    if (max[d]<min[d]) {
      BIASERR("max["<<d<<"] smaller than min["<<d<<"] ("<<max[d]
              <<" < "<<min[d]<<")");
      return -3;
    }
  }
  if (min.Size()!=max.Size() || req_dim!=(int)min.Size()){
    BIASERR("argument vector size mismatch: "<<min.Size()<<" "<<max.Size()
            <<"  "<<req_dim);
    return -1;
  }
  if (num==0){
    res.clear(); 
    return 0;
  }
  double r[7];
  int dimension;
  const int num_runs = req_dim/6 + (req_dim%6==0)?(0):(1);
  const Vector<double> scale = max - min;
  Vector<double> zero(req_dim, 0.0);
    res.resize(num, zero);
  for (int k=0; k<num_runs; k++){
    int k6 = k*6;
    dimension = -1;
    // initialize
    sobseq_(&dimension, &r[0]);
    dimension = (k==(num_runs-1))?(req_dim%6):(6);
    for (unsigned i=0; i<num; i++){
      sobseq_(&dimension, &r[0]);
      for (int d=0; d<dimension; d++){
        res[i][k6+d] = r[d] * scale[k6+d] + min[k6+d];
        //cout << "r[d]: "<<r[d]<<"\tres: "<<res[i][k6+d]<<endl;
        INFNANCHECK(r[d]);
        INFNANCHECK(res[i][k6+d]);
      }
    }
  }
  return 0;
}

int Random::
GetQuasiNormalDistributed(const int dim, Vector<double>& res)
{
  if (dim>6){
    BIASERR("quasi random sequences can be only generated for spaces with "
            "dimension less than 7. Requested dimension is "<<res.Size());
    return -1;
  } 
  if (dim<=0){
    BIASERR("invalid dimension: "<<dim);
    return -2;
  }

  Vector<double> uniform(dim), min(dim), max(dim);
  for (int i=0; i<dim; i++){
    min[i] = -1.0;
    max[i] = +1.0;
  }
 
  int mres = 0;
  GetQuasiUniformDistributed(dim, min, max, uniform);
  for (int i=0; i<dim; i++){
    if ((mres = Uniform2Normal_(uniform[i], res[i]))!=0){
      return mres;
    }
  }
  return 0;
}

int Random::
GetQuasiNormalDistributed(const int dim, const unsigned num,
                          std::vector<Vector<double> >& res)
{
  if (dim<=0){
    BIASERR("invalid dimension: "<<dim);
    return -2;
  }

  Vector<double> min(dim), max(dim);
  for (int i=0; i<dim; i++){
    min[i] = -1.0;
    max[i] = +1.0;
  }
 
  int mres = GetQuasiUniformDistributed(dim, min, max, num, res);
  if (mres!=0) return mres;
  BIASASSERT(res.size()==(unsigned)num);
  for (unsigned i=0; i<num; i++){
    for (int d=0; d<dim; d++){
      //cout << "uniform: "<<res[i][d];
      if ((mres = Uniform2Normal_(res[i][d], res[i][d]))!=0){
        return mres;
      }
      //cout << "\tnormal: "<<res[i][d]<<endl;
    }
  }
  return 0;
}


int Random::
Uniform2Normal_(const double src, double& dst) const
{
  if (src<-1.0 || src>1.0) return -1;
  const double abs_src = fabs(src);
  if (src==0.0) {
    //static const double sqrt_2_pi = 
    //  2.5066282746310005024157652848110452530069867406099383166299235763422937;
    dst = 0.0;
  } else {
    static const double sqrt_2 = 
      1.4142135623730950488016887242096980785696718753769480731766797379907325;

    dst = sqrt(-2.0 * log(abs_src)) / sqrt_2;
    if (src<0.0) dst = -dst;
  }
  INFNANCHECK(dst);
  return 0;
}


int Random::
Uniform2Normal_(const double src1, const double src2, 
                double& dst1, double& dst2) const
{
  if (src1<-1.0 || src1>1.0) return -2;
  if (src2<-1.0 || src2>1.0) return -2;
  double rsq = src1*src1 + src2*src2;
  if (rsq>=1.0 || rsq==0.0) return -1;
  double fac = sqrt(-2.0 * log(rsq) / rsq);
  dst1 = src1 * fac;
  dst2 = src2 * fac;
  INFNANCHECK(dst1);
  INFNANCHECK(dst2);
  return 0;
}

// from www.nr.com copied from pdf
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS DBL_EPSILON
#define RNMX (1.0-EPS)
double Random::
ran1_(int *idum)
{
  int j;
  long k;
  static long iy=0;
  static long iv[NTAB];
  float temp;

  if (*idum <= 0 || !iy){
    if (-(*idum) < 1) *idum = 1;
    else *idum = -(*idum);
    for (j=NTAB+7; j>=0; j--){
      k=(*idum)/IQ;
      *idum=IA*(*idum-k*IQ)-IR*k;
      if (*idum <0) *idum += IM;
      if (j<NTAB) iv[j] = *idum;
    }
    iy=iv[0];
  }
  k=(*idum)/IQ;
  *idum=IA*(*idum-k*IQ)-IR*k;
  if (*idum < 0) *idum += IM;
  j=iy/NDIV;
  iy=iv[j];
  iv[j] = *idum;
  if ((temp=float(AM*iy)) > RNMX) return RNMX;
  else return temp;
}
#undef IA 
#undef IM 
#undef AM 
#undef IQ 
#undef IR 
#undef NTAB 
#undef NDIV 
#undef EPS 
#undef RNMX 

//from D_RAN1.c
#define M1 259200
#define IA1 7141
#define IC1 54773
#define RM1 (1.0/M1)
#define M2 134456
#define IA2 8121
#define IC2 28411
#define RM2 (1.0/M2)
#define M3 243000
#define IA3 4561
#define IC3 51349
double Random::
ran1c_(int *idum)
{
  static long ix1,ix2,ix3;
  static double r[98];
  double temp;
  static int iff=0;
  int j;
  //void nrerror();

  if (*idum < 0 || iff == 0) {
    iff=1;
    ix1=(IC1-(*idum)) % M1;
    ix1=(IA1*ix1+IC1) % M1;
    ix2=ix1 % M2;
    ix1=(IA1*ix1+IC1) % M1;
    ix3=ix1 % M3;
    for (j=1;j<=97;j++) {
      ix1=(IA1*ix1+IC1) % M1;
      ix2=(IA2*ix2+IC2) % M2;
      r[j]=(ix1+ix2*RM2)*RM1;
    }
    *idum=1;
  }
  ix1=(IA1*ix1+IC1) % M1;
  ix2=(IA2*ix2+IC2) % M2;
  ix3=(IA3*ix3+IC3) % M3;
  j=(int)rint(1 + ((97.0*ix3)/M3));
  if (j > 97 || j < 1) BIASERR("RAN1: This cannot happen. "<<j);
  temp=r[j];
  r[j]=(ix1+ix2*RM2)*RM1;
  return temp;
}
#undef M1
#undef IA1
#undef IC1
#undef RM1
#undef M2
#undef IA2
#undef IC2
#undef RM2
#undef M3
#undef IA3
#undef IC3

// from D_RAN2.c
#define M 714025
#define IA 1366
#define IC 150889
double Random::
ran2_(long *idum)
{
  //cout << "idum "<<*idum<<endl;
  static long iy=0, ir[98];
  static int iff=0;
  int j;

  BIASASSERT(idum!=NULL);
  if ((*idum) < 0 || iff == 0) 
  {
      iff=1;
      if (((*idum)=(IC-(*idum)) % M) < 0) *idum = -(*idum);
      for (j=1;j<=97;j++) {
          *idum=(IA*(*idum)+IC) % M;
          ir[j]=(*idum);
      }
      (*idum) = (IA*(*idum)+IC) % M;
      iy=(*idum);
  }
  j = (int)floor(1 + 97.0*iy / M);
  //if (j > 97 || j < 1) BIASERR("RAN2: This should not happen (j="<<j<<")!");
  if (j < 0) j = -j;
  if (j > 97 || j < 1) j = (j % 97) + 1;
  iy=ir[j];
  *idum=(IA*(*idum)+IC) % M;
  ir[j]=(*idum);
  return (double) iy/M;
}
#undef M
#undef IA
#undef IC

// from D_GASDEV.c
double Random::
gasdev_(int *idum)
{
  static int iset=0;
  static double gset;
  double fac,r,v1,v2;

  if  (iset == 0) {
    do {
      v1=2.0*ran1_(idum)-1.0;
      v2=2.0*ran1_(idum)-1.0;
      r=v1*v1+v2*v2;
    } while (r >= 1.0);
    fac=sqrt(-2.0*log(r)/r);
    gset=v1*fac;
    iset=1;
    return v2*fac;
  } else {
    iset=0;
    return gset;
  }
}

double Random::
gasdevf_(long *idum)
{
  static int iset=0;
  static double gset;
  double fac,r,v1,v2;

  if  (iset == 0) {
    do {
      v1=2.0*ran2_(idum)-1.0;
      v2=2.0*ran2_(idum)-1.0;
      r=v1*v1+v2*v2;
    } while (r >= 1.0);
    fac=sqrt(-2.0*log(r)/r);
    gset=v1*fac;
    iset=1;
    return v2*fac;
  } else {
    iset=0;
    return gset;
  }
}

// from makenois.c
double Random::
gasdev3_(int *idum)
{
  static int iset=1;
  static double gset;
  double fac,r,v1,v2;

  if (iset) {
    do {
      v1=2.0*ran1_(idum)-1.0;
      v2=2.0*ran1_(idum)-1.0;
      r=v1*v1+v2*v2;
    } while (r >= 1.0);
    fac=sqrt(-2.0*log(r)/r);
    gset = v1*fac;
    iset = 0;
    return v2*fac;
  } else {
    iset = 1;
    return gset;
  }
}

void Random::GetMeanAndCovariance(const std::vector<double>& samples,
                                  double &mean, double &covariance) const {
  mean = 0;
  for (unsigned int i=0; i<samples.size(); i++) {
    mean += samples[i];
  }
  covariance = 0;
  if (samples.empty()) return;
  mean /= samples.size();
  for (unsigned int i=0; i<samples.size(); i++) {
    covariance += (mean - samples[i])*(mean - samples[i]);
  }
  covariance /= samples.size();
}


void Random::GetMeanAndCovariance(const std::vector<BIAS::Vector<double>
                                  >&samples,
                                  BIAS::Vector<double> &mean, 
                                  BIAS::Matrix<double> &covariance) const {
  BIASASSERT(!samples.empty());
  const int dim = samples[0].Size();
  mean.newsize(dim);
  for (unsigned int i=0; i<samples.size(); i++) {
    mean += samples[i];
  }
  mean /= (double)samples.size();
 
  covariance.newsize(dim,dim);
  covariance.SetZero();
  for (unsigned int i=0; i<samples.size(); i++) {
    covariance += 
      BIAS::Vector<double>(mean - samples[i]).OuterProduct((mean -
                                                            samples[i]));
  }
  covariance *= 1.0/samples.size();
}

#define MAXBIT 30
#define MAXDIM 6
void Random::
sobseq_(int *n, double x[])
{
    int j,k,l;
    unsigned long i,im,ipp;
    static double fac;
    static unsigned long in,ix[MAXDIM+1],*iu[MAXBIT+1];
    static unsigned long mdeg[MAXDIM+1]={0,1,2,3,3,4,4};
    static unsigned long ip[MAXDIM+1]={0,0,1,1,2,1,4};
    static unsigned long iv[MAXDIM*MAXBIT+1]={
        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};

    if (*n < 0) {
    //cout << "initializing\n";
        for (k=1;k<=MAXDIM;k++) ix[k]=0;
        in=0;
        if (iv[1] != 1) return;
        fac=1.0/(1L << MAXBIT);
        for (j=1,k=0;j<=MAXBIT;j++,k+=MAXDIM) iu[j] = &iv[k];
        for (k=1;k<=MAXDIM;k++) {
            for (j=1; j<=(int)mdeg[k]; j++) 
        iu[j][k] <<= (MAXBIT-j);

            for (j=mdeg[k]+1;j<=MAXBIT;j++) {
                ipp=ip[k];
                i=iu[j-mdeg[k]][k];
                i ^= (i >> mdeg[k]);
                for (l=mdeg[k]-1;l>=1;l--) {
                    if (ipp & 1) i ^= iu[j-l][k];
                    ipp >>= 1;
                }
                iu[j][k]=i;
            }
        }
    } else {
        im=in++;
        for (j=1;j<=MAXBIT;j++) {
            if (!(im & 1)) break;
            im >>= 1;
        }
        if (j > MAXBIT) {
      BIASERR("MAXBIT too small in sobseq");
      BIASABORT;
    }
        im=(j-1)*MAXDIM;
        for (k=1;k<=min(*n,MAXDIM);k++) {
            ix[k] ^= iv[im+k];
      //cout << "x["<<k-1<<"] = "<<x[k-1]<<endl;
            x[k-1]=ix[k]*fac;
        }
    }
}
#undef MAXBIT
#undef MAXDIM
#undef NRANSI