[SciPy-dev] New random capability [was: Release of scipy core...]

[Tom Loredo]

Robert,

Wow, thanks for what sounds like some well thought out improvements
in the "random" capability of scipy.  Working on this has been on
my todo list, but at lower priority than other things.  I think 
it's in more capable hands in any case!

> I also think some provision for setting a simple seed, perhaps just a long,
> is useful for testing when it is nice to have a repeatable sequences.

Actually, I think some users need more than this (I certainly speak
at least for myself!).  I think there should be capability to store
the full current state, and restore it.  This would let you not
only repeat an earlier calculation verbatim, but also pick up where
you left off, so to speak, guaranteeing you the properties of the
full chain through many subsequent runs, as well as the ability
to backtrack to some extent.

I've copied below the code I use for this now.  By default it stores
the state in the CWD in ".scipy-rng-state" as a pickle.  When it stores
the state, if there happens to already be a stored state file, it
renames it ".scipy-rng-state.1" (and if that exists, also it is renamed
"...2").  In many years of doing Monte Carlo simulations and
statistical calculations, there have been lots of times I realized I
wanted to rerun a simulation just before the one I just started, so
this backup capability has been very handy.

Though scipy's current RNG has very simple state, I wrote these
anticipating future RNG capability that might let one use
different generators with different state.  So the pickle is
actually a tuple, (id, state), with id a string identifying
the generator.  If indeed there are ever multiple RNGs, 
restore_rng should probably set the current generator to the
one whose state is restored, and not just restore the state
of a RNG that might not actually be used.

You'll probably come up with something better; the code is
attached merely as a model of the capability I'm talking about.

-Tom Loredo

~~~~~~~~~~~~~~~~~~~~~
import os, pickle
import os.path as path

# Note: use 'SSRand' since 'from scipy import *' would overwrite 'rand'
import scipy.stats.rand as SSRand

def restore_rng(fname='.scipy-rng-state'):
    """
    Restore the state of SciPy's RNG from the contents of a file in the CWD
    if the file exists; otherwise use a (fixed) default initialization.
    """
    if os.access(fname, os.R_OK | os.W_OK):
        rng_file = open(fname, 'r')
        id, state = pickle.load(rng_file)
        rng_file.close()
        print 'Recovered RNG state:', id, state
        if id == 'C2LCGX':
            SSRand.set_seeds(*state)
        else:
            raise ValueError, 'Invalid ID for RNG in %s!' % fname
    else:
        print 'No accessible RNG status file; using fixed default initialization.'
        SSRand.set_seeds(1234567890, 123456789)  # From L'Ecuyer & Cote (1991)
        
def save_rng(fname='.scipy-rng-state'):
    """
    Save the state of SciPy's RNG to a file in the CWD.  Backup the
    previous two saved states if present.
    If the RNG state file exists (from a previous save), rename the previous
    one with a '.1' suffix.  If a '.1' file exists, rename it with a 
    '.2' suffix.
    """
    id = 'C2LCGX'  # SciPy's RNG from L'Ecuyer & Cote (1991)
    state = SSRand.get_seeds()
    if path.exists(fname):
        fname1 = fname + '.1'
        if path.exists(fname1):
            fname2 = fname + '.2'
            os.rename(fname1, fname2)
        os.rename(fname, fname1)
    ofile = open(fname, 'w')
    pickle.dump((id, state), ofile)
    ofile.close()