rngstreammaker(rngFun,init)

function rngStream = rngstreammaker(rngFun,init)
%RNGSTREAMMAKER Create a random number stream.
%   RNGSTREAM = RNGSTREAMMAKER(RNGFUN,INIT) returns a function handle to a
%   pseudorandom number generator that is not affected by calls to other
%   random number functions, i.e., a random number stream.  RNGFUN is a
%   random number generator function, specified with @.  INIT is an integer
%   that determines the initial state, as described for RAND or RANDN.
%
%   RNGSTREAM returns random values, generated by calling RNGFUN, using the
%   following syntax:
%
%      S = RNGSTREAM('state') returns the current state of the generator.
%
%      SOLD = RNGSTREAM('state',SNEW) sets the state of the generator to SNEW,
%      and returns the previous state in SOLD.  SNEW must come from a previous
%      call to RNGSTREAM.
%
%      R = RNGSTREAM(...) returns an array of random values generated by
%      calling RNGFUN.  The input in this case typically include distribution
%      parameters and output size arguments, but may be anything that RNGFUN
%      is prepared to accept.
%
%   Examples:
%
%      % A RNG stream for uniform random values.
%      rng = rngstreammaker(@rand);
%      vals = rng([5,1]);
%
%      % A RNG stream for t random values with 5 d.f.
%      t5rng = rngstreammaker(@(size) trnd(5,size));
%      t5vals = t5rng([10,1]);
%
%      % A RNG stream for arbitrary gamma random values, initialized to state 0.
%      gamrng = rngstreammaker(@gamrnd,0);
%      gam12vals = gamrng(1,2,[10,1]);
%      gam23vals = gamrng(2,3,[1,10]);

%   Copyright 2009 The MathWorks, Inc.
%   Revision: 1.0  Date: 2004/12/01
%
%   Requires MATLAB R14.

if nargin < 1
    error('Requires at least one input.');
end

% Initialize rand and randn to specified state.
if nargin > 1
    state0 = rand('state'); staten0 = randn('state');
    try
        rand('state',init); randn('state',init);
    catch
        rand('state',state0); randn('state',staten0);
        error('Invalid initial state: %s',lasterr);
    end
end

% Get current state of rand and randn.
state = rand('state'); staten = randn('state');

% Put rand and randn back where they were.
if nargin > 1
    rand('state',state0); randn('state',staten0);
end


rngStream = @rngStreamTemplate;

    function r = rngStreamTemplate(varargin)
        if nargin > 0 && ischar(varargin{1}) && strcmp(varargin{1},'state')
            if nargin == 1
                % Return the current state.
                r = {state staten};
            elseif nargin == 2
                new = varargin{2};
                if iscell(new) && length(new)==2
                    % Return the existing state if asked for.
                    if nargout > 0, r = {state staten}; end;
                    % Update the state.
                    state0 = rand('state'); staten0 = randn('state');
                    state = new{1}; staten = new{2};
                    try
                        rand('state',state); randn('state',staten);
                    catch
                        rand('state',state0); randn('state',staten0);
                        error('Invalid state: %s',lasterr);
                    end
                    rand('state',state0); randn('state',staten0);
                else
                    error('Invalid state.');
                end
            else
                error('Wrong number of inputs.');
            end
        else
            % Generate the random values and get the updated state.
            state0 = rand('state'); staten0 = randn('state');
            rand('state',state); randn('state',staten);
            try
                r = rngFun(varargin{:});
            catch
                rand('state',state0); randn('state',staten0);
                error('Invalid inputs: %s',lasterr);
            end
            state = rand('state'); staten = randn('state');
            rand('state',state0); randn('state',staten0);
        end
    end
end