;+ ; NAME: ; PERMUTE ; ; PURPOSE: ; This function returns an array containing the numbers ; [0, ..., N-1] in random order. They are useful as indices ; when permuting a dataset, for example in a balanced bootstrap ; Monte Carlo algorithm. ; ; CATEGORY: ; Statistics. ; ; CALLING SEQUENCE: ; ; Result = PERMUTE(N) ; ; INPUTS: ; N: The number of items to be permuted. ; ; OPTIONAL INPUTS: ; SEED: A random number seed, see RANDOMU. ; ; OUTPUTS: ; This function returns an N-element array containing a random ; permutation of the integers from 0 through N-1. ; ; SIDE EFFECTS: ; Unless Seed is specified, IDL's global random number ; seed is changed. ; ; PROCEDURE: ; This is an in-place swapping algorithm. It starts with an ; index array. For each position in the array, it swaps the ; occupant of that position with the occupant of a random ; position from there (inclusive) to the end of the array. The ; last iteration is not necessary to compute, since it swaps ; with itself. ; ; See http://www.techuser.net/randpermgen.html for a proof. The ; 2-line code there has been optimized for IDL's vector ; architecture. This is a linear-time algorithm. ; ; EXAMPLE: ; Show some permutations of 6 numbers: ; print, permute(6) ; 0 2 1 3 4 5 ; print, permute(6) ; 2 4 3 5 1 0 ; print, permute(6) ; 0 4 3 1 2 5 ; ; Permute the array [2, 4, 6, 8] ; a = [2, 4, 6, 8] ; print, a[permute(4)] ; 4 8 6 2 ; ; Test randomness (results should be close to k): ; m = 6l ; k = 10000l ; n = m * k ; a = lonarr(m, n) ; for i = 0l, n-1, 1 do a[*, i] = permute(m) ; for i = 0l, m-1, 1 do print, histogram(a[i, *]) ; 9885 10062 10051 9915 10028 10059 ; 10096 10087 10094 9913 9933 9877 ; 10041 10013 9968 9958 9911 10109 ; 9880 9858 10166 10049 10081 9966 ; 10093 9915 9800 10166 9969 10057 ; 10005 10065 9921 9999 10078 9932 ; ; Time the algorithm: ; maxn = 7 ; t = dblarr(maxn) ; n = 10L^(indgen(maxn)+1) ; for i = 0, maxn-1, 1 do begin &\$ ; t1 = systime(/s) &\$ ; print, n[i] &\$ ; a = permute(n[i]) &\$ ; t2 = systime(/s) &\$ ; t[i] = t2-t1 &\$ ; endfor ; print, ' Elements Seconds Elements Per Second' ; print, transpose([[n], [t], [t/n]]) ; ; Elements Seconds Elements Per Second ; 10.000000 0.00012397766 1.2397766e-05 ; 100.00000 0.00015020370 1.5020370e-06 ; 1000.0000 0.0011651516 1.1651516e-06 ; 10000.000 0.018178225 1.8178225e-06 ; 100000.00 0.13504505 1.3504505e-06 ; 1000000.0 1.3817160 1.3817160e-06 ; 10000000. 14.609985 1.4609985e-06 ; ; These times are for a 2.071 GHz AMD Athlon 2800+ CPU. ; ; MODIFICATION HISTORY: ; Written by: Joseph Harrington, Cornell. 2006-03-22 ; jh@alum.mit.edu ;- function PERMUTE, N, Seed ; Don't stop here! on_error, 2 ; test inputs if n eq 1 then return, 0L if n lt 1 then message, 'N = ' + strtrim(n, 2) + ', must be 1 or more.' ar = lindgen(n) rar = reverse(ar[0 : n - 2]) + 2 r = (n - 1) - long( randomu(seed, n - 1) * rar ) for i = 0L, n - 2, 1 do begin t = ar[i] ar[i] = ar[r[i]] ar[r[i]] = t endfor return, ar end