A set of functions for generating and testing random permutations of the integers (1,2, ..., n).
The main functions provided are:
1. p = GRPfys(n); generates a random permutation of length n.
2. p = GRPmex(n); generates a random permutation of length n.
3. p = GRDrej(n); generates a derangement of length n.
4. p = GRDmex(n); generates a derangement of length n.
5. p = GRDMPP(p); generates a derangement of length n.
6. p = GRPcyc(p); generates a cyclic permutation of length n.
7. Various interrogation functions, e.g., IsCyc(p).
Included is a set of notes on Random Permutations and their implementation and testing.
1.2 | Changed package name from RPGTools to RPGLab. Updated notes to RPGLabNotes1.1.pdf |
Jan Simon (view profile)
This is a neat collection of codes for combinatorics. The consequent mentioning of the referenced publications is extra-ordinary.
Splitting the code to a lot of tiny functions is a good programming practice, but inlinig the subfunction or even appending them as local functions in the same M-file is a little bit faster in Matlab. However, this does not reduce the quality of this submission.
Derek O'Connor (view profile)
Jan,
(1). No, your modification still uses 2 x p memory. So the "clear p" has to stay in p=GRPfys(n) to conserve memory.
(2). Is there a typo in "p=zeros(1, n); p(1) = 1;" ?
p = GRPfysSimon(5) gives [0 0 1 0 0], and moves the '1' randomly in subsequent calls.
All the generators I've written must start with a (any?) permutation p and then do random transpositions on p,
p = Trans(p,i,r), using the three-statement swap: t = p(i); p(i) = p(r); p(r) = t.
The correctness proofs of the algorithms assume this condition.
I do not use Knuth's two-statement short-cut for permutations of [1,2,...,n] and so GRPfys(p) works for any p. For example, if we want to generalize GRPfys to get a permutation of the integers [L,L+1,...,U] then the only change needed in GRPfys is replace
"p=I(n)" with "p = L:U; n = U-L+1" to get p = GRPfysLU(L,U). This works for
p = GRPfysLU(-5,5) and p = GRPfysLU(-6,-2).
Alternatively, "p = GRPfys(U-L+1)+L-1" does the same thing but is error-prone and puts a burden on the user.
(3). I'll put a link to Shuffle in the next update.
(4). GRPfysSimon(3) neatly simulates the Shell or Three-Card-Trick game.
Thanks for you interest,
Derek.
Jan Simon (view profile)
Sorry Derek, my message was too lean: I do not meant this single line, but the general idea.
function p=GRDrej(n)
p = GRPfys(n);
NotDer = HasFP(p);
while NotDer
p(:) = GRPfys(n); % <= here
NotDer = HasFP(p);
end
In GRPfys you do not have to create I(N) at first: "p=zeros(1, n); p(1) = 1;" does the job also - but is not remarkably faster, if I measured correctly.
Perhaps it would be helpful, if you add a link to Shuffle in the FEX.
Kind regards, Jan
Derek O'Connor (view profile)
"p(:)=GRPfys(n);" causes an assignment error:
// ??? In an assignment A(:) = B, the number of elements in A and B must be the same.
//
// Error in ==> GRDrej at 17
// p(:) = GRPfys(n);
Again, this shows that the semantics (meaning) of Matlab statements are not well-defined.
Here is another example of "chaotic" semantics and the "perils of vectorization", which I found this afternoon when writing and testing functions for inverting a permutation p:
function q = InvPCmp(p);
%
% Invert the permutation p so that p(q) = q(p) = I
%
% Derek O'Connor 31 Jan 2011.
%
n = length(p);
q = zeros(1,n);
for i = 1:n
q(p(i)) = i;
end;
return % InvPCmp(n)
%
% Here are two 'slick' one-line alternatives
%
% (1) InvPSrt: [t,q] = sort(p);
%
% (2) InvPVec: q(p) = 1:length(p);
%
Alternative (1) InvPSrt(p), uses the same trick that Matlab's randperm uses.
Alternative (2) InvPVec(p), is what all good Matlab-ers like: a vectorized version of the component-style function InvPCmp(p).
Here are the last lines of timing-memory test results:
InvPCmp: n = 2^29 tg = 195 ti = 70 ti/tg = 0.36 Mem = 8.0 GB
InvPSrt: n = 2^29 tg = 196 ti = 142 ti/tg = 0.72 Mem = 11.5 GB
InvPVec: n = 2^28 tg = 91 ti = 38 ti/tg = 0.42 Mem = 9.5 GB
n = 2^29 tg = 197 started page swapping Mem > 16.0 GB
You can see that the component version InvPCmp is the fastest and uses only the memory necessary to store two vectors p and q of length n = 2^29.
The sort version InvPSrt takes twice as long and uses almost 50% more memory than InvPCmp.
The vectorized version InvPVec was a complete disaster: it took longer on the smaller vectors and ran out of memory for n = 2^29.
Derek O'Connor
Jan Simon (view profile)
Does "p(:)=GRPfys(n);" reduce the memory foot print?
Derek O'Connor (view profile)
If what you say is true, and I believe it is, then this demonstrates very bad memory management by Matlab's interpreter, JIT accelerator, or something else. The very fact that it can be "fixed" so easily confirms this.
The semantics of some Matlab statements seem to be (increasingly?) chaotic. The only way I can determine the meaning (or behavior) of a Matlab statement is to run it, and then speculate about what has happened, as we are doing now. This puts a huge burden of uncertainty on the programmer which leads to bad programs and wasted time.
Derek O'Connor
Bruno Luong (view profile)
I guess the memory arises in the first statement
p = GRPfys(n);
The result of GRPfys(n) return as output of GRPfys(n), and right before the old p is released from memory and assigned by the output, there are temporary two copies: old p and new p.
Derek O'Connor (view profile)
I have noticed that the derangement generator GRPrej(n) uses 2 times the memory that it should, when the while loop is executed 2 or more times:
NotDer = true;
while NotDer
p = GRPfys(n);
NotDer = HasFP(p);
end;
For example, with n = 2^27, the vector p should be 1 GB. Using Windows Task Manager (WTM) I saw memory use increase by 1 GB on the first pass through the while-loop. The, on the second pass memory increased by another 1GB, and stayed at that level during subsequent passes. HasFP(p) does not alter p.
This problem does not arise in GRDmex(n) which uses Jan Simon's mexed Shuffle function.
I have fixed this problem but I don't know why the fix works, or indeed, why I need the fix at all.
NotDer = true;
while NotDer
clear p;
p = GRPfys(n);
NotDer = HasFP(p);
end;
Can anybody give an explanation? I'm using R2008a, 64 bit under Windows 7.
Derek O'Connor
Bruno Luong (view profile)
This set of permutation tools is very useful (derangement + cycling), and the algorithms used in the implementation are carefully selected.