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 first this, I discovered that Matlab 2008b introduced more
% complicated random number state functionality through "streams":
% http://blogs.mathworks.com/loren/2008/11/05/new-ways-with-random-numbers-part-i/
% This function now saves all the states of the classic generators and (if you
% have version >2000b) the default stream. Calling the function it returns will
% restore all the states and the default stream (again, if it exists).

% Iain Murray, October 2008, October 2009

% These are all the classic deprecated generators.
% I believe that they are unlikely to change:
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 the classic 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

try
    % This should work in Matlab >2008b
    stream = RandStream.getDefaultStream();
    saved_state = stream.State;
    resetters{end+1} = @() restore_stream(stream, saved_state);
    fn = @() run_all(resetters);
catch
    % If in older Matlab or Octave 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('Sorry, this routine is confused, possibly by a new version of Matlab.');
end

function restore_stream(stream, saved_state)
RandStream.setDefaultStream(stream);
stream.State = saved_state;