classdef RandStream < matlab.mixin.SetGetExactNames & handle
    %RandStream  GPU random number streams
    %   parallel.gpu.RandStream provides an interface for controlling the
    %   properties of one or more random number streams to be used on the
    %   GPU. Usage is the same as RandStream with the following
    %   restrictions:
    %    * Only the "combRecursive" (MRG32K3A), "Philox" and "Threefry" 
    %      generators are supported.
    %    * The "Inversion" normal transform is supported for all generators.
    %    * The 'BoxMuller' normal transform is supported only for the
    %      "Philox" and "Threefry" generators.
    %    * Setting the substream property is not allowed.
    %
    %   Examples:
    %
    %   % Create a single stream and designate it as the current global stream:
    %   s = parallel.gpu.RandStream("CombRecursive","Seed",1);
    %   parallel.gpu.RandStream.setGlobalStream(s);
    %
    %   % Create three independent streams:
    %   [s1,s2,s3] = parallel.gpu.RandStream.create("CombRecursive","NumStreams",3);
    %   r1 = rand(s1,100000,1);
    %   r2 = rand(s2,100000,1);
    %   r3 = rand(s3,100000,1);
    %   corrcoef(gather([r1,r2,r3]))
    %
    %   RandStream methods:
    %       RandStream/RandStream - Create a random number stream
    %       create                - Create multiple independent random number streams
    %       list                  - List available random number generator algorithms
    %       getGlobalStream       - Get the current global random number stream
    %       setGlobalStream       - Replace the global random number stream
    %       reset                 - Reset a stream to its initial internal state
    %       rand                  - Pseudo-random numbers from a uniform distribution
    %       randn                 - Pseudo-random numbers from a standard normal distribution
    %       randi                 - Pseudo-random integers from a uniform discrete distribution
    %       randperm              - Random permutation
    %
    %   See also gpurng,
    %            gpuArray/rand,
    %            gpuArray/randi,
    %            gpuArray/randn,
    %            gpuArray/randperm
    
    %   Copyright 2011-2020 The MathWorks, Inc.
    
    %% Public properties
    properties (Access = public)
        NormalTransform;              % Which method is being used to convert uniform into normal
        StreamIndex = uint64(1)       % The index of this stream
    end
    
    %% Read-only properties
    properties (SetAccess = private)
        Antithetic = false      % Should values be inverted (i.e. subtracted from 1 in the case of uniform)
        FullPrecision = true    % Whether to use two samples to ensure all bits of a double are set
        NumStreams = uint64(1)  % Number of separate streams being sampled
        Type = 'Threefry4x64-20'% The Random Number Generator algorithm being used
    end
    
    %% Dependent properties
    properties (Dependent, Transient)
        State                   % The current state of the RandStream. Note that the size of this state vector depends on the generator chosen.
        Seed                    % Seed used to initialize
    end
    
    %% Hidden properties (should be private but are needed elsewhere)
    properties (Hidden)
        CurrentState = uint32([])         % Current generator state
        InitialState = uint32([])         % Generator state on initialization (to allow reset)
    end
    
    %% Hidden constants (should be private but are needed elsewhere)
    properties(Hidden, Constant)
        DefaultGenerator = 'Threefry4x64-20'
    end
    
    %% Constants
    properties(GetAccess=private, Constant)
        ValidGenerators = {
            'MRG32K3A', 'MATLAB:RandStream:list:mrg32k3aDescription'
            'Philox4x32_10', 'MATLAB:RandStream:list:philox4x32_10Description'
            'Threefry4x64_20', 'MATLAB:RandStream:list:threefry4x64_20Description'
            }
        ValidRandNAlgs = {'Inversion', 'BoxMuller'}
    end

    %% Private properties
    properties(Access=private)
        Seed_ = uint32(0)        % Seed used to initialize
    end

    
    %% Public methods
    methods
        
        function obj = RandStream( name, varargin )
%RANDSTREAM  Create a GPU random number stream.
%   S = parallel.gpu.RandStream('GENTYPE') creates a random number stream
%   that uses the uniform pseudorandom number generator algorithm specified
%   by GENTYPE. Type "parallel.gpu.RandStream.list" for a list of possible
%   values for GENTYPE, or see <a href="matlab:helpview([docroot '\techdoc\math\math.map'],'choose_random_number_generator')">Choosing a Random Number Generator</a> 
%   for details on these generator algorithms.
%
%   Once you have created a random stream, you can use
%   parallel.gpu.RandStream.setGlobalStream to make it the global stream,
%   so that gpuArray.rand, gpuArray.randi, and gpuArray.randn draw values
%   from it.
%
%   [ ... ] = parallel.gpu.RandStream('GENTYPE','PARAM1',val1,'PARAM2',val2,...) 
%   allows you to specify optional parameter name/value pairs to control
%   creation of the stream.  Parameters are:
%
%      Seed            - a non-negative scalar integer seed with which to
%                        initialize the stream, or 'shuffle' to create a seed
%                        based on the current time.  Default is 0.
%      NormalTransform - the transformation algorithm that RANDN(S, ...) uses
%                        to generate normal pseudorandom values from uniform
%                        pseudorandom values.  Currently the 'Inversion'
%                        method is supported for all generators on the GPU.
%                        The 'BoxMuller' transform is only supported for
%                        the 'Philox4x32-10' and 'Threefry4x64-20'
%                        generators on the GPU.
%
%   The additional parameters "Antithetic" and "FullPrecision" cannot
%   currently be modified.
%
%   Streams created using parallel.gpu.RandStream may not be independent
%   from each other. Use parallel.gpu.RandStream.create to create multiple
%   streams that are independent.
%
%   Examples:
%
%      Create a random number stream, make it the global stream, and save and
%      restore its state to reproduce the output of RANDN:
%         s = parallel.gpu.RandStream('mrg32k3a');
%         parallel.gpu.RandStream.setGlobalStream(s);
%         savedState = s.State;
%         z1 = randn(1,5)
%         s.State = savedState;
%         z2 = randn(1,5) % z2 contains exactly the same values as z1
%
%      Return RAND, RANDI, and RANDN to their default startup settings:
%         s = parallel.gpu.RandStream('mrg32k3a','Seed',0)
%         parallel.gpu.RandStream.setGlobalStream(s);
%
%      Replace the current global random number stream with a stream whose
%      seed is based on the current time, so RAND, RANDI, and RANDN will
%      return different values in different MATLAB sessions.  NOTE: It is
%      usually not desirable to do this more than once per MATLAB session.
%         s = parallel.gpu.RandStream('mrg32k3a','Seed','shuffle');
%         parallel.gpu.RandStream.setGlobalStream(s);
%
%   See also gpurng,
%            parallel.gpu.RandStream.create,
%            parallel.gpu.RandStream.list,
%            parallel.gpu.RandStream.getGlobalStream,
%            parallel.gpu.RandStream.setGlobalStream,
%            parallel.gpu.RandStream/rand,
%            parallel.gpu.RandStream/randi,
%            parallel.gpu.RandStream/randn,
%            parallel.gpu.RandStream/randperm.

            if nargin < 1
                error(message('MATLAB:RandStream:TooFewInputs'));
            end
            obj.Type = name;
            
            % set the default normal transformation
            if strcmp(obj.Type,'MRG32K3A')
                obj.NormalTransform = obj.ValidRandNAlgs{1};
            else
                obj.NormalTransform = obj.ValidRandNAlgs{2};
            end
                
            if nargin>1
                params = varargin(1:2:end);
                values = varargin(2:2:end);
                if numel( params ) ~= numel( values )  ...
                        || any( ~cellfun( @matlab.internal.datatypes.isScalarText, params ) )
                    error(message('parallel:gpu:randstream:BadSyntax'));
                end
                % Check that all names are valid parameters
                parallel.gpu.RandStream.validateParamName( params );
                % Now set the remaining params on this stream. We have to
                % special-case Seed since we use a private version. We also
                % trap any errors so that they appear to come from the
                % constructor, not some helper function.
                try
                    for ii=1:numel( params )
                        if strcmpi( params{ii}, 'Seed' )
                            obj.Seed_ = values{ii};
                        else
                            obj.(params{ii}) = values{ii};
                        end
                    end
                catch err
                    throw(err);
                end
                
                % Check that the RandStream was not created with
                % non-compatible generator and normal transform.
                validateNormalTransform(obj,obj.NormalTransform);
                
            end
        end

    end % Public methods
    
    
    %% Property-access methods
    methods
        function s = get.State( obj )
            if isempty( obj.InitialState ) || isempty( obj.CurrentState )
                % Call a GPU function to force-initialize the state
                parallel.internal.gpu.randInitState( obj );
            end
            s = prepareStateForUser(obj.Type, obj.CurrentState(:), obj.InitialState(:));
        end
        
        function set.State( obj, s )
            % Check that the user-provided state is valid and throw an
            % appropriate error if not.
            checkUserState( obj.Type, s );
            % Before setting the state we must ensure it is properly
            % initialized. The "get" method will do this for us.
            [~] = obj.State;
            % Often the state contains initialization info which shouldn't
            % be updated.
            obj.CurrentState = interpretUserState(obj.Type, s);
        end
        
        function set.Seed( ~, ~ )
            % Throw a custom error for trying to set the seed as there is
            % an alternative way to do it
            error( message( 'MATLAB:RandStream:subsasgn:IllegalSeedAssignment' ) );
        end
        
        function set.Seed_( obj, value )
            if matlab.internal.datatypes.isScalarText( value ) && strcmpi( value, 'shuffle' )
                seed = RandStream.shuffleSeed();
            elseif isnumeric( value ) && isscalar( value ) && isreal( value ) ...
                    && (value >= 0) && (value < 4294967296) %2^32
                seed = cast(floor(value), 'like', uint32(1)); % Make sure we cast away sparsity, GPU-ness, etc.
            else
                error( message( 'MATLAB:RandStream:create:BadSeed' ) );
            end
            obj.Seed_ = seed;
        end
        
        function value = get.Seed( obj )
            value = obj.Seed_;
        end
        
        function set.StreamIndex( obj, value )
            if ~isnumeric(value) || ~isreal(value) || ~isscalar(value)
                error(message('MATLAB:RandStream:create:BadStreamIndex'));
            end
            obj.StreamIndex = cast(value, 'like', uint64(1)); % Make sure we cast away sparsity, GPU-ness, etc.
        end
        
        function set.Type( obj, value )
            % Check the value is valid, and convert to the canonical
            % algorithm name.
            if isstring(value) && isscalar(value)
                value = char(value);
            elseif ~ischar(value)
                error( message( 'MATLAB:RandStream:create:InvalidRNGType' ) );
            end
            % For backwards compatibility, treat "-" in the name as "_"
            value = strrep(value, "-", "_");
            
            algValue = RandStream.algName( value );
            idx = find( strcmpi( obj.ValidGenerators(:,1), algValue ) );
            if isempty( idx )
                error(message('parallel:gpu:randstream:CreateUnknownRNGType', value));
            else
                % Make sure we have the exact matching name
                algValue = obj.ValidGenerators{idx(1),1};
            end
            obj.Type = algValue;
        end
        
        function set.NormalTransform( obj, value )
            if ~ischar( value )
                error(message('MATLAB:RandStream:set:InvalidRandnAlg'));
            end
            whichAlg = strcmpi( value, obj.ValidRandNAlgs );
            if ~any( whichAlg )
                error(message('MATLAB:RandStream:set:UnknownRandnAlg', value));
            end
            alg = obj.ValidRandNAlgs{whichAlg};
            validateNormalTransform(obj,alg);
            obj.NormalTransform = alg; 
            % As and when this can be changed on the fly we will need to
            % update persistent storage
        end
    end
    
    %% Public static methods
    methods (Static)
        list();
        varargout = create( name, varargin );
        oldstream = setGlobalStream( stream );
        stream = getGlobalStream();
        stream = createDefaultStream();
    end % Public static methods
    
    %% Private methods
    methods (Access=private)
        out = isGlobal( obj );
    end % Private methods
    
    %% Private static methods
    methods (Access=private, Static)
        validateParamName( name );
    end % Private methods
    
    %% Hidden private methods
    methods(Access=private, Hidden)
        function validateNormalTransform(obj,normalTransform)
            % Check that the RandStream has a normal transform compatible
            % with the generator.
            if ( strcmp(obj.Type,'MRG32K3A') && strcmp(normalTransform,'BoxMuller') )
                error(message('parallel:gpu:randstream:InvalidNormalTransform',normalTransform,obj.Type));
            end
        end
    end
    
    %% Hidden public methods
    methods(Access=public, Hidden)
        % Methods that we inherit from base handle class, but do not want
        function a = fields(varargin),          throwUndefinedError; end %#ok<STOUT>
        function a = lt(varargin),              throwUndefinedError; end %#ok<STOUT>
        function a = le(varargin),              throwUndefinedError; end %#ok<STOUT>
        function a = ge(varargin),              throwUndefinedError; end %#ok<STOUT>
        function a = gt(varargin),              throwUndefinedError; end %#ok<STOUT>
        function a = permute(varargin),         throwUndefinedError; end %#ok<STOUT>
        function a = reshape(varargin),         throwUndefinedError; end %#ok<STOUT>
        function a = transpose(varargin),       throwUndefinedError; end %#ok<STOUT>
        function a = ctranspose(varargin),      throwUndefinedError; end %#ok<STOUT>
        function [a,b] = sort(varargin),        throwUndefinedError; end %#ok<STOUT>
        
        % Inherit default EQ, NE, ISVALID, FIELDNAMES, FINDPROP,
        % ADDLISTENER, NOTIFY from base handle class
        
        % All of these have to be taken away because they can create
        % non-scalar or empty arrays of objects.
        function a = findobj(varargin),         throwUndefinedError; end %#ok<STOUT>
        function a = cat(varargin),             throwNoCatError();   end %#ok<STOUT>
        function a = horzcat(varargin),         throwNoCatError();   end %#ok<STOUT>
        function a = vertcat(varargin),         throwNoCatError();   end %#ok<STOUT>
        
        % These are allowable but we don't really want to see them in the
        % method list (so that our methods match RandStreams as closely as
        % possible).
        function getdisp(varargin), getdisp@matlab.mixin.SetGetExactNames(varargin{:}); end
        function setdisp(varargin), setdisp@matlab.mixin.SetGetExactNames(varargin{:}), end
        function delete(varargin), end
    end
    
end % classdef


%% Helper functions

function throwNoCatError()
throwAsCaller(MException(message('MATLAB:RandStream:NoCatAllowed')));
end

function throwUndefinedError()
st = dbstack;
name = regexp(st(2).name,'\.','split');
throwAsCaller(MException(message('MATLAB:RandStream:UndefinedFunction', name{2})));
end

function checkUserState( name, state )
% Perform various checks on the supplied state vector
switch upper(name)
    case 'MRG32K3A'
        requiredStateLength = 12;
        requiredStateType = 'uint32';
        valueTest = @checkMRG32K3AStateValues;
    case 'PHILOX4X32_10'
        requiredStateLength = 7;
        requiredStateType = 'uint32';
        valueTest = @checkPhilox4x32StateValues;
    case 'THREEFRY4X64_20'
        requiredStateLength = 17;
        requiredStateType = 'uint32';
        valueTest = @checkThreeFry4x64StateValues;
    otherwise
        assert(false,'Should never get here - managed to have an unsupported RNG type');
end

% Check the class and length.
if ~isa( state, requiredStateType )
    throwAsCaller(MException(message('MATLAB:RandStream:set:InvalidStateClass', name)));
end
if numel(state) ~= requiredStateLength
    throwAsCaller(MException(message('MATLAB:RandStream:set:InvalidStateLen', name)));
end
% Test the values for any generator-specific caps or restrictions
if ~isempty(valueTest)
    err = valueTest(state);
    if ~isempty(err)
        throwAsCaller(MException(message(err)));
    end
end
end

function err = checkMRG32K3AStateValues( state )
% Check the generator state to make sure it is valid.

% We can simply try constructing a client-side RandStream (this is pretty
% fast).
err = '';
try
    rs = RandStream('mrg32k3a','Seed',0);
    rs.State = state;
catch err
    err = err.identifier;
end

end % checkMRG32K3AStateValues

function err = checkPhilox4x32StateValues( state )
% Check the generator state to make sure it is valid.
%
% The Philox4x32 state is composed of six 32-bit values that can
% take their full range. The final element is a marker of where we are
% within a sample if a partial sample has been taken. It therefore
% must have a value less than two (since we get two doubles per
% sample). We must ignore the final bit though.
err = '';
if mod(state(end),2^31) >= 2
    err = 'MATLAB:RandStream:philox4x32_10:InvalidState';
end
end % checkPhilox4x32StateValues

function err = checkThreeFry4x64StateValues( state )
% Check the generator state to make sure it is valid.
%
% The ThreeFry4x64 state is composed of eight 64-bit values that can
% take their full range. The final element is a marker of where we are
% within a sample if a partial sample has been taken. It therefore
% must have a value less than four (since we get four doubles per
% sample). We must ignore the final bit though.
err = '';
if mod(state(end),2^31) >= 4
    err = 'MATLAB:RandStream:threefry4x64_20:InvalidState';
end
end % checkThreeFry4x64StateValues

function out = interpretUserState( type, state )
% Extract a current state from some user-supplied state vector
if strcmpi(type, 'MRG32K3a')
    % For MRG32K3a the user state contains both current and initial states. We
    % only need the "current" bit
    out = state(1:length(state)/2);
else
    % For all of the remaining supported generators the user state is just the
    % current state.
    out = state;
end
end % interpretUserState

function out = prepareStateForUser( type, currentstate, initialstate )
% Extract a current state from some user-supplied state vector
if strcmpi(type, 'MRG32K3a')
    % For MRG32K3a the user state contains both current and initial states.
    out = [currentstate; initialstate];
else
    % For all of the remaining supported generators the user state is just the
    % current state.
    out = currentstate;
end
end % prepareStateForUser