1 /*
2 * See copyright in copyright.h and the accompanying file COPYING
3 */
4
5 /*
6 *========================================================================
7 * This just counts the permutations of n samples. They should
8 * occur n! times each. We count them and do a straight chisq.
9 *========================================================================
10 */
11
12 #include <dieharder/libdieharder.h>
13
14 #define RGB_PERM_KMAX 10
15 uint nperms;
16 double fpipi(int pi1,int pi2,int nkp);
17 uint rgb_permutations_k;
18
rgb_permutations(Test ** test,int irun)19 int rgb_permutations(Test **test,int irun)
20 {
21
22 uint i,k,permindex=0,t;
23 Vtest vtest;
24 double *testv;
25 size_t ps[4096];
26 gsl_permutation** lookup;
27
28
29 MYDEBUG(D_RGB_PERMUTATIONS){
30 printf("#==================================================================\n");
31 printf("# rgb_permutations: Debug with %u\n",D_RGB_PERMUTATIONS);
32 }
33
34 /*
35 * Number of permutations. Note that the minimum ntuple value for a
36 * valid test is 2. If ntuple is less than 2, we choose the default
37 * test size as 5 (like operm5).
38 */
39 if(ntuple<2){
40 test[0]->ntuple = 5;
41 } else {
42 test[0]->ntuple = ntuple;
43 }
44 k = test[0]->ntuple;
45 nperms = gsl_sf_fact(k);
46
47 /*
48 * A vector to accumulate rands in some sort order
49 */
50 testv = (double *)malloc(k*sizeof(double));
51
52 MYDEBUG(D_RGB_PERMUTATIONS){
53 printf("# rgb_permutations: There are %u permutations of length k = %u\n",nperms,k);
54 }
55
56 /*
57 * Create a test, initialize it.
58 */
59 Vtest_create(&vtest,nperms);
60 vtest.cutoff = 5.0;
61 for(i=0;i<nperms;i++){
62 vtest.x[i] = 0.0;
63 vtest.y[i] = (double) test[0]->tsamples/nperms;
64 }
65
66 MYDEBUG(D_RGB_PERMUTATIONS){
67 printf("# rgb_permutations: Allocating permutation lookup table.\n");
68 }
69 lookup = (gsl_permutation**) malloc(nperms*sizeof(gsl_permutation*));
70 for(i=0;i<nperms;i++){
71 lookup[i] = gsl_permutation_alloc(k);
72 }
73 for(i=0;i<nperms;i++){
74 if(i == 0){
75 gsl_permutation_init(lookup[i]);
76 } else {
77 gsl_permutation_memcpy(lookup[i],lookup[i-1]);
78 gsl_permutation_next(lookup[i]);
79 }
80 }
81
82 MYDEBUG(D_RGB_PERMUTATIONS){
83 for(i=0;i<nperms;i++){
84 printf("# rgb_permutations: %u => ",i);
85 gsl_permutation_fprintf(stdout,lookup[i]," %u");
86 printf("\n");
87 }
88 }
89
90 /*
91 * We count the order permutations in a long string of samples of
92 * rgb_permutation_k non-overlapping rands. This is done by:
93 * a) Filling testv[] with rgb_permutation_k rands.
94 * b) Using gsl_sort_index to generate the permutation index.
95 * c) Incrementing a counter for that index (a-c done tsamples times)
96 * d) Doing a straight chisq on the counter vector with nperms-1 DOF
97 *
98 * This test should be done with tsamples > 30*nperms, easily met for
99 * reasonable rgb_permutation_k
100 */
101 for(t=0;t<test[0]->tsamples;t++){
102 /*
103 * To sort into a perm, test vector needs to be double.
104 */
105 for(i=0;i<k;i++) {
106 testv[i] = (double) gsl_rng_get(rng);
107 MYDEBUG(D_RGB_PERMUTATIONS){
108 printf("# rgb_permutations: testv[%u] = %u\n",i,(uint) testv[i]);
109 }
110 }
111
112 gsl_sort_index(ps,testv,1,k);
113
114 MYDEBUG(D_RGB_PERMUTATIONS){
115 for(i=0;i<k;i++) {
116 printf("# rgb_permutations: ps[%u] = %lu\n",i,ps[i]);
117 }
118 }
119
120 for(i=0;i<nperms;i++){
121 if(memcmp(ps,lookup[i]->data,k*sizeof(size_t))==0){
122 permindex = i;
123 MYDEBUG(D_RGB_PERMUTATIONS){
124 printf("# Found permutation: ");
125 gsl_permutation_fprintf(stdout,lookup[i]," %u");
126 printf(" = %u\n",i);
127 }
128 break;
129 }
130 }
131
132 vtest.x[permindex]++;
133 MYDEBUG(D_RGB_PERMUTATIONS){
134 printf("# rgb_permutations: Augmenting vtest.x[%u] = %f\n",permindex,vtest.x[permindex]);
135 }
136
137 }
138
139 MYDEBUG(D_RGB_PERMUTATIONS){
140 printf("# rgb_permutations:==============================\n");
141 printf("# rgb_permutations: permutation count = \n");
142 for(i=0;i<nperms;i++){
143 printf("# count[%u] = %u\n",i,(uint) vtest.x[i]);
144 }
145 }
146
147 Vtest_eval(&vtest);
148 test[0]->pvalues[irun] = vtest.pvalue;
149 MYDEBUG(D_RGB_PERMUTATIONS) {
150 printf("# rgb_permutations(): test[0]->pvalues[%u] = %10.5f\n",irun,test[0]->pvalues[irun]);
151 }
152
153 for(i=0;i<nperms;i++){
154 gsl_permutation_free(lookup[i]);
155 }
156 free(lookup);
157 free(testv);
158 Vtest_destroy(&vtest);
159
160 return(0);
161
162 }
163
164