function result=cosmo_rand(varargin)
% generate uniform pseudo-random numbers, optionally using a seed value
%
% result=cosmo_rand(s1,...,sN,['seed',seed])
%
% Input:
%    s*              scalar or vector indicating dimensions of the result
%    'seed', seed    (optional) if provided, use this seed value for
%                    pseudo-random number generation
%
% Output:
%    result          array of size s1 x s2 x ... sN. If the seed option is
%                    used, repeated calls with the same seed and element
%                    dimensions gives the same result
% Example:
%     % generate 2x2 pseudo-random number matrices twice, just like 'rand'
%     % (repeated calls give different outputs)
%     x1=cosmo_rand(2,2);
%     x2=cosmo_rand(2,2);
%     isequal(x1,x2)
%     %|| false
%     %
%     % as above, but specify a seed; repeated calls give the same output
%     x3=cosmo_rand(2,2,'seed',314);
%     x4=cosmo_rand(2,2,'seed',314);
%     isequal(x3,x4)
%     %|| true
%     %
%     % using a different seed gives a different output
%     x5=cosmo_rand(2,2,'seed',315);
%     isequal(x3,x5)
%     %|| false
%
%
% Notes:
%   - this function behaves identically to the builtin 'rand' function,
%     except that it supports a 'seed' option, which allows for
%     deterministic pseudo-number generation
%   - when using the 'seed' option, this function gives identical output
%     under both matlab and octave. To achieve this, the PRNG is set to a
%     different state for the two platforms
%   - this function uses the Mersenne twister algorithm by default, even
%     when 'seed' is used (unlike Matlab and Octave).
%
% #   For CoSMoMVPA's copyright information and license terms,   #
% #   see the COPYING file distributed with CoSMoMVPA.           #

    [sizes,seed,class_func]=process_input(varargin{:});

    randomizer=@rand; % default
    if seed~=0
        is_matlab=cosmo_wtf('is_matlab');

        if is_matlab
            rng_state=get_mersenne_state_from_seed(seed, is_matlab);

            stream=RandStream('mt19937ar','Seed',rng_state.Seed);
            stream.State=rng_state.State;

            randomizer=@stream.rand;
        else
            % preserve old PRNG state
            orig_rng_state=rand('state');
            cleaner=onCleanup(@()rand('state',orig_rng_state));

            % set random number generation state
            rng_state=get_mersenne_state_from_seed(seed, is_matlab);
            rand('state',rng_state);
        end
    end

    result=class_func(randomizer(sizes));


function rng_state=get_mersenne_state_from_seed(seed, is_matlab)
    % set the PRNG of the mersenne twister based on seed
    %
    % based on pseudo-code from wikipedia:
    % http://en.wikipedia.org/wiki/Mersenne_twister
    persistent cached_seed
    persistent cached_rng_state

    if isequal(cached_seed,seed)
        rng_state=cached_rng_state;
        return;
    end

    max_uint32=2^32-1;
    state=uint64(zeros(625,1));
    state(1)=bitand(uint64(seed),max_uint32);

    mersenne_mult=uint64(1812433253);

    for j=1:623
        v=mersenne_mult.*bitxor(state(j),bitshift(state(j),-30))+uint64(j);
        state(j+1)=bitand(v,max_uint32);
    end

    state(end)=1;

    if is_matlab
        % reverse counter relative to Octave
        % (this is undocumented in both Matlab and Octave)
        state(end)=uint64(625)-state(end);

        % matlab uses a struct to set the state
        rng_state=struct();
        rng_state.State=uint32(state);
        rng_state.Type='twister';
        rng_state.Seed=uint32(0);
    else
        % octave uses a vector to set the state
        rng_state=state;
    end

    cached_rng_state=rng_state;
    cached_seed=seed;

function x=identify_func(x)
    % do nothing

function [sizes,seed,class_func]=process_input(varargin)
    persistent cached_varargin;
    persistent cached_sizes;
    persistent cached_seed;
    persistent cached_class_func;
    if isequal(varargin,cached_varargin)
        sizes=cached_sizes;
        seed=cached_seed;
        class_func=cached_class_func;
        return;
    end

    n=numel(varargin);

    seed=0;
    sizes_cell=cell(1,n);
    class_func=[];

    has_processed_sizes=false;

    % process each argument
    k=0;
    while k<n
        k=k+1;
        arg=varargin{k};
        if isnumeric(arg)
            if has_processed_sizes
                error('size argument not allowed after string argument');
            end

            ensure_positive_vector(k,arg);
            sizes_cell{k}=arg(:)';

        elseif ischar(arg)
            has_processed_sizes=true;
            switch arg
                case 'single'
                    if ~isempty(class_func)
                        error('type can only be set once');
                    end
                    class_func=@single;

                case 'double'
                    if ~isempty(class_func)
                        error('type can only be set once');
                    end
                    class_func=@identify_func;

                case 'seed'
                    k=k+1;
                    if k>n
                        error('missing value after key ''%s''', arg);
                    end
                    value=varargin{k};
                    ensure_positive_scalar(k,value);
                    seed=value;

                otherwise
                    error('unsupported key ''%s''', arg);
            end

        else
            error('illegal input at position %d', k);
        end
    end

    sizes=[sizes_cell{:}];

    % no size provided, output is scalar
    if isempty(sizes)
        sizes=1;
    end

    if isempty(class_func)
        class_func=@identify_func;
    end

    cached_varargin=varargin;
    cached_sizes=sizes;
    cached_seed=seed;
    cached_class_func=class_func;

function ensure_positive_scalar(k,arg)
    ensure_positive_vector(k,arg);
    if ~isscalar(arg)
        error('argument at position %d is not a scalar',k);
    end

function ensure_positive_vector(k,arg)
    if ~isvector(arg) || ~isnumeric(arg) || ~(all(arg>=0))
        error('argument at position %d is not positive',k);
    end
