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