random_generators.h

Go to the documentation of this file.

/*
 *  random_generators.h
 *
 *  This file is part of NEST.
 *
 *  Copyright (C) 2004 The NEST Initiative
 *
 *  NEST 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 2 of the License, or
 *  (at your option) any later version.
 *
 *  NEST 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 NEST.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
 
#ifndef RANDOM_GENERATORS_H
#define RANDOM_GENERATORS_H
 
// C++ includes:
#include <algorithm>
#include <initializer_list>
#include <iterator>
#include <memory>
#include <random>
#include <type_traits>
#include <utility>
 
// libnestutil includes:
#include "randutils.hpp"
 
// nestkernel includes:
#include "node_collection.h"
 
namespace nest
{
 
// Forward declarations
class BaseRandomGenerator;
template < typename DistributionT >
class RandomDistribution;
 
// Type definitions
using RngPtr = BaseRandomGenerator*;
 
using uniform_int_distribution = RandomDistribution< std::uniform_int_distribution< unsigned long > >;
using uniform_real_distribution = RandomDistribution< std::uniform_real_distribution<> >;
using normal_distribution = RandomDistribution< std::normal_distribution<> >;
using lognormal_distribution = RandomDistribution< std::lognormal_distribution<> >;
using binomial_distribution = RandomDistribution< std::binomial_distribution< unsigned long > >;
using gamma_distribution = RandomDistribution< std::gamma_distribution<> >;
using exponential_distribution = RandomDistribution< std::exponential_distribution<> >;
using discrete_distribution = RandomDistribution< std::discrete_distribution< unsigned long > >;
using poisson_distribution = RandomDistribution< std::poisson_distribution< unsigned long > >;
 
 
class BaseRandomGenerator
{
public:
  virtual ~BaseRandomGenerator()
  {
  }
 
  virtual unsigned long operator()( std::uniform_int_distribution< unsigned long >& d ) = 0;
  virtual double operator()( std::uniform_real_distribution<>& d ) = 0;
  virtual double operator()( std::normal_distribution<>& d ) = 0;
  virtual double operator()( std::lognormal_distribution<>& d ) = 0;
  virtual unsigned long operator()( std::binomial_distribution< unsigned long >& d ) = 0;
  virtual double operator()( std::gamma_distribution<>& d ) = 0;
  virtual double operator()( std::exponential_distribution<>& d ) = 0;
  virtual unsigned long operator()( std::discrete_distribution< unsigned long >& d ) = 0;
  virtual unsigned long operator()( std::poisson_distribution< unsigned long >& d ) = 0;
 
  virtual unsigned long operator()( std::uniform_int_distribution< unsigned long >& d,
    std::uniform_int_distribution< unsigned long >::param_type& p ) = 0;
  virtual double operator()( std::uniform_real_distribution<>& d, std::uniform_real_distribution<>::param_type& p ) = 0;
  virtual double operator()( std::normal_distribution<>& d, std::normal_distribution<>::param_type& p ) = 0;
  virtual double operator()( std::lognormal_distribution<>& d, std::lognormal_distribution<>::param_type& p ) = 0;
  virtual unsigned long operator()( std::binomial_distribution< unsigned long >& d,
    std::binomial_distribution< unsigned long >::param_type& p ) = 0;
  virtual double operator()( std::gamma_distribution<>& d, std::gamma_distribution<>::param_type& p ) = 0;
  virtual double operator()( std::exponential_distribution<>& d, std::exponential_distribution<>::param_type& p ) = 0;
  virtual unsigned long operator()( std::discrete_distribution< unsigned long >& d,
    std::discrete_distribution< unsigned long >::param_type& p ) = 0;
  virtual unsigned long operator()( std::poisson_distribution< unsigned long >& d,
    std::poisson_distribution< unsigned long >::param_type& p ) = 0;
 
  virtual double drand() = 0;
 
  virtual unsigned long ulrand( unsigned long N ) = 0;
 
  virtual void stable_sample( NodeCollection::const_iterator first,
    NodeCollection::const_iterator last,
    std::back_insert_iterator< std::vector< NodeIDTriple > > dest,
    size_t n ) = 0;
 
  virtual void partial_shuffle( std::vector< size_t >& data, const size_t n ) = 0;
};
 
template < typename RandomEngineT >
class RandomGenerator final : public BaseRandomGenerator
{
public:
  using result_type = typename RandomEngineT::result_type;
 
  RandomGenerator() = delete;
  RandomGenerator( const RandomEngineT& rng ) = delete;
 
  explicit RandomGenerator( std::initializer_list< std::uint32_t > seed )
    : rng_()
    , uniform_double_dist_0_1_( 0.0, 1.0 )
  {
    // Melissa O'Neill's seed sequence generator which provides better distributed seeds
    // than std::seed_seq; see https://www.pcg-random.org/posts/developing-a-seed_seq-alternative.html
    randutils::seed_seq_fe128 sseq( seed );
    rng_.seed( sseq );
  }
 
  inline unsigned long
  operator()( std::uniform_int_distribution< unsigned long >& d ) override
  {
    return d( rng_ );
  }
 
  inline double
  operator()( std::uniform_real_distribution<>& d ) override
  {
    return d( rng_ );
  }
 
  inline double
  operator()( std::normal_distribution<>& d ) override
  {
    return d( rng_ );
  }
 
  inline double
  operator()( std::lognormal_distribution<>& d ) override
  {
    return d( rng_ );
  }
 
  inline unsigned long
  operator()( std::binomial_distribution< unsigned long >& d ) override
  {
    return d( rng_ );
  }
 
  inline double
  operator()( std::gamma_distribution<>& d ) override
  {
    return d( rng_ );
  }
 
  inline double
  operator()( std::exponential_distribution<>& d ) override
  {
    return d( rng_ );
  }
 
  inline unsigned long
  operator()( std::discrete_distribution< unsigned long >& d ) override
  {
    return d( rng_ );
  }
 
  inline unsigned long
  operator()( std::poisson_distribution< unsigned long >& d ) override
  {
    return d( rng_ );
  }
 
  inline unsigned long
  operator()( std::uniform_int_distribution< unsigned long >& d,
    std::uniform_int_distribution< unsigned long >::param_type& p ) override
  {
    return d( rng_, p );
  }
 
  inline double
  operator()( std::uniform_real_distribution<>& d, std::uniform_real_distribution<>::param_type& p ) override
  {
    return d( rng_, p );
  }
 
  inline double
  operator()( std::normal_distribution<>& d, std::normal_distribution<>::param_type& p ) override
  {
    return d( rng_, p );
  }
 
  inline double
  operator()( std::lognormal_distribution<>& d, std::lognormal_distribution<>::param_type& p ) override
  {
    return d( rng_, p );
  }
 
  inline unsigned long
  operator()( std::binomial_distribution< unsigned long >& d,
    std::binomial_distribution< unsigned long >::param_type& p ) override
  {
    return d( rng_, p );
  }
 
  inline double
  operator()( std::gamma_distribution<>& d, std::gamma_distribution<>::param_type& p ) override
  {
    return d( rng_, p );
  }
 
  inline double
  operator()( std::exponential_distribution<>& d, std::exponential_distribution<>::param_type& p ) override
  {
    return d( rng_, p );
  }
 
  inline unsigned long
  operator()( std::discrete_distribution< unsigned long >& d,
    std::discrete_distribution< unsigned long >::param_type& p ) override
  {
    return d( rng_, p );
  }
 
  inline unsigned long
  operator()( std::poisson_distribution< unsigned long >& d,
    std::poisson_distribution< unsigned long >::param_type& p ) override
  {
    return d( rng_, p );
  }
 
  inline double
  drand() override
  {
    return uniform_double_dist_0_1_( rng_ );
  }
 
  inline unsigned long
  ulrand( unsigned long N ) override
  {
    std::uniform_int_distribution< unsigned long >::param_type param( 0, N - 1 );
    return uniform_ulong_dist_( rng_, param );
  }
 
  inline void
  stable_sample( NodeCollection::const_iterator first,
    NodeCollection::const_iterator last,
    std::back_insert_iterator< std::vector< NodeIDTriple > > dest,
    size_t n ) override
  {
    std::sample( first, last, dest, n, rng_ );
  }
 
  inline void
  partial_shuffle( std::vector< size_t >& data, const size_t n ) override
  {
    const size_t num_elems = std::min( n, data.size() );
 
    // Partial Fisher-Yates-Durstenfeld-Knuth: at each step select an element from the remaining range
    // [ next, data.size() ) and swap it to the front. Draws num_elems random numbers and
    // leaves the selection in random order.
    for ( size_t next = 0, remaining = data.size(); next < num_elems; ++next, --remaining )
    {
      const size_t rnd = next + ulrand( remaining );
      std::swap( data[ next ], data[ rnd ] );
    }
    data.resize( num_elems );  // keep only the selected elements.
  }
 
private:
  RandomEngineT rng_;  
  std::uniform_int_distribution< unsigned long > uniform_ulong_dist_;
  std::uniform_real_distribution<> uniform_double_dist_0_1_;
};
 
class BaseRandomGeneratorFactory
{
public:
  virtual ~BaseRandomGeneratorFactory() {};
 
  virtual RngPtr create( std::initializer_list< std::uint32_t > seed_initializer ) const = 0;
};
 
template < typename RandomEngineT >
class RandomGeneratorFactory final : public BaseRandomGeneratorFactory
{
public:
  inline RngPtr
  create( std::initializer_list< std::uint32_t > seed_initializer ) const override
  {
    return new RandomGenerator< RandomEngineT >( seed_initializer );
  }
};
 
template < typename DistributionT >
class RandomDistribution
{
public:
  using result_type = typename DistributionT::result_type;
  using param_type = typename DistributionT::param_type;
 
  RandomDistribution() = default;
  ~RandomDistribution() = default;
 
  inline result_type
  operator()( RngPtr g )
  {
    static_assert( std::is_same< result_type, unsigned long >::value or std::is_same< result_type, double >::value,
      "result_type of the distribution must be unsigned long or double" );
    return g->operator()( distribution_ );
  }
 
  inline result_type
  operator()( RngPtr g, param_type& params )
  {
    static_assert( std::is_same< result_type, unsigned long >::value or std::is_same< result_type, double >::value,
      "result_type of the distribution must be unsigned long or double" );
    return g->operator()( distribution_, params );
  }
 
  inline void
  param( const param_type& params )
  {
    distribution_.param( params );
  }
 
  inline result_type
  min() const
  {
    return distribution_.min();
  }
 
  inline result_type
  max() const
  {
    return distribution_.max();
  }
 
private:
  DistributionT distribution_;  
};
 
}  // namespace nest
 
#endif /* #ifndef RANDOM_GENERATORS_H */