function C = nextperm(N,K)
%NEXTPERM Loop through permutations without replacement.
% NEXTPERM(N,K), when first called, returns a function handle. This
% function handle, when called, returns the next permutation without
% replacement of K elements taken from the set 1:N. This can be useful
% when the number of such permutations is too large to hold in memory at
% once. If the number of such permutations is not too large, use
% COMBINATOR(N,K,'p') (on the FEX) instead.
%
% The number of permutations without replacement is: N!/(N-K)!
% where N >= 1, N >= K >= 0
%
% Examples:
%
% % To use each permutation one at a time, put it in a loop.
% N = 4; % Length of the set.
% K = 3; % Number of samples taken for each sampling.
% H = nextperm(N,K);
% for ii = 1:((prod(1:N)/(prod(1:(N-K)))))
% A = H();
% % Do stuff with A: use it as an index, etc.
% end
%
%
% % To build all of the permutations, do this (See note below):
% ROWS = (prod(1:N)/(prod(1:(N-K))));
% C = ones(ROWS,K);
% for ii = 1:ROWS
% C(ii,:) = H();
% end
% %Note this is a lot slower than using combinator(N,K,'p')
%
% The function handle will cycle through when the final permutation is
% returned.
%
% See also, nchoosek, perms, combinator, npermutek (both on the FEX)
%
% Author: Matt Fig
% Contact: popkenai@yahoo.com
% Date: 6/9/2009
% Reference: http://mathworld.wolfram.com/BallPicking.html
% Arg checking.
if K>N
error('K must be less than or equal to N.')
end
if isempty(N) || K == 0
C = [];
return
elseif numel(N)~=1 || N<=0 || ~isreal(N) || floor(N) ~= N
error('N should be one real, positive integer. See help.')
elseif numel(K)~=1 || K<0 || ~isreal(K) || floor(K) ~= K
error('K should be one real non-negative integer. See help.')
end
CNT = 0; % Initializations for the nested function.
G = [];
TMP = [];
op = 1; % These are the variables which are passed to subfunc.
blk = 1; % Keeps track of the permutation blocks in the subfunc.
idx = 1:N/2; % Index vectors for subfunc.
idx2 = idx;
WV = 1:K; % Initial WV.
lim = 0; % Sets the limit for working index.
inc = K; % Controls which element of WV is being worked on.
cnt = 1; % Keeps track of which block we are in. See loop below.
BC = prod(1:N)/(prod(1:(N-K))); % Number of permutations.
L = prod(1:K); % Size of the blocks.
CNT = 0; % Counter for blocks.
Floop = BC - prod(1:K); % Length of blocks.
A2 = (N-K+1):N; % Seed for the final block.
B2 = 1:K;
ii = 1;
if K==1
C = @nestfunc2; % Return argument is a func handle.
else
C = @nestfunc;
end
function B = nestfunc2
% A handle to this function is passed back from the main func if K=1.
if WV > N
B = 1;
WV = 2;
return
end
B = WV(1);
WV = WV + 1;
end
function B = nestfunc
% A handle to this function is passed back from the main func.
if ii<=Floop
if CNT == 0
for jj = 1:inc
WV(K + jj - inc) = lim + jj;
end
% This is the first combination. We will permute it
% below for the rest of the calls in this block.
B = WV;
cnt = cnt + L;
if lim<(N-inc)
inc = 0; % Reset this guy.
end
TMP = WV; % This serves as seed for perm index.
inc = inc+1; % Increment the counter.
lim = WV(K+1-inc); % Limit for working index.
CNT = 1;
G = 1:K; % Seed for nextp subfunc
op = 1;
blk = 1;
idx = 1:N/2;
idx2 = idx;
ii = ii + 1; % "Loop" index.
else
% Permute the seed.
[G,op,blk,idx,idx2] = nextp(G,op,blk,idx,idx2,K);
B = TMP(G); % Index into current combination.
CNT = CNT + 1;
if CNT==(L) % Goes back to if for next seed (combin).
CNT = 0;
end
ii = ii + 1;
end
else
if ii == Floop+1 % We are at the last block
op = 1; % Re-initialize for subfunc.
blk = 1;
idx = 1:N/2;
idx2 = idx;
B = A2; % Seed for this last block.
cnt = cnt + 1;
ii = ii + 1;
return
elseif ii == BC + 1 % Time to start over.
WV = 1:K; % Seed for first block.
op = 1; % Re-initialize for subfunc.
blk = 1;
idx = 1:N/2;
idx2 = idx;
inc = 1; % Reset the incrementer.
lim = WV(K+1-inc); % And the lim.
cnt = L + 1; % Reset block counter.
CNT = 1;
ii = 2; % Reset "Loop" index.
B = WV;
TMP = WV;
B2 = 1:K;
G = 1:K;
return
end
% Permute seed the seed.
[B2,op,blk,idx,idx2] = nextp(B2,op,blk,idx,idx2,K);
B = A2(B2);
cnt = cnt + 1;
ii = ii + 1;
end
end
end
function [x,op,blk,idx,idx2] = nextp(x,op,blk,idx,idx2,K)
% Delivers one permutation at a time. This is a modification of an
% algorithm attributed to H. F. Trotter.
% x = 1:4;
% op = 1;
% blk = 1;
% idx = 1:length(x)/2;
% idx2 = idx;
% C(1,:) = x;
% for jj = 2:factorial(length(x))
% [x,op,blk,idx,idx2] = nextp(x,op,blk,idx,idx2,length(x));
% C(jj,:) = x;
% end
if op<K
np = op + 1; % Index of where the 1 goes.
x(op) = x(np); % Here we are just doing the switcheroo on adjacents.
x(np) = 1; % np is the current position of the 1.
op = np; % op is the old position of the 1, for next time.
return
else
x(K) = x(1); % Here we are switching the endpoints of the vect.
x(1) = 1;
op = 1; % Reset to the first position.
blk = blk + 1; % Keep track of this block. Reset every N*(N-1)th call.
if blk<K % Not through with this block yet.
return
end
blk = 1; % If here, we need to mix internal elements of the vect.
low = 2; % Start with the second element.
upp = K - 1; % And the next to last element.
while true % We always break this loop.
% In here, on every N*(N-1)th call, we are mixing internal elems.
% The number of times through the while loop is usually very small,
% even for large N most calls will go through less than 3 times.
cur = idx(low); % Holds the current index to switch.
if cur==upp
np = low;
idx2(low) = idx2(low) + 1;
else
np = cur + 1; % For next iter.
end
tmp = x(cur); % Switcheroo.
x(cur) = x(np);
x(np) = tmp;
idx(low) = np;
if idx2(low)<upp
break % Out of while loop, we've mixed them enough.
end
idx2(low) = low;
low = low + 1; % These two march towards each other.
upp = upp - 1;
end
end
end
