Random.hh

/* 
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
*/


#ifndef __BIAS_Random_hh__
#define __BIAS_Random_hh__

#include <Base/Common/BIASpragmaStart.hh>

#include <ctime>
#include <vector>

#include <Base/Debug/Error.hh>
#include <Base/Debug/Debug.hh>
#include "Vector2.hh"
#include "Matrix2x2.hh"


namespace BIAS{
/** @class Random
    @ingroup g_mathalgo
    @brief class for producing random numbers from different distributions
       
    IMPORTANT: Problems arise when multiple of these instances are used
         since every call to Reset() results in a call to srand() (from 
         stdlib) with a specific seed. Calling srand() is a global action 
         that will affect all other instances as well 
         
    @author woelk 2002 */
  class BIASMathBase_EXPORT Random : public Debug
  {  
  public:
    /** initializes vars/seed with deterministic values, NOT calling reset */
    Random();

    /** calls Reset(seed) and therefore srand() */
    Random(const long int seed);
    
    ~Random();

    /** calls srand() with a seed generated from system call time, 
       also initializes some other variables */
    void Reset();

    /** calls srand with seed, also initializes some other variables */
    void Reset(const long int seed);

    /** on succesive calls return normal distributed random variable
        with mean and standard deviation sigma */
    inline double GetNormalDistributed(const double mean, const double sigma)
    { return (mean + (gasdevf_(&LongIdum_) * sigma)); };

    /** interface to GSL for comparison purposes */
    double GetNormalDistributedGSL(const double mean, const double sigma);

    /** calculates num normal ditributed random variables and stores
        them in res */
    void GetNormalDistributed(const double mean, const double sigma, 
                              const int num, std::vector<double>& res); 
    
    /** on succesive calls return uniform distributed random variable
        between min and max */
    inline double GetUniformDistributed(const double min, const double max)
    { return (min + (ran2_(&LongIdum_) * (max-min))); };

    /** calculates num uniform ditributed random variables and stores
        them in res */
    void GetUniformDistributed(const double min, const double max, 
                               const int num, std::vector<double>& res); 

    /** get uniform distributed random variable 
        including min/max */
    int GetUniformDistributedInt(const int min, const int max);

    /** get uniform distributed random variable 
        including min/max */
    unsigned GetUniformDistributedInt(const unsigned min, const unsigned max);

    /** See GetNormalDistributed().
        @author woelk 07/2005 */
    int GetNormalDistributed(const Vector2<double>& mean, 
                             const Matrix2x2<double>& cov,
                             Vector2<double>& ran);

    /** Generates a vector of length num containing 2D (dependent) normal 
        distributed random variables. The 2D distribution is described by the 
        covariance matrix cov.

        @returns 0 on success , negative on error
        @author woelk 07/2005 */
    int GetNormalDistributed(const Vector2<double>& mean, 
                             const Matrix2x2<double>& cov,
                             const int num, 
                             std::vector<Vector2<double> >& ran);
    
    /** @brief On subsequent calls returns a Sobol' sequence of quasi random 
        numbers 'uniformly distributed' in a N-dimensional cube between min 
        and max. The dimension of the space must not exceed 6. 
        min[i] >= max[i] must hold for all i
        @author woelk 11/2006 */
    int GetQuasiUniformDistributed(const int dim,
                                   const Vector<double>& min, 
                                   const Vector<double>& max,
                                   Vector<double>& res);

    /** @brief Returns a Sobol' sequence of num quasi random numbers 
        'uniformly distributed' in the N-dimensional cube between min and max.
        min[i] >= max[i] must hold for all i
        @author woelk 11/2006 */
    int GetQuasiUniformDistributed(const int dim,
                                   const Vector<double>& min, 
                                   const Vector<double>& max,
                                   const unsigned num, 
                                   std::vector<Vector<double> >& res);

    /** @brief On subsequent calls returns a Sobol' sequence of quasi random 
        numbers 'normal distributed' with mean 0 and Covariance matrix identity
        Input:  dim  - dimension of requested sequqnce, in range [1,6]
        Output: res  - quasi random vector
        See GetQuasiUniformDistributed() and sobseq_() for details.
        @author woelk 11/2006 */
    int GetQuasiNormalDistributed(const int dim, Vector<double>& res);
    
    /** @brief returns a quasi random seuqence of 'normal dsitributed'
        numbres.
        Input:  dim  - dimension of requested sequqnce, in range [1,6]
                num  - length of requested sequence
        Output: res  - quasi random vector sequence
        See GetQuasiUniformDistributed() and sobseq_() for details.
        @author woelk 11/2006 */
    int GetQuasiNormalDistributed(const int dim, const unsigned num,
                                  std::vector<Vector<double> >& res);

    
    /// old deprectade interfaces
    inline void GetNormalDistributed(const double mean, const double sigma,
                                     const int num, double* res)
    { 
      BIASERR("you are using a deprecated interface");
      std::vector<double> mres;
      GetNormalDistributed(mean, sigma, num, mres);
      BIASASSERT((int)mres.size()==num);
      for (int i=0; i<num; i++) res[i] = mres[i];
    }
    inline void GetUniformDistributed(const double min, const double max,
                                      const int num, double* res)
    { 
      BIASERR("you are using a deprecated interface");
      std::vector<double> mres;
      GetUniformDistributed(min, max, num, mres);
      BIASASSERT((int)mres.size()==num);
      for (int i=0; i<num; i++) res[i] = mres[i];
    }
    
    /** @brief check-function to compute mean and cov from samples */
    void GetMeanAndCovariance(const std::vector<double>& samples,
                              double &mean, double &covariance) const;

    /** @brief check-function to compute mean and cov from samples */
    void GetMeanAndCovariance(const std::vector<BIAS::Vector<double> >&samples,
                              BIAS::Vector<double> &mean, 
                              BIAS::Matrix<double> &covariance) const;

  protected:
    unsigned int Seed_;
    int Idum_;
    long LongIdum_;

    /// produces uniform random variable between 0 and 1 exclusive
    double ran1_(int *idum); // from nr pdf
    double ran1c_(int *idum); // from D_RAN1.c
    double ran2_(long *idum); // from D_RAN2.c
    /// produces normal distributed random variable with unit variance and
    /// zero mean
    double gasdev_(int *idum);
    // from makenois.c
    double gasdev3_(int *idum);
    // selbst gebastelt
    double gasdevf_(long *idum);

    /** Generates a sobol sequence of n dimensional 'quasi' random vectors
        on succesive calls.
        When n is negative, internally intializes a set of MAXBIT=6 direction 
        numbers for each MAXDIM different Sobol' sequences. When n is positive
        (but <=MAXDIM), returns as thevector x[0..n-1] the next values from n 
        of theses sequneces. (n must not be changed between initializations. )
        See Press, "Numerical Recipes in C", Chapter 7.7, p. 312
        http://www.nrbook.com/a/bookcpdf.html
        @author woelk 11/2006 */
    void sobseq_(int *n, double x[]);


    /** @brief transfroms a uniform disributed random variable into a 
        normal distributed random variable

        The algiorithm is based on the inversion Principle, see 
        L. Devroye "Non-Uniform Random Variate Generation", 1986, Springer, 
        chapter II.2 for details

        Input:  src  - uniform deviat in range [-1.0, +1.0]
        Output: dst  - normal distributed deviate
 
        returns  0  on success
                -1  when src is out of bounds

        @author woelk 11/2006 */
    int Uniform2Normal_(const double src, double& dst) const;

    /** @brief Box-Muller transform generating 2 normal distributed
        random variables from two uniform distributed random variables.

        Input:  src1, src2  -  uniform deviate in range [-1, +1] with 
                               src1*src1+src2*src2<1.0 and
                               src1*src1+src2*src2!=0.0
        Output: dst1, dst2  -  normal deviates with sigma 1

        returns  0 on success, 
                -2 when src1 or src2 out of bounds
                -1 either src1*src1+src2*src2>=1.0 or 
                   rc1*src1+src2*src2 == 0.0

        See Press, "Numerical Recipes in C", Chapter 7, p. 289
        http://www.nrbook.com/a/bookcpdf.html 
        for details.
        @author woelk 11/2006 */
     int Uniform2Normal_(const double src1, const double src2, 
                         double& dst1, double& dst2) const;
  };


} // namespace

#include <Base/Common/BIASpragmaEnd.hh>

#endif // __Random_hh__