bufrngstreammaker(rngFun,blocksize,init)

function rngStream = bufrngstreammaker(rngFun,blocksize,init)
%BUFRNGSTREAMMAKER Create a buffered random number stream.
%   RNGSTREAM = BUFRNGSTREAMMAKER(RNGFUN,BLOCKSIZE,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 @, that accepts a size
%   vector.  BLOCKSIZE is the number of values maintained internally, and
%   determines how often RNGSTREAM must call RNGFUN.  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(SIZE) returns an array of random values, generated by
%      calling RNGFUN, of size SIZE.
%
%   RNGSTREAM maintains an internal buffer of BLOCKSIZE pseudorandom values,
%   so it does not call RNGFUN each time it is called.
%
%   Example:
%
%      % A RNG stream for uniform random values.
%      rng = bufrngstreammaker(@rand, 10000);
%      vals = rng([5,1]);
%
%      % A RNG stream for t random values with 5 d.f.
%      t5rng = bufrngstreammaker(@(size) trnd(5,size), 1000);
%      t5val = t5rng();

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

if nargin < 2
    error('Requires at least two inputs.');
end

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

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

try
    buf = rngFun([1,blocksize]);
catch
    if nargin > 1
        rand('state',state0); randn('state',staten0);
    end
    error('RNGFUN generated the following error: %s',lasterr);
end
bufptr = 1;

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

rngStream = @rngStreamTemplate;

    function r = rngStreamTemplate(size,newstate)
        if nargin > 0 && ischar(size) && strcmp(size,'state')
            if nargin == 1
                % Return the current state.
                r = {state staten};
            elseif nargin == 2
                if iscell(newstate) && length(newstate)==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 = newstate{1}; staten = newstate{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
        elseif nargin == 1
            len = prod(size);
            buflen = length(buf) - bufptr + 1;
            if len > buflen
                % Generate the random values and get the updated state.
                state0 = rand('state'); staten0 = randn('state');
                rand('state',state); randn('state',staten);
                try
                    buf = [buf(bufptr:end) rngFun([1,len-buflen+blocksize])];
                catch
                    rand('state',state0); randn('state',staten0);
                    error('RNGFUN generated the following error: %s',lasterr);
                end
                state = rand('state'); staten = randn('state');
                rand('state',state0); randn('state',staten0);
                bufptr = 1;
            end
            r = reshape(buf(bufptr:bufptr+len-1),size);
            bufptr = bufptr + len;
        else % faster special case for scalar value
            if bufptr > blocksize
                % Generate the random values and get the updated state.
                state0 = rand('state'); staten0 = randn('state');
                rand('state',state); randn('state',staten);
                try
                    buf = rngFun([1,blocksize]);
                catch
                    rand('state',state0); randn('state',staten0);
                    error('RNGFUN generated the following error: %s',lasterr);
                end
                state = rand('state'); staten = randn('state');
                rand('state',state0); randn('state',staten0);
                bufptr = 1;
            end
            r = buf(bufptr);
            bufptr = bufptr + 1;
        end
    end
end