random.h

/**************************************************************************
 * Copyright (c) 2002-2011 T. M. Murali                                   *
 *                                                                        *
 * This file is part of Biorithm.                                         *
 *                                                                        *
 * Biorithm is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU General Public License as published by   *
 * the Free Software Foundation, either version 3 of the License, or      *
 * (at your option) any later version.                                    *
 *                                                                        *
 * Biorithm 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 General Public License for more details.                           *
 *                                                                        *
 * You should have received a copy of the GNU General Public License      *
 * along with Biorithm.  If not, see <http://www.gnu.org/licenses/>.      *
 *                                                                        *
 **************************************************************************/

#ifndef _RANDOM_H
#define _RANDOM_H

#include <vector>
#include <time.h>
#include <math.h>

#include "constants.h"
// Test for GCC >= 4.3. from http://gcc.gnu.org/onlinedocs/cpp/Common-Predefined-Macros.html. See definition of GCC_VERSION macro in libutil/constants.h
#if GCC_VERSION >= 40300
// see http://gcc.gnu.org/gcc-4.3/changes.html
#include <tr1/unordered_set>
using namespace std::tr1;
#else
#include <ext/hash_map>
using namespace __gnu_cxx;
#define unordered_set hash_set;
#define unordered_map hash_map;
#endif

using namespace std;

#include "global.h"
#include "Mersennetwister.h"

// have to specialise template in the same namespace as its definition.

#ifdef __GNUC__
#if __GNUC__ < 3
namespace Sgi { using ::hash_map; }; // inherit globals
#else
#if __GNUC_MINOR__ == 0
namespace Sgi = std;               // GCC 3.0
#else
namespace Sgi = ::__gnu_cxx;       // GCC 3.1 and later
#endif
#endif
#else      // ...  there are other compilers, right?
namespace Sgi = std;
#endif

#if (GCC_VERSION < 40300)
// i need to define some hash templates.
namespace __gnu_cxx
{
template <> struct hash< double >
{
  // this function is not great but all the hash functions in
  // stl_hash_fun.h that take numbers as parameters simply return
  // their argument so doing so for double may be ok.
  size_t operator()(const double& num) const { return static_cast< size_t >(num); }
};
}
#endif // GCC_VERSION

// dummy class just so that i can define generators for every type i
// want by using partial specialisation and calling the appropriate
// MTRand function.
template <class MyNT> class MyRand : public MTRand
{
public:
  MyRand()
      : MTRand()
    {}
  virtual ~MyRand()
    {}
  
  MyNT operator()();
  MyNT operator()(MyNT num);
  MyNT operator()(MyNT low, MyNT high);
};

template <> int MyRand< int >::operator()();
template <> int MyRand< int >::operator()(int num);
template <> int MyRand< int >::operator()(int low, int high);
template <> double MyRand< double >::operator()();
template <> double MyRand< double >::operator()(double num);
template <> double MyRand< double >::operator()(double low, double high);

enum MyUniqueRandomStatus {MyUniqueRandomStatus_OK, MyUniqueRandomStatus_NONE};


// make this templated? if so, define a function that returns a random
// number of that type.
template < class MyNT > class MyUniqueRandom
{
public:
  // can only be called with l and h. makes them set for the lifetime
  // of the instance. therefore, i dont have to worry about dealing
  // with the user changing low and high in the middle of the lifetime
  // of the instance.
  MyUniqueRandom(MyNT l, MyNT h)
      : low(l), high(h), generator(), unique()
    {
    }
//   MyUniqueRandom(MyUniqueRandom& r)
//       : low(r.low), high(r.high), generator(), unique(r.unique)
//     {
//     }
//   MyUniqueRandom& operator=(const MyUniqueRandom& rhs)
//     {
//       if (this != &rhs)
//         {
//           low = rhs.low;
//           high = rhs.high;
//           unique = rhs.unique;
//         }
//       return(*this);
//     }
  
  virtual ~MyUniqueRandom()
    {}

  // forget all the random numbers you have generated and stored so far.
  void reset()
    {
      // resize() doesn't decrease size of hash_maps.
      unique.clear();
    }

  // eventhough the compiler (gcc) will not inline the generate()
  // functions, i have to put their defns. in the header file.
  // otherwise, i have to use the export keyword either before the
  // MyUniqueRandom class or before the defns. of these functions in
  // random.C. however, gcc 2.96 does not support the export keyword.
  // see sec 18 and sec 16.8.2, in particular, of effective c++.
  MyUniqueRandomStatus generate(MyNT& num)
    {
      if (allGenerated(unique.size()))
        {
          cerr << "ERROR! All unique random numbers between " << low << " and " << high
               << " have been generated." << endl;
          return(MyUniqueRandomStatus_NONE);
        }
      
      while (true)
        {
          num = generator(low, high);
          if (unique.end() == unique.find(num))
            {
              unique.insert(num);
              return(MyUniqueRandomStatus_OK);
        }
        }
    }

  // need to duplicate code in MyUniqueRandomStatus generate(MyNT&
  // num) because i can't specify a default value for numsToAvoid. if
  // i pass a pointer, then i have to check if the pointer is NULL at
  // two places, which is costly. if i pass a hash set of size 0 then
  // a function that just wants a random number has to pass an empty
  // hash set to generate(), which seems wrong.
  //
  // generate 1 random number between this->low and this->high. the
  // random number should not be in numsToAvoid. return it in 1.
  MyUniqueRandomStatus generate(
    MyNT& num,
#if GCC_VERSION >= 40300
    const unordered_set< MyNT, hash< MyNT > >& numsToAvoid)
#else
    const hash_set< MyNT, hash< MyNT > >& numsToAvoid)
#endif
    {
      if (allGenerated(unique.size() + numsToAvoid.size()))
        {
          cerr << "ERROR! All unique random numbers between " << low << " and " << high
               << " have been generated." << endl;
          return(MyUniqueRandomStatus_NONE);
        }
      
      while (true)
        {
          num = generator(low, high);
          if ((unique.end() == unique.find(num)) &&
              (numsToAvoid.end() == numsToAvoid.find(num)))
            {
              unique.insert(num);
              return(MyUniqueRandomStatus_OK);
            }
        }
    }


  // generate count random numbers between this->low and this->high. the
  // random numbers should not be in numsToAvoid. return them in nums.
  MyUniqueRandomStatus generate(
    int count,
#if GCC_VERSION >= 40300
    const unordered_set< MyNT, hash< MyNT > >& numsToAvoid,
#else
    const hash_set< MyNT, hash< MyNT > >& numsToAvoid,
#endif
    vector< MyNT >& nums)
    {
      if (!canGenerate(count + numsToAvoid.size()))
        {
          cerr << "ERROR! MyUniqueRandom::generate(): cannot generate " <<
            count + numsToAvoid.size() << " unique random numbers between "
               << low << " and " << high << "." << endl;
          return(MyUniqueRandomStatus_NONE);
    }
      // TODO: DESIGN. making localUnique static for the same reason
      // as making randomIncides static in MyPointSet::generate(int
      // num, vector< int >& indices, bool forbidChosen = true).
#if GCC_VERSION >= 40300
      static unordered_set< MyNT, hash< MyNT > > localUnique;
#else
      static hash_set< MyNT, hash< MyNT > > localUnique;
#endif
      MyNT num;
      while (0 < count)
        {
          num = generator(low, high);
          // TODO: BUG. what if i have already generated all possible
          // numbers between low and high, especially when MyNT = int? i
          // go into a infinite loop here.
          if ((localUnique.end() == localUnique.find(num)) &&
              (numsToAvoid.end() == numsToAvoid.find(num)))
            {
              localUnique.insert(num);
              nums[--count] = num;
            }
        }
      localUnique.clear();
      return(MyUniqueRandomStatus_OK);
    }

  

  // generate count random numbers from those in numsToUse. i am
  // guaranteeing that each number is between low and high. the random
  // numbers should not be in numsToAvoid. return them in nums.
  MyUniqueRandomStatus generate(
    int count, 
    const vector< int >& numsToUse,
#if GCC_VERSION >= 40300
    const unordered_set< MyNT, hash< MyNT > >& numsToAvoid,
#else
    const hash_set< MyNT, hash< MyNT > >& numsToAvoid,
#endif
    vector< MyNT >& nums)
    {
         // i want a number in the range [0, numsToUse.size() - 1]. 
         if (!canGenerate(0, numsToUse.size() - 1, 
                                   count + numsToAvoid.size()))
        {
          cerr << "ERROR! MyUniqueRandom::generate(): cannot generate " <<
            count + numsToAvoid.size() << " unique random numbers between "
               << 0 << " and " << numsToUse.size() - 1 << "." << endl;
          return(MyUniqueRandomStatus_NONE);
           }
      // TODO: DESIGN. making localUnique static for the same reason
      // as making randomIncides static in MyPointSet::generate(int
      // num, vector< int >& indices, bool forbidChosen = true).
#if GCC_VERSION >= 40300
      static unordered_set< MyNT, hash< MyNT > > localUnique;
#else
      static hash_set< MyNT, hash< MyNT > > localUnique;
#endif
      MyNT num, randomNum;
      while (0 < count)
        {
                //          num = generator(low, high);
          randomNum = generator(0, numsToUse.size() - 1);
          num = numsToUse[randomNum];
          // TODO: BUG. what if i have already generated all possible
          // numbers between low and high, especially when MyNT = int? i
          // go into a infinite loop here.
                //
                // look like i have to first cheque if the hash is empty.
          if ((localUnique.empty() ||
                        (localUnique.end() == localUnique.find(num))) &&
                    (numsToAvoid.empty() ||
                        (numsToAvoid.end() == numsToAvoid.find(num))))
            {
              localUnique.insert(num);
              nums[--count] = num;
            }
        }
      localUnique.clear();
      return(MyUniqueRandomStatus_OK);
    }
  
private:
  int numPossibleUnique() const;
  bool allGenerated(int num) const;
  bool canGenerate(int num) const;
  // don't use MyUniqueRandom::[low,high] but the ones passed.
  bool canGenerate(int newlow, int newhigh, int num) const;
  
private:
  MyNT low;
  MyNT high;
  MyRand< MyNT > generator;
#if GCC_VERSION >= 40300
  unordered_set< MyNT > unique;    
#else
  hash_set< MyNT > unique;  
#endif
};



#endif // _RANDOM_H