1 /*
2  *  random_manager.cpp
3  *
4  *  This file is part of NEST.
5  *
6  *  Copyright (C) 2004 The NEST Initiative
7  *
8  *  NEST is free software: you can redistribute it and/or modify
9  *  it under the terms of the GNU General Public License as published by
10  *  the Free Software Foundation, either version 2 of the License, or
11  *  (at your option) any later version.
12  *
13  *  NEST is distributed in the hope that it will be useful,
14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  *  GNU General Public License for more details.
17  *
18  *  You should have received a copy of the GNU General Public License
19  *  along with NEST.  If not, see <http://www.gnu.org/licenses/>.
20  *
21  */
22 
23 #include "random_manager.h"
24 
25 // C++ includes:
26 #include <random>
27 #include <utility>
28 
29 // Includes from nestkernel:
30 #include "kernel_manager.h"
31 #include "random_generators.h"
32 #include "vp_manager_impl.h"
33 
34 // Includes from libnestutil:
35 #include "Random123/conventional/Engine.hpp"
36 #include "Random123/philox.h"
37 #include "Random123/threefry.h"
38 
39 
40 const std::string nest::RandomManager::DEFAULT_RNG_TYPE_ = "mt19937_64";
41 
42 const std::uint32_t nest::RandomManager::DEFAULT_BASE_SEED_ = 143202461;
43 
44 const std::uint32_t nest::RandomManager::RANK_SYNCED_SEEDER_ = 0xc229212d;
45 const std::uint32_t nest::RandomManager::THREAD_SYNCED_SEEDER_ = 0x37722d5e;
46 const std::uint32_t nest::RandomManager::THREAD_SPECIFIC_SEEDER_ = 0xb84c9bae;
47 
48 
RandomManager()49 nest::RandomManager::RandomManager()
50   : current_rng_type_( DEFAULT_RNG_TYPE_ )
51   , base_seed_( DEFAULT_BASE_SEED_ )
52   , rank_synced_rng_( nullptr )
53 {
54 }
55 
~RandomManager()56 nest::RandomManager::~RandomManager()
57 {
58   finalize();
59 }
60 
61 void
initialize()62 nest::RandomManager::initialize()
63 {
64   register_rng_type< std::mt19937 >( "mt19937" );
65   register_rng_type< std::mt19937_64 >( "mt19937_64" );
66 #ifdef HAVE_RANDOM123
67   register_rng_type< r123::Engine< r123::Philox4x32 > >( "Philox_32" );
68 #if R123_USE_PHILOX_64BIT
69   register_rng_type< r123::Engine< r123::Philox4x64 > >( "Philox_64" );
70 #endif
71   register_rng_type< r123::Engine< r123::Threefry4x32 > >( "Threefry_32" );
72 #if R123_USE_64BIT
73   register_rng_type< r123::Engine< r123::Threefry4x64 > >( "Threefry_64" );
74 #endif
75 #endif
76 
77   current_rng_type_ = DEFAULT_RNG_TYPE_;
78   base_seed_ = DEFAULT_BASE_SEED_;
79 
80   reset_rngs_();
81 }
82 
83 void
finalize()84 nest::RandomManager::finalize()
85 {
86   for ( auto& it : rng_types_ )
87   {
88     delete it.second;
89   }
90 
91   rng_types_.clear();
92   vp_specific_rngs_.clear();
93 }
94 
95 void
reset_rngs_()96 nest::RandomManager::reset_rngs_()
97 {
98   // Delete existing RNGs.
99   delete rank_synced_rng_;
100 
101   auto delete_rngs = []( std::vector< RngPtr >& rng_vec )
102   {
103     for ( auto rng : rng_vec )
104     {
105       delete rng;
106     }
107   };
108 
109   delete_rngs( vp_synced_rngs_ );
110   delete_rngs( vp_specific_rngs_ );
111 
112   // Create new RNGs of the currently used RNG type.
113   rank_synced_rng_ = rng_types_[ current_rng_type_ ]->create( { base_seed_, RANK_SYNCED_SEEDER_ } );
114 
115   vp_synced_rngs_.resize( kernel().vp_manager.get_num_threads() );
116   vp_specific_rngs_.resize( kernel().vp_manager.get_num_threads() );
117 
118 #pragma omp parallel
119   {
120     const auto tid = kernel().vp_manager.get_thread_id();
121     const std::uint32_t vp = kernel().vp_manager.get_vp();
122     vp_synced_rngs_[ tid ] = rng_types_[ current_rng_type_ ]->create( { base_seed_, THREAD_SYNCED_SEEDER_ } );
123     vp_specific_rngs_[ tid ] = rng_types_[ current_rng_type_ ]->create( { base_seed_, THREAD_SPECIFIC_SEEDER_, vp } );
124   }
125 }
126 
127 void
get_status(DictionaryDatum & d)128 nest::RandomManager::get_status( DictionaryDatum& d )
129 {
130   ArrayDatum rng_types;
131   for ( auto rng = rng_types_.begin(); rng != rng_types_.end(); ++rng )
132   {
133     rng_types.push_back( rng->first );
134   }
135 
136   def< ArrayDatum >( d, names::rng_types, rng_types );
137   def< long >( d, names::rng_seed, base_seed_ );
138   def< std::string >( d, names::rng_type, current_rng_type_ );
139 }
140 
141 void
set_status(const DictionaryDatum & d)142 nest::RandomManager::set_status( const DictionaryDatum& d )
143 {
144   long rng_seed;
145   bool rng_seed_updated = updateValue< long >( d, names::rng_seed, rng_seed );
146 
147   if ( rng_seed_updated )
148   {
149     if ( not( 0 < rng_seed and rng_seed < ( 1L << 32 ) ) )
150     {
151       throw BadProperty( "RNG seed must be in (0, 2^32-1)." );
152     }
153 
154     base_seed_ = static_cast< std::uint32_t >( rng_seed );
155   }
156 
157   std::string rng_type;
158   bool rng_type_updated = updateValue< std::string >( d, names::rng_type, rng_type );
159 
160   if ( rng_type_updated )
161   {
162     auto rng = rng_types_.find( rng_type );
163     if ( rng == rng_types_.end() )
164     {
165       std::string msg = "'%1' is not a known RNG type. See /rng_types for available types";
166       throw BadProperty( String::compose( msg, rng_type ) );
167     }
168 
169     current_rng_type_ = rng_type;
170   }
171 
172   // If number of threads has been changed, we need to update the RNGs.
173   bool n_threads_updated = d->known( names::local_num_threads );
174   if ( n_threads_updated or rng_seed_updated or rng_type_updated )
175   {
176     reset_rngs_();
177   }
178 }
179 
180 void
check_rng_synchrony() const181 nest::RandomManager::check_rng_synchrony() const
182 {
183   // Compare more than a single number to avoid false negatives
184   const long NUM_ROUNDS = 5;
185 
186   // We check rank-synchrony even if we are on a single process to keep the code simple.
187   for ( auto n = 0; n < NUM_ROUNDS; ++n )
188   {
189     const auto r = rank_synced_rng_->drand();
190     const auto min = kernel().mpi_manager.min_cross_ranks( r );
191     const auto max = kernel().mpi_manager.max_cross_ranks( r );
192     if ( min != max )
193     {
194       throw KernelException( "Rank-synchronized random number generators are out of sync." );
195     }
196   }
197 
198   // We check thread-synchrony under all circumstances to keep the code simple.
199   for ( auto n = 0; n < NUM_ROUNDS; ++n )
200   {
201     const index num_threads = kernel().vp_manager.get_num_threads();
202     double local_min = std::numeric_limits< double >::max();
203     double local_max = std::numeric_limits< double >::min();
204     for ( index t = 0; t < num_threads; ++t )
205     {
206       const auto r = vp_synced_rngs_[ t ]->drand();
207       local_min = std::min( r, local_min );
208       local_max = std::max( r, local_max );
209     }
210 
211     // Finding the local min and max on each thread and then determining the
212     // global min/max, ensures that all ranks will learn about sync errors.
213     const long min = kernel().mpi_manager.min_cross_ranks( local_min );
214     const long max = kernel().mpi_manager.max_cross_ranks( local_max );
215     if ( min != max )
216     {
217       throw KernelException( "Thread-synchronized random number generators are out of sync." );
218     }
219   }
220 }
221 
222 template < typename RNG_TYPE >
223 void
register_rng_type(std::string name)224 nest::RandomManager::register_rng_type( std::string name )
225 {
226   rng_types_.insert( std::make_pair( name, new RandomGeneratorFactory< RNG_TYPE >() ) );
227 }
228