NYCJUG/2011-05-10/DevelopmentOfSimpleSinglePermutationAlgorithm
Wherein we develop a simple algorithm, suitable for translation into a less array-oriented language like Java or Javascript, for permuting a vector randomly. In doing so, we come to appreciate the value of test-driven development and implement a useful way of associating a test case with a set of algorithms as well as discovering some subtleties of this fairly simple problem.
Some Experiments in J to Help Build a Simple Permuter
We want a way to randomize the order of a vector of items, like the "shuffle" function on a music player. A naive approach to this problem is to pick each item in turn on the basis of a random number. However, this gives no guarantee that we won't pick the same item twice in a row or otherwise soon after we've picked it once. A permutation of items, however, guarantees an even distribution of the items being shuffled. More background on the motivation for this may be found in this discussion at NYCJUG.
Since this is to be implemented in Javascript, we have to restrict the form of the algorithm to avoid major array operations (so we must loop unnecessarily from a J perspective) and the relevant J primitives more advanced than the random number generator restricted to using only the argument of zero in order to emulate the simple Javascript equivalent.
The code explained below is File:PermuteExperiments.ijs.
Initial Attempt
My first attempt was truly awful: it is overly elaborate as it loops through the vector and moves each element to a new position, keeping track of what's been moved already and adjusting the initially-chosen target to avoid moving an item more than once.
NB.* permuteExperiments.ijs: test various simple permutation algorithms.
permute0=: 3 : 0
p0=. <.(#y)*(#y)?@$0
posn=. _1$~#y
for_ii. i.#y do. next=. ii{p0 NB. Find closest if already taken:
if. _1~:next{posn do.
whpos=. _1=posn
possibs=. ((#y)*-.whpos)+whpos#^:_1]|next-I.whpos
next=. (]i.<./)possibs
end.
posn=. (ii{y) next}posn
end.
posn
)
Start of Test-Driven Development
However, this does immediately raise the question of how to verify if this works correctly. Initially, I simply ran it with a small, simple argument and eye-balled the results, something like this:
permute0 i.10 1 3 4 2 6 9 5 8 7 0 permute0 i.10 7 4 9 5 3 1 0 2 8 6 permute0 i.10 2 0 4 5 9 1 6 3 7 8
However, I quickly tired of this and also wondered how to verify the result for an argument too large to "eye-ball". This leads me to write my first test case to verify the goodness of a result.
permute0TC0_TC0_=: 3 : 0
assert. y (([: *./ e.~) *. ([: -. -:) ~: ] */ .= [: {. ]) permute0_base_ y
)
The name is my clumsy attempt to associate the algorithm being tested with a test case, allowing for the possibility that I might have multiple test cases. Also, I put the test case code into its own namespace "TC0" to avoid cluttering the working area.
This tests if the input and the result have all the same elements - ([: *./ e.~) - and the input is not the same as the output - ([: -. -:) ~: ] - except if the input consists only of identical elements. Strictly speaking, this second part of the test is incorrect. Mathematically, one permutation of a vector is itself. However, for the purposes of shuffling items to avoid repetitions, we want a permutation that differs from the original.
Note, too, that by taking an argument to be an arbitrary vector, we make things harder on ourselves than necessary. Eventually, I came to realize that I should only deal with index vectors, as in my "eye-ball" tests above.
One other consideration that eventually came to light is that my test flunks a result that matches its input but the algorithm does nothing to ensure that this does not happen. I didn't notice this with my initial testing as I gave only large arguments to the test case:
permute0TC0_TC0_ &> 1e6$<i.1000
This returns without error most of the time - there's only a one in !1000 chance that I'll return the original as my permuation. However, for smaller arguments - perhaps more germane to my original idea of shuffling a small list of things - the inherent contradiction between the algorithm and the test case shows up:
permute0TC0_TC0_ &>3$<i.10
permute0TC0_TC0_ &>1e6$<i.10
|assertion failure: permute0TC0_TC0_
| y(([:*./e.~)*.([:-.-:)~:]*/ .=[:{.])permute0_base_ y
|permute0TC0_TC0_[0]
The first set of three trials finished without a problem but the larger test on a million instances of ten items uncovers the one in three million chance (since !10 is about three million) that the permutation matches the input, so our assertion fails. Running this latter test several times shows that it fails about one time in three, which is what we'd expect.
There's nothing strictly wrong with my algorithm here - the error is in the test! However, I didn't realize this until much further along in my development of this algorithm, so the appropriate checking doesn't figure into these initial attempts.
In ignorance of this looming problem, I continued with my test case development. Since I allowed for the possibility of multiple test cases per algorithm - in the base name of permute0TC0 - as well as multiple test case suites - in naming the namespace TC0 - I reached for another test case.
This was a reach as it's conceivable - especially in a notation as powerful as J - we could encompass all the necessary criteria in a single test case. However, playing around with permutation vectors and their differences, it occurred to me that the sum of the differences should always be zero.
I don't have a formal proof of this but it seems intuitively obvious. Plus, it held true for a few million experiments and it's trivially obvious for the special cases of when the permutation matches the starting vector or when it equals the "flip" of the original ( |.original ). So, on this shaky foundation, I wrote my next test case:
permute0TC1_TC0_=: 3 : 0 assert. 0=+/y-permute0_base_ y )
As it turned out, this is a bad test case even though it is strictly correct. It's bad because it's a weak assertion, it gets a domain error for non-numeric arguments (which are still allowed at this point in the development of this algorithm) and it has comparison-tolerance issues for floating-point arguments.
Initial Evolution of Test Case Format
As I began work on my second attempt at this algorithm expression, I began to see problems with how I'd set up the test cases. It seems fragile to assume that the function being tested resides in the base namespace. Perhaps the test case should take the original vector and its permutation as arguments? This would push the namespace problem out of the test case, perhaps something like this (using my bad test case as an example):
permute0TC1_TC1_=: 4 : 0 assert. 0=+/x-y )
So, we can run this test case like this:
(] permute0TC1_TC1_ permute0) &> 1e3$<i.100
However, though this is an improvement, my second algorithm attempt revealed a more serious shortcoming with the format of my test cases.
A Simpler Algorithm and More General Test Case
Whereas the first attempt suffered from baroque elaborations, the second one seems almost too simple:
permute1=: 3 : 0
for_ii. i.#y do. rnd=. <.(#y)*?0
y=. (y{~ii,rnd) (rnd,ii)}y end. NB. Swap randomly
)
Here, we just run through the vector, swapping willy-nilly. Once this passes the "eye-ball" test, it's time to get more formal and run our test case - uh oh - our good one has the name of our first attempt hard-wired into it. We could copy and modify this but we see that this does not scale well - we'd have to do this every time and this is patently a bad idea.
Aha! We see that the general solution requires our test cases to be adverbs - the testing code should remain constant as we apply it to different algorithm implementations. This also obviates our dyadic test case idea from before as we can now use any verb within the test, as seen here:
permutes_TC1_=: 1 : 0
assert. y (([: *./ e.~) *. ([: -. -:) ~: ] */ .= [: {. ]) u y
)
permute1 permutes_TC1_ &> 1e3$<i.100
We name the adverb permutes as the plural of permute since it can work with any number of "permute" verbs.
Greater Efficiency Requires Test Cases for Scoring
Now that we have a decent algorithm and test case structure, the next step is to look toward a more efficient algorithm. We need to consider efficiency along two dimensions - run-time and effectiveness of the permuter. This latter criterion implies a new kind of test case - one that rates or scores the performance of an algorithm.
Here's a first cut at scoring a permutation, based on the idea that something like "4 3 0 1 2" is better than "0 1 2 4 3" as a permutation of i.5.
NB.* permutesDiff_TC1_: score how permuted result is from 0 to 1 (best). permutesDiff_TC1_=: 1 : 0 1-(#y)%~+/y=u y NB. 0: no diffs; 1: all different. )
So, applying this to the specific permutations mentioned here ranks the slightly different one as worse than the more different one.
1-5%~+/|0 1 2 4 3 = i.5 0.4 1- 5%~+/|4 3 0 1 2 = i.5 1
Here's an attempt at a more efficient algorithm that works like permute1 but only makes swaps through the vector to its middle. The hope here is that enough of the swaps will randomly hit the second half to permute it as thoroughly as if we had worked through to the end of the vector.
permute2=: 3 : 0
for_ii. i. <.>:-:#y do. rnd=. <.(#y)*?0
y=. (y{~ii,rnd) (rnd,ii)}y end. NB. Swap half randomly
)
The timing comparison is simple enough:
arg=. 1e6$<i.100 6!:2 'permute1 &> arg' 1080.6088 6!:2 'permute2 &> arg' 516.03242
We see that permute2 is about twice as fast as permute1 - at least in J. This makes sense as permute2 does about half the work as permute1.
So, it's faster but how do its results compare? Using this our evaluation test case, we see that permute2 does not do as good a job as permute1:
rr=. (permute1 permutesDiff_TC1_ &> arg);permute2 permutesDiff_TC1_ &> arg=. 1e3$<i.100 usus=: (<./,>./,mean,stddev) NB. Usual Stats: min, max, mean, standard deviation usus &>rr 0.94 1 0.99072 0.010353202 0.61 0.8 0.70125 0.02888203
Note that the mean score for permute2 is only 0.7, compared to about 1.0 for permute1, so our hope that half the work would be sufficient was not fulfilled.
Some Variations on the Algorithm
Further variations on the idea of only iterating through part of the vector leads us to the following variant:
permute3=: 3 : 0
for_ii. i. <.>:-:#y do.
rnd2=. (<:#y) (],-) ii,<.(#y)*?0
y=. (y{~rnd2) (|.rnd2)}y end. NB. Swap 2*half randomly
)
This works through the vector from both ends to the middle, so it does twice the work for half the iterations. How does it compare?
usus TT=: 10+1e4?@$1e3 NB. Mix up the argument sizes 10 1009 504.1635 288.16641 TT=: i.&.>TT 6!:2 'permute3 TT' 19.123175 6!:2 'permute1 TT' 40.703114
As with permute2, permute3 is about twice as fast as permute1 (at least in J - Javascript mileage may differ). How well does this new one work?
sc3=. permute3_base_ permutesDiff_TC1_ &> 1e3$<i.10 NB. Same as scoring argument for permute1 usus sc3 0 1 0.9054 0.1425172
It looks like it's not quite as good but perhaps this is good enough for our purposes.
One file variant, permute4, generalizes the idea from permute3 of doing multiple swaps each time through: in this case, it does four at a time.
permute4=: 3 : 0
bases=. <.sz*4%~i.4 [ sz=. #y
for_ii. i. <.>:4%~sz do.
rnd4=. ~.(sz|bases+ii),<.sz*4?@$0
y=. (y{~rnd4) (|.rnd4)}y end. NB. Swap 4*1/4 randomly
)
Here are some results for these two more-advanced permuters.
load 'C:\amisc\J\NYCJUG\201105\PermuteExperiments.ijs' 6!:2 'sc3=. ;permute3 permutesDiff_TC1_&.>TT' [ TT=: 1e3?@$&.>1e3$1e3 4.8283016 6!:2 'sc4=. ;permute4 permutesDiff_TC1_&.>TT' 2.6690535
So, permute4 is about twice as fast as permute3. How well do these two work?
usus&>sc3;sc4 0.99 1 0.998141 0.0018106345 0.99 1 0.99644 0.0018530607
They work equally well according to the measure we're using, so timing may be the tie-breaker.
Just to check the robustness of our timings, let's time a few variations of arguments.
6!:2 'sc4=. ;permute4 permutesDiff_TC1_&.>TT' [ TT=: 1e3?@$&.>1e3$1e6 2.6797247 6!:2 'sc3=. ;permute3 permutesDiff_TC1_&.>TT' 4.8199215 usus&>sc3;sc4 0.99 1 0.999092 0.0014972568 0.991 1 0.997458 0.0015731213
Run it again in different order (permute3 before permute4) in case there's some kind of caching effect throwing off the relative timings:
6!:2 'sc3=. ;permute3 permutesDiff_TC1_&.>TT' [ TT=: 1e4?@$&.>1e3?@$1e6 47.568377 6!:2 'sc4=. ;permute4 permutesDiff_TC1_&.>TT' 25.919693 usus&>sc3;sc4 0.9983 1 0.9999032 0.00016503635 0.9978 1 0.9997363 0.00018320992
So, permute4 is about twice as fast as permute3 and the two score equally well.
Comment: Uniformly Random Permutations
Note that permute1 doesn't give a uniformly random permutation - e.g. if (3 = #y) then there are 6 possible permutations, but permute1 has 3^3 = 27 possibilities. As 27 isn't a multiple of 6, the permutations can't be equally likely:
b=. permute1"1 ] 10000#,:0 1 2 +/ b-:"1 ] 0 1 2 14861 +/ b-:"1 ] 0 2 1 18534 +/ b-:"1 ] 1 0 2 18220 +/ b-:"1 ] 1 2 0 18668 +/ b-:"1 ] 2 0 1 14957 +/ b-:"1 ] 2 1 0 14760 NB. enumeration shows that above perms occur in ratio 4:5:5:5:4:4
However, if you swap successive elements of y with random elements from that point onwards, then the permutation is uniformly random (by induction):
permute1a=: 3 : 0
for_ii. i.<.#y do. rnd=. ii+<.(ii-~#y)*?0
y=. (y{~ii,rnd) (rnd,ii)}y end. NB. Swap randomly
)
b=. permute1a"1 ] 100000#,:0 1 2
+/"1 (i.6) =/ ((i.6) A. i.3) i. b
16578 16684 16721 16667 16635 16715
-- Ewart Shaw <<DateTime(2011-06-14T10:43:41+0100)>>
A Javascript Implementation
Here's an example of implementing an algorithm developed in J in Javascript. This implements permute1 one of the earlier versions, not the improved one suggested by Ewart Shaw.
First, we assign the URLs from which we'll want to choose randomly.
<script type='text/javascript'> var links=new Array(); links[0]="<li><a href='http://www.aprogramminglanguage.com' target='_blank'>The APL Movie</a></li> "; links[1]="<li><a href='http://www.benbest.com/computer/oopapl.html' target='_blank_'>Ben Best: Object-Oriented Programming and APL</a></li>"; links[2]="<li><a href='http://elliscave.com/APL_J' target='_blank_'>Skip Cave: Matrix Language History Page</a></li>"; links[3]="<li><a href='http://hakank.org/apl/' target='_blank_'>H. Kjellerstrand's APL page</a></li>"; links[4]="<li><a href='http://aplwiki.com/CategoryAplTree' target='_blank_'>APL Tools and Utility Library</a></li>"; links[5]="<li><a href='http://webhome.idirect.com/~kornal' target='_blank_'>Alex Kornilovski</a></li>"; links[6]="<li><a href='http://www.jsoftware.com/papers/mmf.htm' target='_blank_'><i>Memory Mapped Files in J</i> by Chris Burke</a></li>";
We also define the vector of random indexes into "links" outside the function to have global scope.
var ixctr=[];
Note the intention of code as a basic courtesy for those who may follow.
// Randomly permute indexes of "links" so we can choose // group of random links w/o repeats.
The main function:
function disp_links(number)
{ var ubermax=links.length;
We have to loop to do something as simple as "i." .
for (i=0;i<ubermax;i++) {ixctr[i]=i;}
for (i=0;i<ubermax;i++)
{ var k = Math.floor(Math.random()*ubermax);
Random swap requires temporary variable? (I don't know JavaScript very well)
tmp=ixctr[k]; ixctr[k]=ixctr[i]; ixctr[i]=tmp;
}
The function also outputs the random links, doesn't just do the permutation.
for (i=0;i<number;i++)
{ document.writeln(links[ixctr[i]]);
}
}
</script>
Show how to use this.
// Example use: pick 3 links at random. <script type='text/javascript'>disp_links(3);</script></p>
-- Devon McCormick <<DateTime(2013-12-29T19:01:30-0200)>>