Mercurial > hg > octave-terminal
changeset 7537:a2950622f070
make octave_rand a proper singleton class
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Wed, 27 Feb 2008 04:33:39 -0500 |
parents | 4dda6fbc8ba6 |
children | 2c4b0cbda85a |
files | liboctave/ChangeLog liboctave/oct-rand.cc liboctave/oct-rand.h |
diffstat | 3 files changed, 380 insertions(+), 263 deletions(-) [+] |
line wrap: on
line diff
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,7 @@ +2008-02-27 John W. Eaton <jwe@octave.org> + + * oct-rand.cc (class octave_rand): Make it a proper singleton class. + 2008-02-26 John W. Eaton <jwe@octave.org> * oct-rand.cc (get_dist_id): Fix typo.
--- a/liboctave/oct-rand.cc +++ b/liboctave/oct-rand.cc @@ -39,23 +39,6 @@ #include "randgamma.h" #include "mach-info.h" -static const int unknown_dist = 0; -static const int uniform_dist = 1; -static const int normal_dist = 2; -static const int expon_dist = 3; -static const int poisson_dist = 4; -static const int gamma_dist = 5; - -// Current distribution of random numbers. -static int current_distribution = uniform_dist; - -// Has the seed/state been set yet? -static bool old_initialized = false; -static bool new_initialized = false; -static bool use_old_generators = false; - -std::map<int, ColumnVector> rand_states; - extern "C" { F77_RET_T @@ -86,160 +69,39 @@ F77_FUNC (setcgn, SETCGN) (const octave_idx_type&); } -static octave_idx_type -force_to_fit_range (octave_idx_type i, octave_idx_type lo, octave_idx_type hi) -{ - assert (hi > lo && lo >= 0 && hi > lo); - - i = i > 0 ? i : -i; - - if (i < lo) - i = lo; - else if (i > hi) - i = i % hi; - - return i; -} - -// Make the random number generator give us a different sequence every -// time we start octave unless we specifically set the seed. The -// technique used below will cycle monthly, but it it does seem to -// work ok to give fairly different seeds each time Octave starts. - -static void -do_old_initialization (void) -{ - octave_localtime tm; - int stored_distribution = current_distribution; - F77_FUNC (setcgn, SETCGN) (uniform_dist); +octave_rand *octave_rand::instance = 0; - int hour = tm.hour() + 1; - int minute = tm.min() + 1; - int second = tm.sec() + 1; - - octave_idx_type s0 = tm.mday() * hour * minute * second; - octave_idx_type s1 = hour * minute * second; - - s0 = force_to_fit_range (s0, 1, 2147483563); - s1 = force_to_fit_range (s1, 1, 2147483399); - - F77_FUNC (setall, SETALL) (s0, s1); - F77_FUNC (setcgn, SETCGN) (stored_distribution); +octave_rand::octave_rand (void) + : current_distribution (uniform_dist), use_old_generators (false), + rand_states () +{ + initialize_ranlib_generators (); - old_initialized = true; -} - -static ColumnVector -get_internal_state (void) -{ - ColumnVector s (MT_N + 1); - - OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1); - - oct_get_state (tmp); - - for (octave_idx_type i = 0; i <= MT_N; i++) - s.elem (i) = static_cast<double> (tmp [i]); - - return s; + initialize_mersenne_twister (); } -static inline void -save_state (void) -{ - rand_states[current_distribution] = get_internal_state ();; -} - -static void -initialize_rand_states (void) -{ - if (! new_initialized) - { - oct_init_by_entropy (); - - ColumnVector s = get_internal_state (); - - rand_states[uniform_dist] = s; - rand_states[normal_dist] = s; - rand_states[expon_dist] = s; - rand_states[poisson_dist] = s; - rand_states[gamma_dist] = s; - - new_initialized = true; - } -} - -static inline void -maybe_initialize (void) +bool +octave_rand::instance_ok (void) { - if (use_old_generators) - { - if (! old_initialized) - do_old_initialization (); - } - else - { - if (! new_initialized) - initialize_rand_states (); - } -} + bool retval = true; + + if (! instance) + instance = new octave_rand (); -static int -get_dist_id (const std::string& d) -{ - int retval = unknown_dist; + if (! instance) + { + (*current_liboctave_error_handler) + ("unable to create octave_rand object!"); - if (d == "uniform" || d == "rand") - retval = uniform_dist; - else if (d == "normal" || d == "randn") - retval = normal_dist; - else if (d == "exponential" || d == "rande") - retval = expon_dist; - else if (d == "poisson" || d == "randp") - retval = poisson_dist; - else if (d == "gamma" || d == "randg") - retval = gamma_dist; - else - (*current_liboctave_error_handler) - ("rand: invalid distribution `%s'", d.c_str ()); + retval = false; + } return retval; } -static void -set_internal_state (const ColumnVector& s) +double +octave_rand::do_seed (void) { - octave_idx_type len = s.length (); - octave_idx_type n = len < MT_N + 1 ? len : MT_N + 1; - - OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1); - - for (octave_idx_type i = 0; i < n; i++) - tmp[i] = static_cast<uint32_t> (s.elem(i)); - - if (len == MT_N + 1 && tmp[MT_N] <= MT_N && tmp[MT_N] > 0) - oct_set_state (tmp); - else - oct_init_by_array (tmp, len); -} - -static inline void -switch_to_generator (int dist) -{ - if (dist != current_distribution) - { - current_distribution = dist; - - set_internal_state (rand_states[dist]); - } -} - -double -octave_rand::seed (void) -{ - if (! old_initialized) - do_old_initialization (); - union d2i { double d; octave_idx_type i[2]; }; union d2i u; @@ -258,13 +120,26 @@ return u.d; } +static octave_idx_type +force_to_fit_range (octave_idx_type i, octave_idx_type lo, octave_idx_type hi) +{ + assert (hi > lo && lo >= 0 && hi > lo); + + i = i > 0 ? i : -i; + + if (i < lo) + i = lo; + else if (i > hi) + i = i % hi; + + return i; +} + void -octave_rand::seed (double s) +octave_rand::do_seed (double s) { use_old_generators = true; - maybe_initialize (); - int i0, i1; union d2i { double d; octave_idx_type i[2]; }; union d2i u; @@ -288,21 +163,16 @@ } ColumnVector -octave_rand::state (const std::string& d) +octave_rand::do_state (const std::string& d) { - if (! new_initialized) - initialize_rand_states (); - return rand_states[d.empty () ? current_distribution : get_dist_id (d)]; } void -octave_rand::state (const ColumnVector& s, const std::string& d) +octave_rand::do_state (const ColumnVector& s, const std::string& d) { use_old_generators = false; - maybe_initialize (); - int old_dist = current_distribution; int new_dist = d.empty () ? current_distribution : get_dist_id (d); @@ -321,12 +191,10 @@ } std::string -octave_rand::distribution (void) +octave_rand::do_distribution (void) { std::string retval; - maybe_initialize (); - switch (current_distribution) { case uniform_dist: @@ -359,7 +227,7 @@ } void -octave_rand::distribution (const std::string& d) +octave_rand::do_distribution (const std::string& d) { int id = get_dist_id (d); @@ -393,50 +261,40 @@ } void -octave_rand::uniform_distribution (void) +octave_rand::do_uniform_distribution (void) { - maybe_initialize (); - switch_to_generator (uniform_dist); F77_FUNC (setcgn, SETCGN) (uniform_dist); } void -octave_rand::normal_distribution (void) +octave_rand::do_normal_distribution (void) { - maybe_initialize (); - switch_to_generator (normal_dist); F77_FUNC (setcgn, SETCGN) (normal_dist); } void -octave_rand::exponential_distribution (void) +octave_rand::do_exponential_distribution (void) { - maybe_initialize (); - switch_to_generator (expon_dist); F77_FUNC (setcgn, SETCGN) (expon_dist); } void -octave_rand::poisson_distribution (void) +octave_rand::do_poisson_distribution (void) { - maybe_initialize (); - switch_to_generator (poisson_dist); F77_FUNC (setcgn, SETCGN) (poisson_dist); } void -octave_rand::gamma_distribution (void) +octave_rand::do_gamma_distribution (void) { - maybe_initialize (); - switch_to_generator (gamma_dist); F77_FUNC (setcgn, SETCGN) (gamma_dist); @@ -444,10 +302,8 @@ double -octave_rand::scalar (double a) +octave_rand::do_scalar (double a) { - maybe_initialize (); - double retval = 0.0; if (use_old_generators) @@ -526,6 +382,167 @@ return retval; } +Matrix +octave_rand::do_matrix (octave_idx_type n, octave_idx_type m, double a) +{ + Matrix retval; + + if (n >= 0 && m >= 0) + { + retval.resize (n, m); + + if (n > 0 && m > 0) + fill (retval.capacity(), retval.fortran_vec(), a); + } + else + (*current_liboctave_error_handler) ("rand: invalid negative argument"); + + return retval; +} + +NDArray +octave_rand::do_nd_array (const dim_vector& dims, double a) +{ + NDArray retval; + + if (! dims.all_zero ()) + { + retval.resize (dims); + + fill (retval.capacity(), retval.fortran_vec(), a); + } + + return retval; +} + +Array<double> +octave_rand::do_vector (octave_idx_type n, double a) +{ + Array<double> retval; + + if (n > 0) + { + retval.resize (n); + + fill (retval.capacity (), retval.fortran_vec (), a); + } + else if (n < 0) + (*current_liboctave_error_handler) ("rand: invalid negative argument"); + + return retval; +} + +// Make the random number generator give us a different sequence every +// time we start octave unless we specifically set the seed. The +// technique used below will cycle monthly, but it it does seem to +// work ok to give fairly different seeds each time Octave starts. + +void +octave_rand::initialize_ranlib_generators (void) +{ + octave_localtime tm; + int stored_distribution = current_distribution; + F77_FUNC (setcgn, SETCGN) (uniform_dist); + + int hour = tm.hour() + 1; + int minute = tm.min() + 1; + int second = tm.sec() + 1; + + octave_idx_type s0 = tm.mday() * hour * minute * second; + octave_idx_type s1 = hour * minute * second; + + s0 = force_to_fit_range (s0, 1, 2147483563); + s1 = force_to_fit_range (s1, 1, 2147483399); + + F77_FUNC (setall, SETALL) (s0, s1); + F77_FUNC (setcgn, SETCGN) (stored_distribution); +} + +void +octave_rand::initialize_mersenne_twister (void) +{ + oct_init_by_entropy (); + + ColumnVector s = get_internal_state (); + + rand_states[uniform_dist] = s; + rand_states[normal_dist] = s; + rand_states[expon_dist] = s; + rand_states[poisson_dist] = s; + rand_states[gamma_dist] = s; +} + +ColumnVector +octave_rand::get_internal_state (void) +{ + ColumnVector s (MT_N + 1); + + OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1); + + oct_get_state (tmp); + + for (octave_idx_type i = 0; i <= MT_N; i++) + s.elem (i) = static_cast<double> (tmp [i]); + + return s; +} + +void +octave_rand::save_state (void) +{ + rand_states[current_distribution] = get_internal_state ();; +} + +int +octave_rand::get_dist_id (const std::string& d) +{ + int retval = unknown_dist; + + if (d == "uniform" || d == "rand") + retval = uniform_dist; + else if (d == "normal" || d == "randn") + retval = normal_dist; + else if (d == "exponential" || d == "rande") + retval = expon_dist; + else if (d == "poisson" || d == "randp") + retval = poisson_dist; + else if (d == "gamma" || d == "randg") + retval = gamma_dist; + else + (*current_liboctave_error_handler) + ("rand: invalid distribution `%s'", d.c_str ()); + + return retval; +} + +void +octave_rand::set_internal_state (const ColumnVector& s) +{ + octave_idx_type len = s.length (); + octave_idx_type n = len < MT_N + 1 ? len : MT_N + 1; + + OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1); + + for (octave_idx_type i = 0; i < n; i++) + tmp[i] = static_cast<uint32_t> (s.elem(i)); + + if (len == MT_N + 1 && tmp[MT_N] <= MT_N && tmp[MT_N] > 0) + oct_set_state (tmp); + else + oct_init_by_array (tmp, len); +} + +void +octave_rand::switch_to_generator (int dist) +{ + if (dist != current_distribution) + { + current_distribution = dist; + + set_internal_state (rand_states[dist]); + } +} + #define MAKE_RAND(len) \ do \ { \ @@ -539,11 +556,9 @@ } \ while (0) -static void -fill_rand (octave_idx_type len, double *v, double a) +void +octave_rand::fill (octave_idx_type len, double *v, double a) { - maybe_initialize (); - if (len < 1) return; @@ -630,58 +645,6 @@ return; } -Matrix -octave_rand::matrix (octave_idx_type n, octave_idx_type m, double a) -{ - Matrix retval; - - if (n >= 0 && m >= 0) - { - retval.resize (n, m); - - if (n > 0 && m > 0) - fill_rand (retval.capacity(), retval.fortran_vec(), a); - } - else - (*current_liboctave_error_handler) ("rand: invalid negative argument"); - - return retval; -} - -NDArray -octave_rand::nd_array (const dim_vector& dims, double a) -{ - NDArray retval; - - if (! dims.all_zero ()) - { - retval.resize (dims); - - fill_rand (retval.capacity(), retval.fortran_vec(), a); - } - - return retval; -} - -Array<double> -octave_rand::vector (octave_idx_type n, double a) -{ - maybe_initialize (); - - Array<double> retval; - - if (n > 0) - { - retval.resize (n); - - fill_rand (retval.capacity(), retval.fortran_vec(), a); - } - else if (n < 0) - (*current_liboctave_error_handler) ("rand: invalid negative argument"); - - return retval; -} - /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/oct-rand.h +++ b/liboctave/oct-rand.h @@ -23,59 +23,209 @@ #if !defined (octave_rand_h) #define octave_rand_h 1 +#include <map> #include <string> #include "dColVector.h" #include "dMatrix.h" #include "dNDArray.h" +#include "lo-ieee.h" -struct +class OCTAVE_API octave_rand { +protected: + + octave_rand (void); + +public: + + ~octave_rand (void) { } + + static bool instance_ok (void); + // Return the current seed. - static double seed (void); + static double seed (void) + { + return instance_ok () ? instance->do_seed () : octave_NaN; + } // Set the seed. - static void seed (double s); + static void seed (double s) + { + if (instance_ok ()) + instance->do_seed (s); + } // Return the current state. - static ColumnVector state (const std::string& d = std::string ()); + static ColumnVector state (const std::string& d = std::string ()) + { + return instance_ok () ? instance->do_state (d) : ColumnVector (); + } // Set the current state/ static void state (const ColumnVector &s, - const std::string& d = std::string ()); + const std::string& d = std::string ()) + { + if (instance_ok ()) + instance->do_state (s, d); + } // Return the current distribution. - static std::string distribution (void); + static std::string distribution (void) + { + return instance_ok () ? instance->do_distribution () : std::string (); + } // Set the current distribution. May be either "uniform" (the // default), "normal", "exponential", "poisson", or "gamma". - static void distribution (const std::string& d); + static void distribution (const std::string& d) + { + if (instance_ok ()) + instance->do_distribution (d); + } - static void uniform_distribution (void); + static void uniform_distribution (void) + { + if (instance_ok ()) + instance->do_uniform_distribution (); + } - static void normal_distribution (void); + static void normal_distribution (void) + { + if (instance_ok ()) + instance->do_normal_distribution (); + } - static void exponential_distribution (void); + static void exponential_distribution (void) + { + if (instance_ok ()) + instance->do_exponential_distribution (); + } - static void poisson_distribution (void); + static void poisson_distribution (void) + { + if (instance_ok ()) + instance->do_poisson_distribution (); + } - static void gamma_distribution (void); + static void gamma_distribution (void) + { + if (instance_ok ()) + instance->do_gamma_distribution (); + } // Return the next number from the sequence. - static double scalar (double a = 1.); + static double scalar (double a = 1.0) + { + return instance_ok () ? instance->do_scalar (a) : octave_NaN; + } // Return a matrix of numbers from the sequence, filled in column // major order. - static Matrix matrix (octave_idx_type r, octave_idx_type c, double a = 1.); + static Matrix matrix (octave_idx_type r, octave_idx_type c, double a = 1.0) + { + return instance_ok () ? instance->do_matrix (r, c, a) : Matrix (); + } // Return an N-dimensional array of numbers from the sequence, // filled in column major order. - static NDArray nd_array (const dim_vector& dims, double a = 1.); + static NDArray nd_array (const dim_vector& dims, double a = 1.0) + { + return instance_ok () ? instance->do_nd_array (dims, a) : NDArray (); + } // Return an array of numbers from the sequence. - static Array<double> vector (octave_idx_type n, double a = 1.); + static Array<double> vector (octave_idx_type n, double a = 1.0) + { + return instance_ok () ? instance->do_vector (n, a) : Array<double> (); + } + +private: + + static octave_rand *instance; + + enum + { + unknown_dist, + uniform_dist, + normal_dist, + expon_dist, + poisson_dist, + gamma_dist + }; + + // Current distribution of random numbers. + int current_distribution; + + // If TRUE, use old RANLIB generators. Otherwise, use Mersenne + // Twister generator. + bool use_old_generators; + + // Saved MT states. + std::map<int, ColumnVector> rand_states; + + // Return the current seed. + double do_seed (void); + + // Set the seed. + void do_seed (double s); + + // Return the current state. + ColumnVector do_state (const std::string& d); + + // Set the current state/ + void do_state (const ColumnVector &s, const std::string& d); + + // Return the current distribution. + std::string do_distribution (void); + + // Set the current distribution. May be either "uniform" (the + // default), "normal", "exponential", "poisson", or "gamma". + void do_distribution (const std::string& d); + + void do_uniform_distribution (void); + + void do_normal_distribution (void); + + void do_exponential_distribution (void); + + void do_poisson_distribution (void); + + void do_gamma_distribution (void); + + // Return the next number from the sequence. + double do_scalar (double a = 1.); + + // Return a matrix of numbers from the sequence, filled in column + // major order. + Matrix do_matrix (octave_idx_type r, octave_idx_type c, double a = 1.); + + // Return an N-dimensional array of numbers from the sequence, + // filled in column major order. + NDArray do_nd_array (const dim_vector& dims, double a = 1.); + + // Return an array of numbers from the sequence. + Array<double> do_vector (octave_idx_type n, double a = 1.); + + // Some helper functions. + + void initialize_ranlib_generators (void); + + void initialize_mersenne_twister (void); + + ColumnVector get_internal_state (void); + + void save_state (void); + + int get_dist_id (const std::string& d); + + void set_internal_state (const ColumnVector& s); + + void switch_to_generator (int dist); + + void fill (octave_idx_type len, double *v, double a); }; #endif