function y = rand_pair(a,b,n,varargin)
% RAND_PAIR Uniformly-sampled random pair combinations
%
% RAND_PAIR(A,B,N) is a 2-column matrix of length N, where each row is a unique pair of numbers
% chosen at random from the arrays A and B.
%
% RAND_PAIR(A,B,N,S) starts the pseudo-random number generator using seed S.
%
% RAND_PAIR(A,B,N,S,'bi') ensures that each pair is unique in both directions
% (i.e. if [2,3] already chosen then [3,2] rejected). If random seed not required then
% use the empty matrix [] instead of S.
%
% RAND_PAIR(....,'diff') ensures that no pair is chosen which consists of the same numbers (i.e. [2,2] rejected);
% set S and 'bi' flag to [] if not required.
%
% Mark Humphries 15/10/2004
if nargin > 3 & ~isempty(varargin{1})
rand('state',varargin{1});
end
if n > length(a) * length(b)
error('Cannot produce that many permutations of pairs')
end
if nargin == 6 & strcmp('diff',varargin{3})
diff = 1;
else
diff = 0;
end
y = zeros(n,2);
if nargin >= 5 & strcmp('bi',varargin{2})
% bidirectional uniqueness - generate pair first then check
for loop = 1:n
isunique = 0;
while ~isunique
first_idx = ceil(rand * length(a));
second_idx = ceil(rand * length(b));
% check for uniqueness if either (a) not checking for sameness OR (b) are checking, and not the same
if ~diff | (diff & a(first_idx) ~= b(second_idx))
% test both combinations
n1 = find(y(:,1) == a(first_idx));
exists = find(b(second_idx) == y(n1,2)); % is second number in second column at same row number?
if exists
% this pair exists
isunique = 0;
else
% otherwise test reversed pair
n2 = find(y(:,2) == a(first_idx));
exists = find(b(second_idx) == y(n2,1)); % is second number in *first* column at same row number?
if exists
isunique = 0;
else
% pair does not exist in either configuration
isunique = 1;
end
end
end
end
y(loop,:) = [a(first_idx) b(second_idx)];
end
elseif diff
% sameness checking without bidirection checking...
% check for fast version if pairs drawn from same distribution and
% there are many links to find
if length(a)==length(b) & a == b % & n > 8000
% all possible pair numbers
possible_pair_nos = 1:length(a)^2;
% pair numbers which are duplicates
duplicates = 1:length(a)+1:length(a)^2;
% remove duplicates
possible_pair_nos(duplicates) = [];
% create random permutation sequence
seq = randperm(length(possible_pair_nos));
% select first n of that sequence to create pairs
selected_pairs = possible_pair_nos(seq(1:n));
first_idx = ceil(selected_pairs ./ length(b));
second_idx = rem(selected_pairs,length(b));
second_idx(second_idx==0) = length(b);
y = [a(first_idx)',b(second_idx)'];
else
old_pairs = zeros(n,1);
% number of unique permutations
perms = length(a) * length(b);
% fastish version for no reverse duplicate checking
for loop = 1:n
isunique = 0;
while ~isunique
pair_no = ceil(rand * perms); % assumes pairs generated thus [a1,b1; a1,b2; a1,b3; a2,b1; a2,b2....]
if all(old_pairs ~= pair_no)
isunique = 1;
first_idx = ceil(pair_no / length(b)); % index for first column
second_idx = rem(pair_no,length(b));
if second_idx == 0 second_idx = length(b); end
first_num = a(first_idx);
second_num = b(second_idx); % number for second column
if first_num == second_num % if same then no longer unique
is_unique = 0;
end
end
end
old_pairs(loop) = pair_no;
y(loop,:) = [first_num second_num];
end
end % version if..then...
else
old_pairs = zeros(n,1);
% number of unique permutations
perms = length(a) * length(b);
% fast version for no reverse duplicate or sameness checking
for loop = 1:n
isunique = 0;
while ~isunique
pair_no = ceil(rand * perms); % assumes pairs generated thus [a1,b1; a1,b2; a1,b3; a2,b1; a2,b2....]
if all(old_pairs ~= pair_no)
isunique = 1;
end
end
old_pairs(loop) = pair_no;
first_idx = ceil(pair_no / length(b)); % number for first column
second_idx = rem(pair_no,length(b));
if second_idx == 0 second_idx = length(b); end
first_num = a(first_idx);
second_num = b(second_idx); % number for second column
y(loop,:) = [first_num second_num];
end
end