function rand_resetter = make_rand_resetter()
%MAKE_RAND_RESETTER returns a function to return all RNGs to their current state
%
% Usage:
%     rand_resetter = make_rand_resetter();
%     ... stuff ...
%     rand_resetter();
%     % Now random numbers are generated as though stuff didn't happen.
%
% The states of all of Matlab/Octave's built in RNGs are reinstated
% by rand_resetter().
%
% This is necessary for robust library routines. One could do this:
%     state = rand('twister');
%     ... stuff ...
%     rand('twister', state);
% But this requires knowing that 'twister' was the correct generator. Matlab
% provides no documented way of probing the generator that is in use. And randn
% needs to be updated too. This function sorts all that out.
%
% Since writing this, I discovered that Matlab 2008b introduces more complicated
% random number state functionality through "streams":
% http://blogs.mathworks.com/loren/2008/11/05/new-ways-with-random-numbers-part-i/
% Only use this function if your code uses the "classic" interface for RNG
% states, because I haven't got around to taking the new stuff into account, and
% (at time of writing) I didn't have Matlab 2008b.

% Iain Murray, October 2008

% May need to add to lists of generator types for future Matlab releases:
fn1 = make_family_resetter(@rand, {'seed', 'state', 'twister'});
fn2 = make_family_resetter(@randn, {'seed', 'state'});

% Return a function that will run both resetters:
rand_resetter = @() run_all({fn1, fn2});

function fn = make_family_resetter(randfn, generators)
% Make a RNG resetter for one family of generators (e.g. rand or randn)

% Save the states of all generators in this family:
% -------------------------------------------------
% one-liner version isn't robust to errors:
%states = cellfun(@(gg) randfn(gg), generators, 'UniformOutput', false);
states = cell(size(generators));
keep = logical(ones(size(generators)));
for ii = 1:numel(generators)
    try
        states{ii} = randfn(generators{ii});
    catch
        % Assume using old version of Matlab that doesn't know this generator
        keep(ii) = false;
    end
end
states = states(keep);
generators = generators(keep);

% bug in Octave 3.0.2 doesn't allow the one-liner:
%resetters = cellfun(@(gg, ss) @() randfn(gg, ss), generators, states, 'UniformOutput', false);
resetters = cell(size(generators));
for ii = 1:numel(generators)
    resetters{ii} = @() randfn(generators{ii}, states{ii});
end

% We work out which generator is in use by comparing some random draws with
% those generated by each generator with their saved states.
sampler = @() randfn(8, 1);
next = sampler();

for ii = 1:numel(generators)
    gg = generators{ii};
    resetters{ii}();
    next_with_gg = sampler();
    if isequal(next_with_gg, next)
        resetters{ii}();
        % Make the last resetter be the one for the current generator so that it
        % will be reinstated.
        last = resetters{end};
        resetters{end} = resetters{ii};
        resetters{ii} = last;

        fn = @() run_all(resetters);
        return;
    end
end

error('Matlab must have implemented a new random generator since this was written.');
% If you get this error, update the list of generators {'seed', 'state', 'twister'}
% at the top of this m-file