oct-rand.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2003-2022 The Octave Project Developers
//
// See the file COPYRIGHT.md in the top-level directory of this
// distribution or <https://octave.org/copyright/>.
//
// This file is part of Octave.
//
// Octave 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.
//
// Octave 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 Octave; see the file COPYING.  If not, see
// <https://www.gnu.org/licenses/>.
//
////////////////////////////////////////////////////////////////////////
 
#if defined (HAVE_CONFIG_H)
#  include "config.h"
#endif
 
#include <cassert>
#include <cstdint>
 
#include <limits>
 
#include "lo-error.h"
#include "lo-ieee.h"
#include "lo-mappers.h"
#include "lo-ranlib-proto.h"
#include "mach-info.h"
#include "oct-locbuf.h"
#include "oct-rand.h"
#include "oct-time.h"
#include "quit.h"
#include "randgamma.h"
#include "randmtzig.h"
#include "randpoisson.h"
#include "singleton-cleanup.h"
 
namespace octave
{
  rand *rand::m_instance = nullptr;
 
  rand::rand (void)
    : m_current_distribution (uniform_dist), m_use_old_generators (false),
      m_rand_states ()
  {
    initialize_ranlib_generators ();
 
    initialize_mersenne_twister ();
  }
 
  bool rand::instance_ok (void)
  {
    bool retval = true;
 
    if (! m_instance)
      {
        m_instance = new rand ();
        singleton_cleanup_list::add (cleanup_instance);
      }
 
    return retval;
  }
 
  double rand::do_seed (void)
  {
    union d2i { double d; int32_t i[2]; };
    union d2i u;
 
    mach_info::float_format ff = mach_info::native_float_format ();
 
    switch (ff)
      {
      case mach_info::flt_fmt_ieee_big_endian:
        F77_FUNC (getsd, GETSD) (u.i[1], u.i[0]);
        break;
 
      default:
        F77_FUNC (getsd, GETSD) (u.i[0], u.i[1]);
        break;
      }
 
    return u.d;
  }
 
  static int32_t
  force_to_fit_range (int32_t i, int32_t lo, int32_t hi)
  {
    assert (hi > lo && lo >= 0);
 
    i = (i > 0 ? i : -i);
 
    if (i < lo)
      i = lo;
    else if (i > hi)
      i = i % hi;
 
    return i;
  }
 
  void rand::do_seed (double s)
  {
    m_use_old_generators = true;
 
    int i0, i1;
    union d2i { double d; int32_t i[2]; };
    union d2i u;
    u.d = s;
 
    mach_info::float_format ff = mach_info::native_float_format ();
 
    switch (ff)
      {
      case mach_info::flt_fmt_ieee_big_endian:
        i1 = force_to_fit_range (u.i[0], 1, 2147483563);
        i0 = force_to_fit_range (u.i[1], 1, 2147483399);
        break;
 
      default:
        i0 = force_to_fit_range (u.i[0], 1, 2147483563);
        i1 = force_to_fit_range (u.i[1], 1, 2147483399);
        break;
      }
 
    F77_FUNC (setsd, SETSD) (i0, i1);
  }
 
  void rand::do_reset (void)
  {
    m_use_old_generators = true;
    initialize_ranlib_generators ();
  }
 
  uint32NDArray rand::do_state (const std::string& d)
  {
    return m_rand_states[d.empty () ? m_current_distribution : get_dist_id (d)];
  }
 
  void rand::do_state (const uint32NDArray& s, const std::string& d)
  {
    m_use_old_generators = false;
 
    int old_dist = m_current_distribution;
 
    int new_dist = (d.empty () ? m_current_distribution : get_dist_id (d));
 
    uint32NDArray saved_state;
 
    if (old_dist != new_dist)
      saved_state = get_internal_state ();
 
    set_internal_state (s);
 
    m_rand_states[new_dist] = get_internal_state ();
 
    if (old_dist != new_dist)
      m_rand_states[old_dist] = saved_state;
  }
 
  void rand::do_reset (const std::string& d)
  {
    m_use_old_generators = false;
 
    int old_dist = m_current_distribution;
 
    int new_dist = (d.empty () ? m_current_distribution : get_dist_id (d));
 
    uint32NDArray saved_state;
 
    if (old_dist != new_dist)
      saved_state = get_internal_state ();
 
    init_mersenne_twister ();
    m_rand_states[new_dist] = get_internal_state ();
 
    if (old_dist != new_dist)
      m_rand_states[old_dist] = saved_state;
  }
 
  std::string rand::do_distribution (void)
  {
    std::string retval;
 
    switch (m_current_distribution)
      {
      case uniform_dist:
        retval = "uniform";
        break;
 
      case normal_dist:
        retval = "normal";
        break;
 
      case expon_dist:
        retval = "exponential";
        break;
 
      case poisson_dist:
        retval = "poisson";
        break;
 
      case gamma_dist:
        retval = "gamma";
        break;
 
      default:
        (*current_liboctave_error_handler)
          ("rand: invalid distribution ID = %d", m_current_distribution);
        break;
      }
 
    return retval;
  }
 
  void rand::do_distribution (const std::string& d)
  {
    int id = get_dist_id (d);
 
    switch (id)
      {
      case uniform_dist:
        rand::uniform_distribution ();
        break;
 
      case normal_dist:
        rand::normal_distribution ();
        break;
 
      case expon_dist:
        rand::exponential_distribution ();
        break;
 
      case poisson_dist:
        rand::poisson_distribution ();
        break;
 
      case gamma_dist:
        rand::gamma_distribution ();
        break;
 
      default:
        (*current_liboctave_error_handler)
          ("rand: invalid distribution ID = %d", id);
        break;
      }
  }
 
  void rand::do_uniform_distribution (void)
  {
    switch_to_generator (uniform_dist);
 
    F77_FUNC (setcgn, SETCGN) (uniform_dist);
  }
 
  void rand::do_normal_distribution (void)
  {
    switch_to_generator (normal_dist);
 
    F77_FUNC (setcgn, SETCGN) (normal_dist);
  }
 
  void rand::do_exponential_distribution (void)
  {
    switch_to_generator (expon_dist);
 
    F77_FUNC (setcgn, SETCGN) (expon_dist);
  }
 
  void rand::do_poisson_distribution (void)
  {
    switch_to_generator (poisson_dist);
 
    F77_FUNC (setcgn, SETCGN) (poisson_dist);
  }
 
  void rand::do_gamma_distribution (void)
  {
    switch_to_generator (gamma_dist);
 
    F77_FUNC (setcgn, SETCGN) (gamma_dist);
  }
 
  template <>
  OCTAVE_API double rand::uniform<double> (void)
  {
    double retval;
 
    if (m_use_old_generators)
      F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval);
    else
      retval = rand_uniform<double> ();
 
    return retval;
  }
 
  template <>
  OCTAVE_API double rand::normal<double> (void)
  {
    double retval;
 
    if (m_use_old_generators)
      F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval);
    else
      retval = rand_normal<double> ();
 
    return retval;
  }
 
  template <>
  OCTAVE_API double rand::exponential<double> (void)
  {
    double retval;
 
    if (m_use_old_generators)
      F77_FUNC (dgenexp, DGENEXP) (1.0, retval);
    else
      retval = rand_exponential<double> ();
 
    return retval;
  }
 
  template <>
  OCTAVE_API double rand::poisson<double> (double a)
  {
    double retval;
 
    if (m_use_old_generators)
      {
        if (a < 0.0 || ! math::isfinite (a))
          retval = numeric_limits<double>::NaN ();
        else
          {
            // workaround bug in ignpoi, by calling with different Mu
            F77_FUNC (dignpoi, DIGNPOI) (a + 1, retval);
            F77_FUNC (dignpoi, DIGNPOI) (a, retval);
          }
      }
    else
      retval = rand_poisson<double> (a);
 
    return retval;
  }
 
  template <>
  OCTAVE_API double rand::gamma<double> (double a)
  {
    double retval;
 
    if (m_use_old_generators)
      {
        if (a <= 0.0 || ! math::isfinite (a))
          retval = numeric_limits<double>::NaN ();
        else
          F77_FUNC (dgengam, DGENGAM) (1.0, a, retval);
      }
    else
      retval = rand_gamma<double> (a);
 
    return retval;
  }
 
  template <>
  OCTAVE_API float rand::uniform<float> (void)
  {
    float retval;
 
    if (m_use_old_generators)
      F77_FUNC (fgenunf, FGENUNF) (0.0f, 1.0f, retval);
    else
      retval = rand_uniform<float> ();
 
    return retval;
  }
 
  template <>
  OCTAVE_API float rand::normal<float> (void)
  {
    float retval;
 
    if (m_use_old_generators)
      F77_FUNC (fgennor, FGENNOR) (0.0f, 1.0f, retval);
    else
      retval = rand_normal<float> ();
 
    return retval;
  }
 
  template <>
  OCTAVE_API float rand::exponential<float> (void)
  {
    float retval;
 
    if (m_use_old_generators)
      F77_FUNC (fgenexp, FGENEXP) (1.0f, retval);
    else
      retval = rand_exponential<float> ();
 
    return retval;
  }
 
  template <>
  OCTAVE_API float rand::poisson<float> (float a)
  {
    float retval;
 
    if (m_use_old_generators)
      {
        if (a < 0.0f || ! math::isfinite (a))
          retval = numeric_limits<float>::NaN ();
        else
          {
            // workaround bug in ignpoi, by calling with different Mu
            F77_FUNC (fignpoi, FIGNPOI) (a + 1, retval);
            F77_FUNC (fignpoi, FIGNPOI) (a, retval);
          }
      }
    else
      {
        // Keep poisson distribution in double precision for accuracy
        retval = rand_poisson<double> (a);
      }
 
    return retval;
  }
 
  template <>
  OCTAVE_API float rand::gamma<float> (float a)
  {
    float retval;
 
    if (m_use_old_generators)
      {
        if (a <= 0.0f || ! math::isfinite (a))
          retval = numeric_limits<float>::NaN ();
        else
          F77_FUNC (fgengam, FGENGAM) (1.0f, a, retval);
      }
    else
      retval = rand_gamma<float> (a);
 
    return retval;
  }
 
  template <typename T>
  T rand::do_scalar (T a)
  {
    T retval = 0;
 
    switch (m_current_distribution)
      {
      case uniform_dist:
        retval = uniform<T> ();
        break;
 
      case normal_dist:
        retval = normal<T> ();
        break;
 
      case expon_dist:
        retval = exponential<T> ();
        break;
 
      case poisson_dist:
        retval = poisson<T> (a);
        break;
 
      case gamma_dist:
        retval = gamma<T> (a);
        break;
 
      default:
        (*current_liboctave_error_handler)
          ("rand: invalid distribution ID = %d", m_current_distribution);
        break;
      }
 
    if (! m_use_old_generators)
      save_state ();
 
    return retval;
  }
 
  template OCTAVE_API double rand::do_scalar<double> (double);
  template OCTAVE_API float rand::do_scalar<float> (float);
 
  template <typename T>
  Array<T>
  rand::do_vector (octave_idx_type n, T a)
  {
    Array<T> retval;
 
    if (n > 0)
      {
        retval.clear (n, 1);
 
        fill (retval.numel (), retval.fortran_vec (), a);
      }
    else if (n < 0)
      (*current_liboctave_error_handler) ("rand: invalid negative argument");
 
    return retval;
  }
 
  template OCTAVE_API Array<double>
  rand::do_vector<double> (octave_idx_type, double);
  template OCTAVE_API Array<float>
  rand::do_vector<float> (octave_idx_type, float);
 
  NDArray rand::do_nd_array (const dim_vector& dims, double a)
  {
    NDArray retval;
 
    if (! dims.all_zero ())
      {
        retval.clear (dims);
 
        fill (retval.numel (), retval.fortran_vec (), a);
      }
 
    return retval;
  }
 
  FloatNDArray rand::do_float_nd_array (const dim_vector& dims, float a)
  {
    FloatNDArray retval;
 
    if (! dims.all_zero ())
      {
        retval.clear (dims);
 
        fill (retval.numel (), retval.fortran_vec (), a);
      }
 
    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 does seem to
  // work ok to give fairly different seeds each time Octave starts.
 
  void rand::initialize_ranlib_generators (void)
  {
    sys::localtime tm;
    int stored_distribution = m_current_distribution;
    F77_FUNC (setcgn, SETCGN) (uniform_dist);
 
    int hour = tm.hour () + 1;
    int minute = tm.min () + 1;
    int second = tm.sec () + 1;
 
    int32_t s0 = tm.mday () * hour * minute * second;
    int32_t 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 rand::initialize_mersenne_twister (void)
  {
    uint32NDArray s;
 
    init_mersenne_twister ();
    s = get_internal_state ();
    m_rand_states[uniform_dist] = s;
 
    init_mersenne_twister ();
    s = get_internal_state ();
    m_rand_states[normal_dist] = s;
 
    init_mersenne_twister ();
    s = get_internal_state ();
    m_rand_states[expon_dist] = s;
 
    init_mersenne_twister ();
    s = get_internal_state ();
    m_rand_states[poisson_dist] = s;
 
    init_mersenne_twister ();
    s = get_internal_state ();
    m_rand_states[gamma_dist] = s;
 
    // All of the initializations above have messed with the internal state.
    // Restore the state of the currently selected distribution.
    set_internal_state (m_rand_states[m_current_distribution]);
  }
 
  uint32NDArray rand::get_internal_state (void)
  {
    uint32NDArray s (dim_vector (MT_N + 1, 1));
 
    get_mersenne_twister_state (reinterpret_cast<uint32_t *> (s.fortran_vec ()));
 
    return s;
  }
 
  void rand::save_state (void)
  {
    m_rand_states[m_current_distribution] = get_internal_state ();;
  }
 
  int 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 rand::set_internal_state (const uint32NDArray& s)
  {
    octave_idx_type len = s.numel ();
 
    const uint32_t *sdata = reinterpret_cast <const uint32_t *> (s.data ());
 
    if (len == MT_N + 1 && sdata[MT_N] <= MT_N && sdata[MT_N] > 0)
      set_mersenne_twister_state (sdata);
    else
      init_mersenne_twister (sdata, len);
  }
 
  void rand::switch_to_generator (int dist)
  {
    if (dist != m_current_distribution)
      {
        m_current_distribution = dist;
 
        set_internal_state (m_rand_states[dist]);
      }
  }
 
  void rand::fill (octave_idx_type len, double *v, double a)
  {
    if (len < 1)
      return;
 
    switch (m_current_distribution)
      {
      case uniform_dist:
        if (m_use_old_generators)
          std::generate_n (v, len, [](void) { double x; F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, x); return x; });
        else
          rand_uniform<double> (len, v);
        break;
 
      case normal_dist:
        if (m_use_old_generators)
          std::generate_n (v, len, [](void) { double x; F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, x); return x; });
        else
          rand_normal<double> (len, v);
        break;
 
      case expon_dist:
        if (m_use_old_generators)
          std::generate_n (v, len, [](void) { double x; F77_FUNC (dgenexp, DGENEXP) (1.0, x); return x; });
        else
          rand_exponential<double> (len, v);
        break;
 
      case poisson_dist:
        if (m_use_old_generators)
          {
            if (a < 0.0 || ! math::isfinite (a))
              std::fill_n (v, len, numeric_limits<double>::NaN ());
            else
              {
                // workaround bug in ignpoi, by calling with different Mu
                double tmp;
                F77_FUNC (dignpoi, DIGNPOI) (a + 1, tmp);
                std::generate_n (v, len, [a](void) { double x; F77_FUNC (dignpoi, DIGNPOI) (a, x); return x; });
              }
          }
        else
          rand_poisson<double> (a, len, v);
        break;
 
      case gamma_dist:
        if (m_use_old_generators)
          {
            if (a <= 0.0 || ! math::isfinite (a))
              std::fill_n (v, len, numeric_limits<double>::NaN ());
            else
              std::generate_n (v, len, [a](void) { double x; F77_FUNC (dgengam, DGENGAM) (1.0, a, x); return x; });
          }
        else
          rand_gamma<double> (a, len, v);
        break;
 
      default:
        (*current_liboctave_error_handler)
          ("rand: invalid distribution ID = %d", m_current_distribution);
        break;
      }
 
    save_state ();
 
    return;
  }
 
  void rand::fill (octave_idx_type len, float *v, float a)
  {
    if (len < 1)
      return;
 
    switch (m_current_distribution)
      {
      case uniform_dist:
        if (m_use_old_generators)
          std::generate_n (v, len, [](void) { float x; F77_FUNC (fgenunf, FGENUNF) (0.0f, 1.0f, x); return x; });
        else
          rand_uniform<float> (len, v);
        break;
 
      case normal_dist:
        if (m_use_old_generators)
          std::generate_n (v, len, [](void) { float x; F77_FUNC (fgennor, FGENNOR) (0.0f, 1.0f, x); return x; });
        else
          rand_normal<float> (len, v);
        break;
 
      case expon_dist:
        if (m_use_old_generators)
          std::generate_n (v, len, [](void) { float x; F77_FUNC (fgenexp, FGENEXP) (1.0f, x); return x; });
        else
          rand_exponential<float> (len, v);
        break;
 
      case poisson_dist:
        if (m_use_old_generators)
          {
            if (a < 0.0f || ! math::isfinite (a))
              std::fill_n (v, len, numeric_limits<float>::NaN ());
            else
              {
                // workaround bug in ignpoi, by calling with different Mu
                float tmp;
                F77_FUNC (fignpoi, FIGNPOI) (a + 1, tmp);
                std::generate_n (v, len, [a](void) { float x; F77_FUNC (fignpoi, FIGNPOI) (a, x); return x; });
              }
          }
        else
          rand_poisson<float> (a, len, v);
        break;
 
      case gamma_dist:
        if (m_use_old_generators)
          {
            if (a <= 0.0f || ! math::isfinite (a))
              std::fill_n (v, len, numeric_limits<float>::NaN ());
            else
              std::generate_n (v, len, [a](void) { float x; F77_FUNC (fgengam, FGENGAM) (1.0f, a, x); return x; });
          }
        else
          rand_gamma<float> (a, len, v);
        break;
 
      default:
        (*current_liboctave_error_handler)
          ("rand: invalid distribution ID = %d", m_current_distribution);
        break;
      }
 
    save_state ();
 
    return;
  }
}