Random.h

/*
 *   Repast for High Performance Computing (Repast HPC)
 *
 *   Copyright (c) 2010 Argonne National Laboratory
 *   All rights reserved.
 *
 *   Redistribution and use in source and binary forms, with
 *   or without modification, are permitted provided that the following
 *   conditions are met:
 *
 *     Redistributions of source code must retain the above copyright notice,
 *     this list of conditions and the following disclaimer.
 *
 *     Redistributions in binary form must reproduce the above copyright notice,
 *     this list of conditions and the following disclaimer in the documentation
 *     and/or other materials provided with the distribution.
 *
 *     Neither the name of the Argonne National Laboratory nor the names of its
 *     contributors may be used to endorse or promote products derived from
 *     this software without specific prior written permission.
 *
 *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 *   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
 *   PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE TRUSTEES OR
 *   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 *   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 *   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 *   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 *   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 *   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
 *   EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 *
 *  Random.h
 *
 *  Created on: Jun 9, 2009
 *      Author: nick
 */
 
#ifndef RANDOM_H_
#define RANDOM_H_
 
#include <vector>
#include <set>
#include <map>
#include <string>
 
 
#include <boost/random/variate_generator.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/triangle_distribution.hpp>
#include <boost/random/cauchy_distribution.hpp>
#include <boost/random/exponential_distribution.hpp>
#include <boost/random/geometric_distribution.hpp>
#include <boost/random/normal_distribution.hpp>
#include <boost/random/lognormal_distribution.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/cstdint.hpp>
 
 
namespace repast {
 
typedef boost::variate_generator<boost::mt19937&, boost::uniform_real<> > _RealUniformGenerator;
typedef boost::variate_generator<boost::mt19937&, boost::uniform_int<> > _IntUniformGenerator;
 
typedef boost::variate_generator<boost::mt19937&, boost::triangle_distribution<> > _TriangleGenerator;
typedef boost::variate_generator<boost::mt19937&, boost::cauchy_distribution<> > _CauchyGenerator;
typedef boost::variate_generator<boost::mt19937&, boost::exponential_distribution<> > _ExponentialGenerator;
typedef boost::variate_generator<boost::mt19937&, boost::geometric_distribution<boost::uniform_real<> > >
        _GeometricGenerator;
typedef boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > _NormalGenerator;
typedef boost::variate_generator<boost::mt19937&, boost::lognormal_distribution<> > _LogNormalGenerator;
 
class NumberGenerator {
 
public:
    virtual ~NumberGenerator() {
    }
 
    virtual double next() = 0;
};
 
template<typename T>
class DefaultNumberGenerator: public NumberGenerator {
 
private:
    T gen;
 
public:
    DefaultNumberGenerator(T generator);
    double next();
};
 
template<typename T>
DefaultNumberGenerator<T>::DefaultNumberGenerator(T generator) :
    gen(generator) {
}
 
template<typename T>
double DefaultNumberGenerator<T>::next() {
    return gen();
}
 
typedef DefaultNumberGenerator<_IntUniformGenerator> IntUniformGenerator;
typedef DefaultNumberGenerator<_RealUniformGenerator> DoubleUniformGenerator;
typedef DefaultNumberGenerator<_TriangleGenerator> TriangleGenerator;
typedef DefaultNumberGenerator<_CauchyGenerator> CauchyGenerator;
typedef DefaultNumberGenerator<_ExponentialGenerator> ExponentialGenerator;
typedef DefaultNumberGenerator<_NormalGenerator> NormalGenerator;
typedef DefaultNumberGenerator<_LogNormalGenerator> LogNormalGenerator;
 
class Random {
 
private:
    static Random* instance_;
 
    boost::mt19937 rng;
    _RealUniformGenerator uniGen;
 
    std::map<std::string, NumberGenerator*> generators;
 
protected:
    Random(boost::uint32_t seed);
    Random(boost::mt19937 generator);
 
public:
 
    static void initialize(boost::uint32_t seed);
 
    static void initialize(boost::mt19937 generator);
 
    static Random* instance();
    virtual ~Random();
 
    void putGenerator(const std::string& id, NumberGenerator* generator);
 
    NumberGenerator* getGenerator(const std::string& id);
 
    boost::mt19937& engine() {
        return rng;
    }
 
    double nextDouble();
 
    DoubleUniformGenerator createUniDoubleGenerator(double from, double to);
 
    IntUniformGenerator createUniIntGenerator(int from, int to);
 
    TriangleGenerator createTriangleGenerator(double lowerBound, double mostLikely, double upperBound);
 
    CauchyGenerator createCauchyGenerator(double median, double sigma);
 
    ExponentialGenerator createExponentialGenerator(double lambda);
 
    NormalGenerator createNormalGenerator(double mean, double sigma);
 
    LogNormalGenerator createLogNormalGenerator(double mean, double sigma);
};
 
ptrdiff_t uni_random(ptrdiff_t i);
 
template<typename I>
int countOf(I iteratorStart, I iteratorEnd){
  I iterator = iteratorStart;
  int c = 0;
  while(iterator != iteratorEnd){ iterator++; c++; }
  return c;
}
 
template<typename T>
void shuffleList(std::vector<T*>& elementList){
  if(elementList.size() <= 1) return;
  IntUniformGenerator rnd = Random::instance()->createUniIntGenerator(0, elementList.size() - 1);
  T* swap;
  for(size_t i = 0, sz = elementList.size(); i < sz; i++){
    int other = rnd.next();
    swap = elementList[i];
    elementList[i] = elementList[other];
    elementList[other] = swap;
  }
 
  // OR
    // see: http://en.wikipedia.org/wiki/Fisher-Yates_shuffle#The_modern_algorithm
 
  //    DoubleUniformGenerator rnd = Random::instance()->createUniDoubleGenerator(0, 1);
  //    T* swap;
  //    for(int pos = elementList.size() - 1; pos > 0; pos--){
  //        int range = pos + 1;
  //        int other = (int)(rnd.next() * (range));
  //        swap = elementList[pos];
  //        elementList[pos] = elementList[other];
  //        elementList[other] = swap;
  //    }
}
 
template<typename T>
void shuffleSet(std::set<T*>& elementSet, std::vector<T*>& elementList){
  elementList.assign(elementSet.begin(), elementSet.end());
  shuffleList(elementList);
 
  //    // Or- probably faster... ??
  //
  //    elementList.reserve(elementSet.size());
  //    typename std::set<T*>::iterator setIterator = elementSet.begin();
  //    DoubleUniformGenerator rnd = Random::instance()->createUniDoubleGenerator(0,1);
  //    for(int i = 0; i < elementSet.size(); i++){
  //        int j = (int)(rnd.next() * (i + 1));
  //        elementList.push_back(elementList[j]);
  //        elementList[j] = *setIterator;
  //        setIterator++;
  //    }
}
 
template<typename I>
class RandomAccess{
private:
  I it;
  I begin;
 
  int interval;
  int maxLandmark;
  std::vector<std::pair<int, I > > landmarks;
 
public:
 
  RandomAccess(I beginning, int size) :
    interval((int)(sqrt((double)size))), maxLandmark(0), it(beginning), begin(beginning)  {
        landmarks.push_back(std::pair<int, I >(0, beginning));
  }
 
  I get(int index){
    bool place = (index > (maxLandmark + interval));
    typename std::vector<std::pair<int, I > >::iterator lm = landmarks.end();
    while((--lm)->first > index);
    int c = lm->first;
    it    = lm->second;
    if(place){
      while(c != index){
        c++;
        it++;
        if(c % interval == 0){
          landmarks.push_back(std::pair<int, I >(c, it));
          maxLandmark = c;
        }
      }
    }
    else{
      while(c != index){
        c++;
        it++;
      }
    }
    return it;
  }
 
};
 
 
template<typename T, typename I>
void selectNElementsAtRandom(I iterator, int size, unsigned int count, std::set<T*>&selectedElements, bool remove = false){
  // When 'selectedElements' is large in comparison to size, or 'count' is
  // large in comparison to size, or both, a potential problem arises; this is both
  // a performance issue and a potential error (infinite loop). The number
  // of valid selections remaining in the original population will be
  // the size of the population minus the number of members of the original
  // population already in the set of 'selectedElements'. There is a performance
  // problem that can arise if the number of valid selections approaches
  // zero, and an infinite loop arises if that number reaches zero. However, a
  // complication is that we do not assume that all members in
  // 'selectedElements' are drawn from the same population that 'iterator' describes.
 
  // At a general level, the solution to this is that when the number of valid
  // selections is expected to become small, the algorithm will switch to
  // selecting the elements that will _not_ be added to the final set, rather
  // than those that will.
 
  // However, the more immediate issue is determining whether the problem
  // will arise at all, and a key question is how many members of 'selectedElements'
  // are drawn from the population. This, however, is costly to check, so
  // it should only be checked if it seems possible to lead to a problem.
 
  // A ceiling; the actual value may be lower
  int maxAlreadySelected = selectedElements.size();
 
  // We are not certain _any_ elements in the selectedElements set are from this population
  int knownSelectedElementsFromPop = 0;
 
  // This is a floor; there are at least this many available, but possibly more
  int minAvailable = size - maxAlreadySelected;
 
  // If you are requesting more than might be available, you must determine
  // the actual number available.
  if(count > minAvailable){
    if(maxAlreadySelected > 0){
      I tempIt = iterator;
      for(int i = 0; i < size; i++){
        T* ptr = (&**tempIt);
        typename std::set<T*>::iterator found = selectedElements.find(ptr);
        if(found != selectedElements.end()) knownSelectedElementsFromPop++;
        tempIt++;
      }
      maxAlreadySelected = knownSelectedElementsFromPop;
      minAvailable = size - maxAlreadySelected;
    }
  }
  // If removing the original elements, will need a copy of them
  typename std::vector<T*> tempToRemove;
  if(remove){
    tempToRemove.assign(selectedElements.begin(), selectedElements.end());
  }
 
  // There is no way to satisfy the request; copy all and return
  if(count > minAvailable){
    I it = iterator;
    for(int i = 0; i < size; i++){
      selectedElements.insert(&**it);
      it++;
    }
  }
  else{
    RandomAccess<I> ra(iterator, size);
    IntUniformGenerator rnd = Random::instance()->createUniIntGenerator(0, size - 1);
    T* ptr;
 
    // If the count of elements is very high in proportion to the population, faster to
    // choose agents that will _not_ be in the final set and then switch...
    // First the normal case, in which the number of agents is low
    if(count <= ((double)minAvailable)*0.6){ // Tests suggest that 2/3 is probably too high, so adjusted to .6; TODO optimize or analytically solve
      for (unsigned int i = 0; i < count; i++)
                do {  ptr = &**ra.get(rnd.next()); } while(!(selectedElements.insert(ptr).second));
    }
    else{ // Now the other case
      std::set<T*> elementsThatWillNotBeAdded;
 
      if(selectedElements.size() > 0){
        if(selectedElements.size() == knownSelectedElementsFromPop){ // Maybe we already checked and they ALL belong
          elementsThatWillNotBeAdded.insert(selectedElements.begin(), selectedElements.end());
        }
        else{
          I tempIt = iterator;
          for(int i = 0; i < size; i++){
            ptr = (&**tempIt);
            if(selectedElements.find(ptr) != selectedElements.end()) elementsThatWillNotBeAdded.insert(ptr);
            tempIt++;
          }
        }
      }
      // Note: duplicate element will not be inserted; failure will leave size of elementsThatWillNotBeAdded unchanged.
      while((size - elementsThatWillNotBeAdded.size()) > count) elementsThatWillNotBeAdded.insert(&**ra.get(rnd.next()));
 
      // Once the set of elements that will not be added is complete, cycle through the iterator and add
      // all of those elements that are not in elementsThatWillNotBeAdded
      I toAdd = iterator;
      typename std::set<T*>::iterator notFound = elementsThatWillNotBeAdded.end();
      for(int i = 0; i < size; i++){
        ptr = &**toAdd;
        if(elementsThatWillNotBeAdded.find(ptr) == notFound) selectedElements.insert(ptr);
        toAdd++;
      }
    }
  }
  // Before returning, remove the elements from the set, if the user requested
  if(remove){
    typename std::vector<T*>::iterator toRemove = tempToRemove.begin();
    while(toRemove != tempToRemove.end()){
      selectedElements.erase(*toRemove);
      toRemove++;
    }
  }
}
 
 
template<typename T, typename I>
void selectNElementsAtRandom(I iteratorStart, I iteratorEnd, int count, std::set<T*>&selectedElements, bool remove = false){
  selectNElementsAtRandom(iteratorStart, countOf(iteratorStart, iteratorEnd), count, selectedElements, remove);
}
 
template<typename T, typename I>
void selectNElementsInRandomOrder(I iterator, int size, int count, std::vector<T*>& selectedElements, bool remove = false){
  // Transfer all elements from the vector to a set
  std::set<T*> selectedElementSet;
  selectedElementSet.insert(selectedElements.begin(), selectedElements.end());
  selectedElements.clear();
  selectNElementsAtRandom(iterator, size, count, selectedElementSet, remove);
  shuffleSet(selectedElementSet, selectedElements);
}
 
template<typename T, typename I>
void selectNElementsInRandomOrder(I iteratorStart, I iteratorEnd, int count, std::vector<T*>& selectedElements, bool remove = false){
  selectNElementsInRandomOrder(iteratorStart, countOf(iteratorStart, iteratorEnd), count, selectedElements, remove);
}
 
 
}
 
#endif /* RANDOM_H_ */