1 #include "mrilib.h"
2
make_two_random_subsets(int bot,int top,int n1,int n2)3 int ** make_two_random_subsets( int bot , int top , int n1 , int n2 )
4 {
5 int num=top-bot+1 , jj,nn , bad ;
6 int *ar1 , *ar2 , *all , **arr ;
7
8 if( num < 2 || n1 < 1 || n2 < 1 || n1+n2 > num ) return NULL ;
9
10 arr = (int **)malloc(sizeof(int *)*2);
11 ar1 = (int * )malloc(sizeof(int)*n1) ; arr[0] = ar1 ;
12 ar2 = (int * )malloc(sizeof(int)*n2) ; arr[1] = ar2 ;
13 all = (int * )malloc(sizeof(int)*num);
14 for( jj=0 ; jj < num ; jj++ ) all[jj] = bot+jj ;
15
16 bad = bot - 666 ;
17 for( nn=0 ; nn < n1 ; nn++ ){
18 while(1){
19 jj = lrand48() % num ;
20 if( all[jj] >= bot ){
21 ar1[nn] = all[jj] ; all[jj] = bad ; break ;
22 }
23 }
24 }
25
26 if( n1+n2 == num ){ /* all the rest into ar2 */
27 for( nn=jj=0 ; jj < num ; jj++ ){
28 if( all[jj] >= bot ) ar2[nn++] = all[jj] ;
29 }
30 } else { /* random subset of what's left */
31 for( nn=0 ; nn < n2 ; nn++ ){
32 while(1){
33 jj = lrand48() % num ;
34 if( all[jj] >= bot ){
35 ar2[nn] = all[jj] ; all[jj] = bad ; break ;
36 }
37 }
38 }
39 }
40
41 free(all) ; qsort_int(n1,ar1) ; qsort_int(n2,ar2) ; return arr ;
42 }
43
44 /*----------------------------------------------------------------------------*/
45
main(int argc,char * argv[])46 int main( int argc , char *argv[] )
47 {
48 int jj , bot , top , n1 , n2 , *ar1 , *ar2 , **arr , iarg , docomma=0 ;
49 char *prefix = "AFNIroolz" , *fnam ; FILE *fp ;
50
51 if( argc < 3 || strcasecmp(argv[1],"-help") == 0 ){
52 printf(
53 "Usage: 2perm [-prefix PPP] [-comma] bot top [n1 n2]\n"
54 "\n"
55 "This program creates 2 random non-overlapping subsets of the set of\n"
56 "integers from 'bot' to 'top' (inclusive). The first subset is of\n"
57 "length 'n1' and the second of length 'n2'. If those values are not\n"
58 "given, then equal size subsets of length (top-bot+1)/2 are used.\n"
59 "\n"
60 "This program is intended for use in various simulation and/or\n"
61 "randomization scripts, or for amusement/hilarity.\n"
62 "\n"
63 "OPTIONS:\n"
64 "========\n"
65 " -prefix PPP == Two output files are created, with names PPP_A and PPP_B,\n"
66 " where 'PPP' is the given prefix. If no '-prefix' option\n"
67 " is given, then the string 'AFNIroolz' will be used.\n"
68 " ++ Each file is a single column of numbers.\n"
69 " ++ Note that the filenames do NOT end in '.1D'.\n"
70 "\n"
71 " -comma == Write each file as a single row of comma-separated numbers.\n"
72 "\n"
73 "EXAMPLE:\n"
74 "========\n"
75 "This illustration shows the purpose of 2perm -- for use in permutation\n"
76 "and/or randomization tests of statistical significance and power.\n"
77 "Given a dataset with 100 sub-bricks (indexed 0..99), split it into two\n"
78 "random halves and do a 2-sample t-test between them.\n"
79 "\n"
80 " 2perm -prefix Q50 0 99\n"
81 " 3dttest++ -setA dataset+orig\"[1dcat Q50_A]\" \\\n"
82 " -setB dataset+orig\"[1dcat Q50_B]\" \\\n"
83 " -no1sam -prefix Q50\n"
84 " \\rm -f Q50_?\n"
85 "\n"
86 "Alternatively:\n"
87 "\n"
88 " 2perm -prefix Q50 -comma 0 99\n"
89 " 3dttest++ -setA dataset+orig\"[`cat Q50_A`]\" \\\n"
90 " -setB dataset+orig\"[`cat Q50_B`]\" \\\n"
91 " -no1sam -prefix Q50\n"
92 " \\rm -f Q50_?\n"
93 "\n"
94 "Note the combined use of the double quote \" and backward quote `\n"
95 "shell operators in this second approach.\n"
96 "\n"
97 "AUTHOR: (no one want to admit they wrote this trivial code).\n"
98 ) ;
99 PRINT_COMPILE_DATE ; exit(0) ;
100 }
101
102 machdep() ;
103
104 iarg = 1 ;
105 while( iarg < argc && argv[iarg][0] == '-' ){
106
107 if( strcasecmp(argv[iarg],"-prefix") == 0 ){
108 if( ++iarg >= argc )
109 ERROR_exit("2perm: need an argument after option %s",argv[iarg-1]) ;
110 prefix = argv[iarg] ;
111 if( !THD_filename_ok(prefix) ) ERROR_exit("2perm: illegal prefix") ;
112 iarg++ ; continue ;
113 }
114
115 if( strcasecmp(argv[iarg],"-comma") == 0 ){
116 docomma = 1 ; iarg++ ; continue ;
117 }
118
119 ERROR_exit("Unknown option '%s'",argv[iarg]) ;
120
121 }
122
123 fnam = (char *)malloc(sizeof(char)*(strlen(prefix)+8)) ;
124
125 if( iarg+2 > argc )
126 ERROR_exit("2perm: need at least 2 arguments for 'bot' and 'top'") ;
127
128 bot = (int)strtod(argv[iarg++],NULL) ;
129 top = (int)strtod(argv[iarg++],NULL) ;
130 if( bot >= top )
131 ERROR_exit("2perm: bot=%d is not less than top=%d",bot,top) ;
132
133 if( iarg+2 <= argc ){
134 n1 = (int)strtod(argv[iarg++],NULL) ;
135 n2 = (int)strtod(argv[iarg++],NULL) ;
136 if( n1 < 1 || n2 < 1 )
137 ERROR_exit("2perm: n1=%d and/or n2=%d are not both positive",n1,n2) ;
138 if( n1+n2 > top-bot+1 )
139 ERROR_exit("2perm: n1=%d + n2=%d is bigger than number of points from %d to %d",
140 n1,n2 , bot,top ) ;
141 } else {
142 n2 = n1 = (top-bot+1)/2 ;
143 }
144
145 arr = make_two_random_subsets(bot,top,n1,n2) ;
146 if( arr == NULL ) ERROR_exit("2perm: bad inputs -- can't create subsets") ;
147 ar1 = arr[0] ; ar2 = arr[1] ;
148
149 sprintf(fnam,"%s_A",prefix) ;
150 fp = fopen(fnam,"w") ;
151 if( fp == NULL ) ERROR_exit("2perm: Can't open file %s for output!",fnam) ;
152 if( docomma ){
153 for( jj=0 ; jj < n1-1 ; jj++ ) fprintf(fp,"%d,",ar1[jj]) ;
154 fprintf(fp,"%d\n",ar1[n1-1]) ;
155 } else {
156 for( jj=0 ; jj < n1 ; jj++ ) fprintf(fp,"%d\n",ar1[jj]) ;
157 }
158 fclose(fp) ;
159
160 sprintf(fnam,"%s_B",prefix) ;
161 fp = fopen(fnam,"w") ;
162 if( fp == NULL ) ERROR_exit("2perm: Can't open file %s for output!",fnam) ;
163 if( docomma ){
164 for( jj=0 ; jj < n2-1 ; jj++ ) fprintf(fp,"%d,",ar2[jj]) ;
165 fprintf(fp,"%d\n",ar2[n2-1]) ;
166 } else {
167 for( jj=0 ; jj < n2 ; jj++ ) fprintf(fp,"%d\n",ar2[jj]) ;
168 }
169 fclose(fp) ;
170
171 exit(0) ;
172 }
173