Re: cellfun and parcellfun

[David Bateman]

Robert I A Patterson wrote:
> >     
> I understand the reasons, although it leads to quite awkward behaviour
> from a Monte Carlo point of view. Imagine I run 8 sample paths using 8
> processes via parcellfunc.  If I have already used some random numbers,
> all 8 paths will be identical, but if this happens to be the first thing
> I do when I start Octave I will get 8 different samples (provided each
> process gets different data from /dev/urandom).  I wonder if some kind
> of note in the parcellfun documentation might be appropriate:
> "Note that each process will inherit a copy of the random number
> generator state that existed before the call to parcellfun.  If that
> state is uninitialised, each process will perform the default
> initialisation of its random number generator separately.  In the latter
> case it is possible, but not certain that each process will initialise
> its random number generator differently."
>   
This is why I suggested that a feature like rand("state", "reset") that 
reinitialized the state of the generator would be a good feature to 
have. Using this as the first call in each forked process would ensure 
that the state is different always in all forked processes. It should be 
relatively easy to do this except for the care needed with the fact that 
the randlib generators all use the same state and the mersenne twister 
based generators have separate states and so rand("state", "reset") 
should change only the uniform generators state, but rand("seed", 
"reset") should change all of them.

Something like the attached should do it. Anyone mind if I commit it?

D.

diff --git a/liboctave/oct-rand.cc b/liboctave/oct-rand.cc
--- a/liboctave/oct-rand.cc
+++ b/liboctave/oct-rand.cc
@@ -165,6 +165,13 @@
   F77_FUNC (setsd, SETSD) (i0, i1);
 }
 
+void
+octave_rand::do_reset (void)
+{
+  use_old_generators = true;
+  initialize_ranlib_generators ();
+}
+
 ColumnVector
 octave_rand::do_state (const std::string& d)
 {
@@ -193,6 +200,27 @@
     rand_states[old_dist] = saved_state;
 }
 
+void
+octave_rand::do_reset (const std::string& d)
+{
+  use_old_generators = false;
+
+  int old_dist = current_distribution;
+
+  int new_dist = d.empty () ? current_distribution : get_dist_id (d);
+
+  ColumnVector saved_state;
+
+  if (old_dist != new_dist)
+    saved_state = get_internal_state ();
+
+  oct_init_by_entropy ();
+  rand_states[new_dist] = get_internal_state ();
+
+  if (old_dist != new_dist)
+    rand_states[old_dist] = saved_state;
+}
+
 std::string
 octave_rand::do_distribution (void)
 {
diff --git a/liboctave/oct-rand.h b/liboctave/oct-rand.h
--- a/liboctave/oct-rand.h
+++ b/liboctave/oct-rand.h
@@ -58,6 +58,13 @@
       instance->do_seed (s);
   }
 
+  // Reset the seed.
+  static void reset (void)
+  {
+    if (instance_ok ())
+      instance->do_reset ();
+  }
+
   // Return the current state.
   static ColumnVector state (const std::string& d = std::string ())
   {
@@ -72,6 +79,13 @@
       instance->do_state (s, d);
   }
   
+  // Reset the current state/
+  static void reset (const std::string& d)
+  {
+    if (instance_ok ())
+      instance->do_reset (d);
+  }
+  
   // Return the current distribution.
   static std::string distribution (void)
   {
@@ -172,12 +186,18 @@
   // Set the seed.
   void do_seed (double s);
 
+  // Reset the seed.
+  void do_reset ();
+
   // 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);
   
+  // Reset the current state/
+  void do_reset (const std::string& d);
+  
   // Return the current distribution.
   std::string do_distribution (void);
 
diff --git a/src/DLD-FUNCTIONS/rand.cc b/src/DLD-FUNCTIONS/rand.cc
--- a/src/DLD-FUNCTIONS/rand.cc
+++ b/src/DLD-FUNCTIONS/rand.cc
@@ -243,16 +243,25 @@
                     if (! error_state)
                       octave_rand::seed (d);
                   }
+                else if (args(idx+1).is_string () 
+                         && args(idx+1).string_value() == "reset")
+                  octave_rand::reset ();
                 else
                   error ("%s: seed must be a real scalar", fcn);
               }
             else if (ts == "state" || ts == "twister")
               {
-                ColumnVector s = 
-                  ColumnVector (args(idx+1).vector_value(false, true));
+                if (args(idx+1).is_string ()
+                    && args(idx+1).string_value() == "reset")
+                  octave_rand::reset (fcn);
+                else
+                  {
+                    ColumnVector s = 
+                      ColumnVector (args(idx+1).vector_value(false, true));
 
-                if (! error_state)
-                  octave_rand::state (s, fcn);
+                    if (! error_state)
+                      octave_rand::state (s, fcn);
+                  }
               }
             else
               error ("%s: unrecognized string argument", fcn);
@@ -314,7 +323,9 @@
 @deftypefn  {Loadable Function} {} rand (@var{x})\n\
 @deftypefnx {Loadable Function} {} rand (@var{n}, @var{m})\n\
 @deftypefnx {Loadable Function} {} rand (\"state\", @var{x})\n\
address@hidden {Loadable Function} {} rand (\"state\", \"reset\")\n\
 @deftypefnx {Loadable Function} {} rand (\"seed\", @var{x})\n\
address@hidden {Loadable Function} {} rand (\"seed\", \"reset\")\n\
 Return a matrix with random elements uniformly distributed on the\n\
 interval (0, 1).  The arguments are handled the same as the arguments\n\
 for @code{eye}.\n\
@@ -376,6 +387,9 @@
 @code{rand} to use the old generators, only setting the seed will.\n\
 To cause @code{rand} to once again use the new generators, the\n\
 keyword \"state\" should be used to reset the state of the @code{rand}.\n\
+\n\
+The state or seed of the generator can be reset to a new random value\n\
+using the \"reset\" keyword.\n\
 @seealso{randn, rande, randg, randp}\n\
 @end deftypefn")
 {
@@ -474,7 +488,9 @@
 @deftypefn  {Loadable Function} {} randn (@var{x})\n\
 @deftypefnx {Loadable Function} {} randn (@var{n}, @var{m})\n\
 @deftypefnx {Loadable Function} {} randn (\"state\", @var{x})\n\
address@hidden {Loadable Function} {} randn (\"state\", \"reset\")\n\
 @deftypefnx {Loadable Function} {} randn (\"seed\", @var{x})\n\
address@hidden {Loadable Function} {} randn (\"seed\", \"reset\")\n\
 Return a matrix with normally distributed random\n\
 elements having zero mean and variance one.  The arguments are\n\
 handled the same as the arguments for @code{rand}.\n\
@@ -547,7 +563,9 @@
 @deftypefn  {Loadable Function} {} rande (@var{x})\n\
 @deftypefnx {Loadable Function} {} rande (@var{n}, @var{m})\n\
 @deftypefnx {Loadable Function} {} rande (\"state\", @var{x})\n\
address@hidden {Loadable Function} {} rande (\"state\", \"reset\")\n\
 @deftypefnx {Loadable Function} {} rande (\"seed\", @var{x})\n\
address@hidden {Loadable Function} {} rande (\"seed\", \"reset\")\n\
 Return a matrix with exponentially distributed random elements.  The\n\
 arguments are handled the same as the arguments for @code{rand}.\n\
 \n\
@@ -620,7 +638,9 @@
 @deftypefn  {Loadable Function} {} randg (@var{a}, @var{x})\n\
 @deftypefnx {Loadable Function} {} randg (@var{a}, @var{n}, @var{m})\n\
 @deftypefnx {Loadable Function} {} randg (\"state\", @var{x})\n\
address@hidden {Loadable Function} {} randg (\"state\", \"reset\")\n\
 @deftypefnx {Loadable Function} {} randg (\"seed\", @var{x})\n\
address@hidden {Loadable Function} {} randg (\"seed\", \"reset\")\n\
 Return a matrix with @code{gamma(@var{a},1)} distributed random elements.\n\
 The arguments are handled the same as the arguments for @code{rand},\n\
 except for the argument @var{a}.\n\
@@ -870,7 +890,9 @@
 @deftypefn  {Loadable Function} {} randp (@var{l}, @var{x})\n\
 @deftypefnx {Loadable Function} {} randp (@var{l}, @var{n}, @var{m})\n\
 @deftypefnx {Loadable Function} {} randp (\"state\", @var{x})\n\
address@hidden {Loadable Function} {} randp (\"state\", \"reset\")\n\
 @deftypefnx {Loadable Function} {} randp (\"seed\", @var{x})\n\
address@hidden {Loadable Function} {} randp (\"seed\", \"reset\")\n\
 Return a matrix with Poisson distributed random elements with mean value\n\
 parameter given by the first argument, @var{l}.  The arguments\n\
 are handled the same as the arguments for @code{rand}, except for the\n\