random(varargin)

% RANDOM    Uniformly distributed random numbers on a non-continuous 
%                    domain defined  by several intervals
%
% It takes the advantage of RAND basic function to generate random entries, 
% chosen from a uniform distribution on the non-continuous domain defined 
% by several intervals.
%
% RANDOM  with no input arguments is the same as RAND.
% RANDOM(D)  returns one random number on the domain D, where D is (Nx2)-real 
%                      matrix defining N non-intersected intervals.
% RANDOM(D, N) returns an N-by-N matrix with random entries, chosen from
%                          a uniform distribution on the domain defined by D.
% RANDOM(D, M,N) and RANDOM(D, [M,N]) return M-by-N matrices with random
%                        entries on the domain defined by D.
% RANDOM(D, M,N,P,...) or RANDOM(D, [M,N,P,...]) generate random arrays on the domain
%                        defined by D
%
% See RAND function to display or reset current state of the uniform generator.
%
% See also RAND

% Version 1.1
% Alex Podgaetsky, November 2002
% [e-mail]
%
% Revisions:
%       Version 1.2  -  December, 2004
%                       Bug fix: Crash using RANDOM(D) was corrected
%                       Thanks to Amin.
%       Version 1.1  -  run time improving 
%       Version 1.0 -   initial version

function  Y = random(varargin)

if nargin == 0
    Y = rand; return;
end
S = size(varargin{1});
if length(S)>2  |   S(2) ~= 2
    error('First input argument must be (N x 2)-matrix!');
end
try
    varargin_2_end = [varargin{2:end}];
    if size(varargin_2_end,1)>1, error;  end
    if isempty(varargin_2_end) varargin_2_end = 1; end
catch
    error('Input arguments starting from the second one must be scalar!')
end
if ~all(isfinite([varargin{1}(:)' varargin_2_end ])) 
    error('NaN and Inf not allowed!');
end

if S(1) == 1
    if varargin{1}(1) == 0  &  varargin{1}(2)==1
        Y = rand(varargin_2_end);
    elseif  varargin{1}(1) == 0
        Y = varargin{1}(2) * rand(varargin_2_end);
    else
        Y = varargin{1}(1) + (varargin{1}(2)-varargin{1}(1)) * rand(varargin_2_end);
    end
    return;
end

if nargin == 1
    dim_size = [1,  1];
elseif nargin == 2    
    if length(varargin{2}) > 1
        dim_size = varargin{2};
    else
        dim_size = [varargin{2},  varargin{2}];
    end
else
    dim_size = varargin_2_end;
end
len = prod(dim_size);

Intervals = sort(varargin{1},2);
if any( diff(reshape(sortrows(Intervals)', 2*S(1), 1))  < 0 )
    error('Domain is badly defined!');
end

D = Intervals(:,2) - Intervals(:,1);
tmp = 1 + sum(repmat(cumsum(D/sum(D)),1,len) < repmat(rand(1,len), S(1),1),  1);
Y = reshape(Intervals(tmp,1)+rand(len,1).*D(tmp'), dim_size);
%%%%%   END OF FILE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%