coconut_random.h

Go to the documentation of this file.

// Random for COCONUT implementation -*- C++ -*-

// $Id: coconut_random.h 862 2009-01-07 15:30:51Z dferi $
// Copyright (C) 2003-2008 Hermann Schichl, Mihaly Csaba Markot
//
// This file is part of the COCONUT API.  This library
// is free software; you can redistribute it and/or modify it under the
// terms of the Library GNU General Public License as published by the
// Free Software Foundation; either version 2, or (at your option)
// any later version.

// This library 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
// Library GNU General Public License for more details.

// As a special exception, you may use this file as part of a free software
// library without restriction.  Specifically, if other files instantiate
// templates or use macros or inline functions from this file, or you compile
// this file and link it with other files to produce an executable, this
// file does not by itself cause the resulting executable to be covered by
// the Library GNU General Public License.  This exception does not however
// invalidate any other reasons why the executable file might be covered by
// the Library GNU General Public License.

#ifndef __COCONUT_RANDOM_H__
#define __COCONUT_RANDOM_H__

#include <cmath>
#include <cstdlib>
#include <api_debug.h>
#include <interval.h>
#if defined(COCONUT_WINDOWS)
#include <windows_utils.h>
#endif

namespace coco {

#if defined(COCONUT_WINDOWS)
#define coconut_random  rand
#else
#define coconut_random  random
#endif

#if defined(COCONUT_WINDOWS)
#define coconut_seed(A) srand(A)
#else
#define coconut_seed(A) srandom(A)
#endif

#define COCONUT_RAND_MAX RAND_MAX

#define COCONUT_RRAND_MIN -1.e08

#define COCONUT_RRAND_MAX +1.e08

#if DEBUG_RANDOM
#define INIT_SEED    1
#else
#define INIT_SEED    coconut_random()
#endif

#define coconut_init_random() coconut_seed(INIT_SEED)

typedef long int rand_t;

inline double d_random()
{
#if (DEBUG_RANDOM > 3)
  double help = (coconut_random()+0.)/(COCONUT_RAND_MAX+0.);
  std::cout << "Random generated: " << help << std::endl;
  return help;
#else
  return (coconut_random()+0.)/(COCONUT_RAND_MAX+0.);
#endif
}

#define COCONUT_RRAND_MIN_BETA 0.01

#define COCONUT_RRAND_MIN_ALPHA 0.5

#define COCONUT_RRAND_MAX_ALPHA 0.99

// Evaluates h(t) discussed below.
inline double r_random_h_eval(double t, double b, double a)
{
  double ret;
  if (t>b)
    ret = (-a*b + t) / ((1./a-2.)*b + t);
  else if (t<-b)
    ret = (a*b + t) / ((1./a-2.)*b - t);
  else
    ret = a*t/b;
  return ret;
}

// Evaluates the inverse of h(t) discussed below.
inline double r_random_hinv_eval(double x, double b, double a)
{
  double ret;
  if (x>a)
    ret = b*(a + (1./a-2.)*x) / (1-x);
  else if (x<-a)
    ret = b*(-a + (1./a-2.)*x) / (1+x);
  else
    ret = b*x/a;
  return ret;
}

inline double r_random(double l, double u, double beta = 10.0,
                       double alpha = 0.9)
{
  l = std::max(l,COCONUT_RRAND_MIN);
  u = std::min(u,COCONUT_RRAND_MAX);
  if(l == u) return l;

  // adjust alpha and beta if necessary
  alpha = std::max(alpha,COCONUT_RRAND_MIN_ALPHA);
  alpha = std::min(alpha,COCONUT_RRAND_MAX_ALPHA);
  beta = std::max(beta,COCONUT_RRAND_MIN_BETA);

  // evaluate h(l) and h(u)
  double hl = r_random_h_eval(l,beta,alpha);
  double hu = r_random_h_eval(u,beta,alpha);
  
  // obtain a random number (form uniform distribution) in [hl,hu]
  double r = d_random()*(hu-hl) + hl;
  
  // evaluate h^{-1}(r)
  double result = r_random_hinv_eval(r,beta,alpha);
  
  // adjust result if it is outside [l,u] due to bad rounding
  result = std::max(result,l);
  result = std::min(result,u);
  
#if (DEBUG_RANDOM > 3)
  std::cout << "r_random: " << l << ", " << u << ", " << rndnum << std::endl;
  std::cout << "r_random: result: " << result << std::endl;
#endif
  return result;
}

inline double r_random(const interval& _i, double beta = 10.0,
                       double alpha = 0.9)
{
  return r_random(_i.inf(), _i.sup(), beta, alpha);
}

inline interval i_random(double l, double r, double beta = 10.0,
                         double alpha = 0.9)
{
  double d(r_random(l,r,beta,alpha)), e(r_random(l,r,beta,alpha));
  return interval(std::min(d,e),std::max(d,e));
}

inline interval i_random(const interval& _i, double beta = 10.0,
                         double alpha = 0.9)
{
  return i_random(_i.inf(), _i.sup(), beta, alpha);
}

inline interval i_random()
{
  return i_random(0., 1., true);
}

struct coconut_random_f
{
  double beta;
  double alpha;
  
  interval operator()(const interval& i) const
    { return i_random(i, beta, alpha); }

  double operator()(double l, double u) const
    { return r_random(l, u, beta, alpha); }

  double operator()(double u) const
    { return r_random(0., u, beta, alpha); }

  double operator()() const
    { return d_random(); }

  int operator()(int n) const
    { double r(d_random()); return (int) lrint(floor(r*n)); }

  int operator()(long int n) const
    { double r(d_random()); return (int) lrint(floor(r*n)); }

  size_t operator()(size_t n) const
    { double r(d_random()); return (size_t) lrint(floor(r*n)); }

  coconut_random_f(double _b = 10.0, double _a = 0.9) : beta(_b), alpha(_a) {}
};

} // namespace coco

#endif /* __COCONUT_RANDOM_H__ */

---