Mercurial > octave
changeset 28106:c08d7b53f1a2
new function, rng
* rng.m: New function.
* scripts/general/module.mk: Update.
author | Guillaume Flandin |
---|---|
date | Wed, 19 Feb 2020 10:02:09 -0500 |
parents | ff3b8b21a890 |
children | d320728d5d06 |
files | scripts/general/module.mk scripts/general/rng.m |
diffstat | 2 files changed, 270 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/general/module.mk Tue Feb 18 21:57:14 2020 -0800 +++ b/scripts/general/module.mk Wed Feb 19 10:02:09 2020 -0500 @@ -60,6 +60,7 @@ %reldir%/repelem.m \ %reldir%/repmat.m \ %reldir%/rescale.m \ + %reldir%/rng.m \ %reldir%/rot90.m \ %reldir%/rotdim.m \ %reldir%/shift.m \
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/general/rng.m Wed Feb 19 10:02:09 2020 -0500 @@ -0,0 +1,269 @@ +######################################################################## +## +## Copyright (C) 2020 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/>. +## +######################################################################## + +## -*- texinfo -*- +## @deftypefn {} {} rng (@var{seed}) +## @deftypefnx {} {} rng (@var{seed}, @var{generator}) +## @deftypefnx {} {} rng ("shuffle") +## @deftypefnx {} {} rng ("shuffle", @var{generator}) +## @deftypefnx {} {} rng ("default") +## @deftypefnx {} {@var{s} =} rng () +## @deftypefnx {} {} rng (@var{s}) +## @deftypefnx {} {@var{s} =} rng (@dots{}) +## Set or query the seed of the random number generator used by rand and randn. +## +## Input scalar numeric value @qcode{seed} is used to initialize the state +## vector of the random number generator. +## +## The optional string @var{generator} specifies the type of random number +## generator to be used. Its value can be @qcode{"twister"}, @qcode{"v5uniform"} +## or @qcode{"v5normal"}. The @qcode{"twister"} keyword is described below. +## @qcode{"v5uniform"} and @qcode{"v5normal"} refer to older versions of Octave +## that used to use a different random number generator. +## +## The state or seed of the random number generator can be reset to a new random +## value using the @qcode{"shuffle"} keyword. +## +## The state or seed of the random number generator can be reset to its default +## values using the @qcode{"default"} keyword. The default values are to use the +## Mersenne Twister with a seed of 0. +## +## The optional return value @var{s} contains the state of the random number +## generator at the time the function is called (i.e. before it might be +## modified according to the input arguments). It is encoded as a structure +## variable with three fields: @qcode{"Type"}, @qcode{"Seed"} and +## @qcode{"State"}. +## Later, the random number generator can be restored to the state @var{s} using +## rng (s). +## +## By default, and with the @qcode{"twister"} option, pseudo-random sequences +## are computed using the Mersenne Twister with a period of @math{2^{19937}-1} +## (See @nospell{M. Matsumoto and T. Nishimura}, +## @cite{Mersenne Twister: A 623-dimensionally equidistributed uniform +## pseudorandom number generator}, +## @nospell{ACM} Trans.@: on Modeling and Computer Simulation Vol.@: 8, No.@: 1, +## pp.@: 3--30, January 1998, +## @url{http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html}). +## Do @strong{not} use for cryptography without securely hashing several +## returned values together, otherwise the generator state can be learned after +## reading 624 consecutive values. +## +## @seealso{rand, randn} +## @end deftypefn + +function retval = rng (varargin) + + if (nargin > 2 || nargout > 1) + print_usage (); + endif + + ## Store current settings of random number generator + ## FIXME: there doesn't seem to be a way to query the type of generator + ## currently used in Octave - assume "twister". + ## FIXME: there doesn't seem to be a way to query the seed initialization + ## value - use "Not applicable". + ## FIXME: rand and randn use different generators - storing both states. + ## For older Matlab generators (v4, v5), the settings are stored like this: + ## struct ("Type","Legacy", "Seed", "Not applicable", "State",{[],[],...}) + s = struct (... + "Type", "twister",... # generator name + "Seed", "Not applicable",... # seed initialization value + "State", {{rand("state"), randn("state")}}); # internal state of the generator + + if (! nargin) + retval = s; + return; + endif + + if (isscalar (varargin{1}) && isnumeric (varargin{1}) && isreal (varargin{1}) && varargin{1} >= 0) + s_rand = s_randn = varargin{1}; + generator = check_generator (varargin(2:end)); + + elseif (ischar (varargin{1}) && strcmpi (varargin{1}, "shuffle")) + ## Seed the random number generator based on the current time + s_rand = s_randn = "reset"; # or sum (1000*clock) + generator = check_generator (varargin(2:end)); + + elseif (ischar (varargin{1}) && strcmpi (varargin{1}, "default") && nargin == 1) + generator = "twister"; + s_rand = s_randn = 0; # In Matlab, seed 0 corresponds to 5489 + + elseif (isstruct (varargin{1}) && isscalar (varargin{1}) && nargin == 1) + if (numfields (varargin{1}) != 3 || ! isfield (varargin{1}, "Type") + || ! isfield (varargin{1}, "Seed") || ! isfield (varargin{1}, "State")) + error ("Input structure not compatible with the one returned by rng ()"); + endif + ## Only the internal state "State" and generator type "Type" are needed + generator = varargin{1}.Type; + if (iscell (varargin{1}.State)) + [s_rand, s_randn] = deal (varargin{1}.State{:}); + else + s_rand = s_randn = varargin{1}.State; + endif + + else + print_usage (); + endif + + ## Set the type of random number generator and seed it + if (isempty (generator)) + generator = s.Type; + endif + switch generator + case "twister" + rand ("state", s_rand); + randn ("state", s_randn); + + case "legacy" + rand ("seed", s_rand); + randn ("seed", s_randn); + + case "v5uniform" + rand ("seed", s_rand); + + case "v5normal" + randn ("seed", s_randn); + + otherwise + error ("Unknown type of random number generator"); + + endswitch + + if (nargout > 0) + retval = s; + endif + +endfunction + + +function gen = check_generator (val) + if (isempty (val)) + gen = ""; + return; + elseif (! iscellstr (val)) + error ("Second input must be a type of random number generator"); + endif + gen = tolower (char (val)); + if (ismember (gen, {"simdtwister", "combrecursive", "philox", "threefry", "multfibonacci", "v4"})) + error ("This random number generator is not available in Octave"); + elseif (! ismember (gen, {"twister", "v5uniform", "v5normal"})) + error ("This type of random number generator is unknown"); + endif +endfunction + + +%!test +%! rng (42); +%! ru1 = rand (); +%! rn1 = randn (); +%! rng (42); +%! ru2 = rand (); +%! rn2 = randn (); +%! assert (ru2, ru1); +%! assert (rn2, rn1); +%! s1 = rng (); +%! s2 = rng (42); +%! assert (isequal (s1, s2)); +%! ru1 = rand (); +%! rn1 = randn (); +%! s3 = rng (42); +%! ru2 = rand (); +%! rn2 = randn (); +%! assert (ru2, ru1); +%! assert (rn2, rn1); + +%!test +%! rng (42, "twister"); +%! ru1 = rand (); +%! rn1 = randn (); +%! rng (42, "twister"); +%! ru2 = rand (); +%! rn2 = randn (); +%! assert (ru2, ru1); +%! assert (rn2, rn1); +%! s1 = rng (); +%! s2 = rng (42, "twister"); +%! assert (isequal (s1, s2)); +%! ru1 = rand (); +%! rn1 = randn (); +%! s3 = rng (42, "twister"); +%! ru2 = rand (); +%! rn2 = randn (); +%! assert (ru2, ru1); +%! assert (rn2, rn1); + +%!test +%! rng ("shuffle"); +%! rng ("shuffle", "twister"); +%! s = rng ("shuffle"); +%! assert (! isequal (s, rng ("shuffle"))); +%! s = rng ("shuffle", "twister"); +%! assert (! isequal (s, rng ("shuffle", "twister"))); + +%!test +%! rng ("default"); +%! ru1 = rand (); +%! rn1 = randn (); +%! rng (0, "twister"); +%! ru2 = rand (); +%! rn2 = randn (); +%! assert (ru2, ru1); +%! assert (rn2, rn1); +%! rng (0, "twister"); +%! s = rng ("default"); +%! assert (isequal (s, rng ())); + +%!test +%! s = rng (); +%! ru1 = rand (); +%! rn1 = randn (); +%! rng (s); +%! ru2 = rand (); +%! rn2 = randn (); +%! assert (ru2, ru1); +%! assert (rn2, rn1); +%! rng (42); rand (1,2); x = rand (1,2); +%! rng (42); rand (1,2); s = rng (); y = rand (1,2); +%! assert (x, y); +%! rng (s); z = rand (1,2); +%! assert (x, z); +%! s1 = rng (); +%! s2 = rng (s1); +%! assert (isequal (s1, s2)); + +## Test input validation +%!error <Invalid call> rng (1, 2, 3) +%!error <Invalid call> [a, b] = rng () +%!error <Invalid call> rng ({}) +%!error <Invalid call> rng ("unknown") +%!error <Invalid call> rng (eye (2)) +%!error <Invalid call> rng (struct ("Seed", {1 ,2})) +%!error <Invalid call> rng (-1) +%!error <Invalid call> rng (struct ("Type",[],"State",[],"Seed",[]), 2) +%!error <Invalid call> rng ("default", "twister") +%!error <Second input must be a type of random number generator> rng (0, struct ()) +%!error <This type of random number generator is unknown> rng (0, "unknown") +%!error <Second input must be a type of random number generator> rng ("shuffle", struct ()) +%!error <This type of random number generator is unknown> rng ("shuffle", "unknown")