1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2001, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*
8 This program performs the nonparametric Wilcoxon-Mann-Whitney two-sample
9 test for determining whether the observations come from the same population.
10 The output consists of an AFNI 'fizt' dataset; the first sub-brick contains
11 an estimate of the treatment effect, the second sub-brick contains the
12 normalized Wilcoxon rank-sum statistic.
13
14 File: 3dMannWhitney.c
15 Author: B. Douglas Ward
16 Date: 23 July 1997
17
18 Mod: Added changes for incorporating History notes.
19 Date: 08 September 1999
20
21 Mod: Replaced dataset input code with calls to THD_open_dataset,
22 to allow operator selection of individual sub-bricks for input.
23 Date: 03 December 1999
24
25 Mod: Added call to AFNI_logger.
26 Date: 15 August 2001
27
28 Mod: Modified routine write_afni_fizt of NPstats.c so that all output
29 subbricks will now have the scaled short integer format.
30 Date: 14 March 2002
31
32 Mod: Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
33 Date: 02 December 2002
34
35 Mod: Use node_allatonce() to hopefully speed things up.
36 Date: 26 Oct 2007 -- RWCox
37 */
38
39 /*---------------------------------------------------------------------------*/
40
41 #define PROGRAM_NAME "3dMannWhitney" /* name of this program */
42 #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */
43 #define PROGRAM_INITIAL "23 July 1997" /* date of initial program release */
44 #define PROGRAM_LATEST "26 Oct 2007" /* date of latest program revision */
45
46 /*---------------------------------------------------------------------------*/
47
48 #include <stdio.h>
49 #include <math.h>
50 #include "mrilib.h"
51
52 #define MAX_OBSERVATIONS 666 /* max. number of observations per cell */
53 #define MAX_NAME_LENGTH THD_MAX_NAME /* max. string length for file names */
54 #define MEGA 1048576 /* one megabyte */
55
56 /*---------------------------------------------------------------------------*/
57
58 typedef struct NP_options
59 {
60 int datum; /* data type for "intensity" data subbrick */
61 char session[MAX_NAME_LENGTH]; /* name of output directory */
62
63
64 int nvoxel; /* number of voxel for special output */
65
66 int m; /* number of X observations */
67 int n; /* number of Y observations */
68
69 char *** xname; /* names of the input data files */
70 char * first_dataset; /* name of the first data set */
71
72 int nx, ny, nz; /* data set dimensions */
73 int nxyz; /* number of voxels per image */
74
75 int workmem; /* working memory */
76
77 char * outfile; /* name of output file */
78
79
80 } NP_options;
81
82
83 #include "NPstats.c"
84
85 #define USE_ARRAY
86 #ifdef USE_ARRAY
87 static int ntar = 0 ;
88 static float *tar = NULL ;
89 #endif
90
91
92 /*---------------------------------------------------------------------------*/
93 /*
94 Routine to display 3dMannWhitney help menu.
95 */
96
display_help_menu()97 void display_help_menu()
98 {
99 printf
100 (
101 "This program performs nonparametric Mann-Whitney two-sample test. \n\n"
102 "Usage: \n"
103 "3dMannWhitney \n"
104 "-dset 1 filename data set for X observations \n"
105 " . . . . . . \n"
106 "-dset 1 filename data set for X observations \n"
107 "-dset 2 filename data set for Y observations \n"
108 " . . . . . . \n"
109 "-dset 2 filename data set for Y observations \n"
110 " \n"
111 "[-workmem mega] number of megabytes of RAM to use \n"
112 " for statistical workspace \n"
113 "[-voxel num] screen output for voxel # num \n"
114 "-out prefixname estimated population delta and \n"
115 " Wilcoxon-Mann-Whitney statistics \n"
116 " written to file prefixname \n"
117 "\n");
118
119 printf
120 (
121 "\n"
122 "N.B.: For this program, the user must specify 1 and only 1 sub-brick \n"
123 " with each -dset command. That is, if an input dataset contains \n"
124 " more than 1 sub-brick, a sub-brick selector must be used, e.g.: \n"
125 " -dset 2 'fred+orig[3]' \n"
126 );
127
128 printf("\n" MASTER_SHORTHELP_STRING ) ;
129
130 PRINT_COMPILE_DATE ; exit(0);
131 }
132
133
134 /*---------------------------------------------------------------------------*/
135 /*
136 Routine to initialize input options.
137 */
138
initialize_options(NP_options * option_data)139 void initialize_options (NP_options * option_data)
140 {
141 int i; /* index */
142
143 option_data->datum = ILLEGAL_TYPE;
144 strcpy (option_data->session, "./");
145
146
147 option_data->nvoxel = -1;
148
149 option_data->m = 0;
150 option_data->n = 0;
151
152 option_data->workmem = 266;
153
154 /*----- allocate memory for storing data file names -----*/
155 option_data->xname = (char ***) malloc (sizeof(char **) * 2);
156 for (i = 0; i < 2; i++)
157 option_data->xname[i]
158 = (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
159
160 option_data->first_dataset = NULL;
161
162 option_data->nx = 0;
163 option_data->ny = 0;
164 option_data->nz = 0;
165 option_data->nxyz = 0;
166
167 option_data->outfile = NULL;
168
169 }
170
171
172 /*---------------------------------------------------------------------------*/
173 /*
174 Routine to get user specified Mann-Whitney options.
175 */
176
get_options(int argc,char ** argv,NP_options * option_data)177 void get_options (int argc, char ** argv, NP_options * option_data)
178 {
179 int nopt = 1; /* input option argument counter */
180 int ival; /* integer input */
181 int nijk; /* count of data files */
182 float fval; /* float input */
183 THD_3dim_dataset * dset=NULL; /* test whether data set exists */
184 char message[MAX_NAME_LENGTH]; /* error message */
185
186
187 /*----- does user request help menu? -----*/
188 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
189
190
191 /*----- add to program log -----*/
192 AFNI_logger (PROGRAM_NAME,argc,argv);
193
194
195 /*----- initialize the input options -----*/
196 initialize_options (option_data);
197
198
199 /*----- main loop over input options -----*/
200 while (nopt < argc)
201 {
202
203
204 /*----- -datum type -----*/
205 if( strncmp(argv[nopt],"-datum",6) == 0 ){
206 if( ++nopt >= argc ) NP_error("need an argument after -datum!") ;
207
208 if( strcmp(argv[nopt],"short") == 0 ){
209 option_data->datum = MRI_short ;
210 } else if( strcmp(argv[nopt],"float") == 0 ){
211 option_data->datum = MRI_float ;
212 } else {
213 char buf[256] ;
214 sprintf(buf,"-datum of type '%s' is not supported in 3dMannWhitney.",
215 argv[nopt] ) ;
216 NP_error(buf) ;
217 }
218 nopt++ ; continue ; /* go to next arg */
219 }
220
221
222 /*----- -session dirname -----*/
223 if( strncmp(argv[nopt],"-session",6) == 0 ){
224 nopt++ ;
225 if( nopt >= argc ) NP_error("need argument after -session!") ;
226 strcpy(option_data->session , argv[nopt++]) ;
227 continue ;
228 }
229
230
231 /*----- -voxel num -----*/
232 if (strncmp(argv[nopt], "-voxel", 6) == 0)
233 {
234 nopt++;
235 if (nopt >= argc) NP_error ("need argument after -voxel ");
236 sscanf (argv[nopt], "%d", &ival);
237 if (ival <= 0)
238 NP_error ("illegal argument after -voxel ");
239 option_data->nvoxel = ival;
240 nopt++;
241 continue;
242 }
243
244
245 /*----- -workmem megabytes -----*/
246
247 if( strncmp(argv[nopt],"-workmem",6) == 0 ){
248 nopt++ ;
249 if( nopt >= argc ) NP_error ("need argument after -workmem!") ;
250 sscanf (argv[nopt], "%d", &ival);
251 if( ival <= 0 ) NP_error ("illegal argument after -workmem!") ;
252 option_data->workmem = ival ;
253 nopt++ ; continue ;
254 }
255
256
257 /*----- -dset level filename -----*/
258 if (strncmp(argv[nopt], "-dset", 5) == 0)
259 {
260 nopt++;
261 if (nopt+1 >= argc) NP_error ("need 2 arguments after -dset ");
262 sscanf (argv[nopt], "%d", &ival);
263 if ((ival <= 0) || (ival > 2))
264 NP_error ("illegal argument after -dset ");
265
266 if (ival == 1)
267 {
268 option_data->m += 1;
269 nijk = option_data->m;
270 }
271 else
272 {
273 option_data->n += 1;
274 nijk = option_data->n;
275 }
276 if (nijk > MAX_OBSERVATIONS)
277 NP_error ("too many data files");
278
279 /*--- check whether input files exist ---*/
280 nopt++;
281 dset = THD_open_dataset( argv[nopt] ) ;
282 CHECK_OPEN_ERROR(dset,argv[nopt]) ;
283
284 /*--- check number of selected sub-bricks ---*/
285 if (DSET_NVALS(dset) != 1)
286 {
287 sprintf(message,"Must specify exactly 1 sub-brick for file %s\n",
288 argv[nopt]);
289 NP_error (message);
290 }
291
292 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
293
294 option_data->xname[ival-1][nijk-1]
295 = malloc (sizeof(char) * MAX_NAME_LENGTH);
296 strcpy (option_data->xname[ival-1][nijk-1], argv[nopt]);
297 nopt++;
298 continue;
299 }
300
301
302 /*----- -out filename -----*/
303 if (strncmp(argv[nopt], "-out", 4) == 0)
304 {
305 nopt++;
306 if (nopt >= argc) NP_error ("need argument after -out ");
307 option_data->outfile = malloc (sizeof(char) * MAX_NAME_LENGTH);
308 strcpy (option_data->outfile, argv[nopt]);
309 nopt++;
310 continue;
311 }
312
313
314 /*----- unknown command -----*/
315 NP_error ("unrecognized command line option ");
316 }
317
318 if( option_data->m < 2 || option_data->n < 2 )
319 ERROR_exit("Need at least 2 datasets in each of '-dset 1' and '-dset 2'!") ;
320 }
321
322
323 /*---------------------------------------------------------------------------*/
324 /*
325 Routine to check for valid inputs.
326 */
327
check_for_valid_inputs(NP_options * option_data)328 void check_for_valid_inputs (NP_options * option_data)
329 {
330
331 if (option_data->m < 1)
332 NP_error ("too few data sets for X-observations ");
333
334 if (option_data->n < 1)
335 NP_error ("too few data sets for Y-observations ");
336
337 if (option_data->nvoxel > option_data->nxyz)
338 NP_error ("argument of -voxel is too large");
339
340 }
341
342
343 /*---------------------------------------------------------------------------*/
344 /*
345 Routine to perform all Mann-Whitney initialization.
346 */
347
initialize(int argc,char ** argv,NP_options ** option_data,float ** delta,float ** zvar)348 void initialize
349 (
350 int argc, /* number of input arguments */
351 char ** argv, /* array of input arguments */
352 NP_options ** option_data, /* user input options */
353 float ** delta, /* estimated shift parameter */
354 float **zvar /* normalized Mann-Whitney statistic */
355 )
356
357 {
358
359
360 /*----- allocate memory space for input data -----*/
361 *option_data = (NP_options *) malloc(sizeof(NP_options));
362 if (*option_data == NULL)
363 NP_error ("memory allocation error");
364
365 /*----- get command line inputs -----*/
366 get_options(argc, argv, *option_data);
367
368 /*----- use first data set to get data set dimensions -----*/
369 (*option_data)->first_dataset = (*option_data)->xname[0][0];
370 get_dimensions (*option_data);
371 printf ("Data set dimensions: nx = %d ny = %d nz = %d nxyz = %d \n",
372 (*option_data)->nx, (*option_data)->ny,
373 (*option_data)->nz, (*option_data)->nxyz);
374
375
376 /*----- check for valid inputs -----*/
377 check_for_valid_inputs (*option_data);
378
379 /*----- check whether output files already exist -----*/
380 if( THD_deathcon() ) check_one_output_file (*option_data, (*option_data)->outfile);
381
382 /*----- allocate memory -----*/
383 *delta = (float *) malloc(sizeof(float) * (*option_data)->nxyz);
384 if (*delta == NULL)
385 NP_error ("memory allocation error");
386 *zvar = (float *) malloc(sizeof(float) * (*option_data)->nxyz);
387 if (*zvar == NULL)
388 NP_error ("memory allocation error");
389
390
391 }
392
393
394 /*---------------------------------------------------------------------------*/
395 /*
396 Calculate the normalized Wilcoxon-Mann-Whitney statistic.
397 */
398
calc_stat(int nvox,int m,int n,float * xarray,float * yarray,float * zvar)399 void calc_stat
400 (
401 int nvox, /* flag for voxel output */
402 int m, /* number of control subjects */
403 int n, /* number of treatment subjects */
404 float * xarray, /* array of control data (x data) */
405 float * yarray, /* array of treatment data (y data) */
406 float * zvar /* normalized Mann-Whitney statistic */
407 )
408
409 {
410 const float EPSILON = 1.0e-10; /* minimum variance limit */
411 int i, j; /* array indices */
412 node * head = NULL; /* points to head of list */
413 node * ptr = NULL; /* points to current position in list */
414 int NN; /* total number of sample points */
415 float wy; /* rank sum statistic */
416 float ewy; /* expected value of wy */
417 float varwy; /* variance of wy */
418 int d; /* count of number of ties */
419 float corr; /* correction to variance to account for ties */
420 float rank; /* rank of data point */
421
422 #ifdef USE_ARRAY
423 if( ntar == 0 ){
424 d = m+n+1 ; ntar = m*n+1 ; if( ntar < d ) ntar = d ;
425 tar = (float *)malloc(sizeof(float)*ntar) ;
426 if( tar == NULL ) ERROR_exit("Can't allocate workspace 'tar[%d]'",ntar) ;
427 }
428 #endif
429
430
431 #ifdef USE_ARRAY
432 memcpy( tar , xarray, sizeof(float)*m ) ;
433 memcpy( tar+m , yarray , sizeof(float)*n ) ;
434 node_allatonce( &head , n+m , tar ) ;
435 #else
436 /*----- enter and sort x-array -----*/
437 for (i = 0; i < m; i++) node_addvalue (&head, xarray[i]);
438
439 /*----- enter and sort y-array -----*/
440 for (j = 0; j < n; j++) node_addvalue (&head, yarray[j]);
441 #endif
442
443
444 /*----- if display voxel, write the ranks of the input data -----*/
445 if (nvox > 0)
446 {
447 printf ("\n");
448 printf ("X ranks: \n");
449 for (i = 0; i < m; i++)
450 {
451 rank = node_get_rank (head, xarray[i]);
452 printf (" %6.1f", rank);
453 if (((i+1) % 10 == 0) && (i < m-1))
454 printf ("\n");
455 }
456
457 printf ("\n\n");
458 printf ("Y ranks: \n");
459 for (i = 0; i < n; i++)
460 {
461 rank = node_get_rank (head, yarray[i]);
462 printf (" %6.1f", rank);
463 if (((i+1) % 10 == 0) && (i < n-1))
464 printf ("\n");
465 }
466
467 printf ("\n\n");
468 }
469
470
471 /*----- calculate rank sum of y-values -----*/
472 wy = 0.0;
473 for (j = 0; j < n; j++)
474 wy += node_get_rank (head, yarray[j]);
475 if (nvox > 0) printf ("Wy = %f \n", wy);
476
477
478 /*----- calculate expected value of rank sum statistic -----*/
479 NN = m + n;
480 ewy = 0.5 * n * (NN+1);
481 if (nvox > 0) printf ("E(Wy) = %f \n", ewy);
482
483 /*----- calculate variance of rank sum statisitc -----*/
484 varwy = m * n * (NN+1) / 12.0;
485
486 /*----- correction to variance due to ties in data -----*/
487 corr = 0.0;
488 ptr = head;
489 while (ptr != NULL)
490 {
491 d = ptr->d;
492 corr += d*d*d - d;
493 ptr = ptr->next;
494 }
495 corr *= m * n / (12.0 * NN * (NN - 1));
496 varwy -= corr;
497 if (nvox > 0) printf ("Var(Wy) = %f \n", varwy);
498
499
500 /*----- calculate normalized Wilcoxon-Mann-Whitney rank-sum statistic -----*/
501 if (varwy < EPSILON)
502 *zvar = 0.0;
503 else
504 *zvar = (wy - ewy) / sqrt(varwy);
505 if (nvox > 0) printf ("Z = %f \n", *zvar);
506
507
508 /*----- deallocate memory -----*/
509 list_delete (&head);
510 }
511
512
513 /*---------------------------------------------------------------------------*/
514 /*
515 Calculate the estimated shift parameter, using the median of the Walsh
516 differences.
517 */
518
calc_shift(int nvox,int m,int n,float * xarray,float * yarray,float * delta_hat)519 void calc_shift
520 (
521 int nvox, /* flag for voxel output */
522 int m, /* number of control subjects */
523 int n, /* number of treatment subjects */
524 float * xarray, /* array of control data (x data) */
525 float * yarray, /* array of treatment data (y data) */
526 float * delta_hat /* median of differences */
527 )
528
529 {
530 register int i, j; /* array indices */
531 int mn; /* number of Walsh differences */
532 node * head = NULL; /* points to head of list */
533 int count; /* list print counter */
534
535
536 /*----- enter and sort array of differences -----*/
537 #ifdef USE_ADDARRAY
538 { register float xi ; register int k ;
539 for( k=i=0 ; i < m ; i++ ){
540 xi = xarray[i] ;
541 for( j=0 ; j < n ; j++ ) tar[k++] = yarray[j] - xi ;
542 }
543 node_allatonce( &head , m*n , tar ) ;
544 }
545 #else
546 for (i = 0; i < m; i++)
547 for (j = 0; j < n; j++)
548 node_addvalue (&head, yarray[j] - xarray[i]);
549 #endif
550
551
552 /*----- if output requested, write the array of ordered differences -----*/
553 if (nvox > 0)
554 {
555 printf ("\n");
556 printf ("Ordered differences: \n");
557 count = 0;
558 list_print (head, &count);
559 printf ("\n");
560 }
561
562
563 /*----- find median of differences -----*/
564 mn = m * n;
565 *delta_hat = node_get_median (head, mn);
566
567 if (nvox > 0)
568 {
569 printf ("\n");
570 printf ("Delta hat = %f \n\n", *delta_hat);
571 }
572
573
574 /*----- deallocate memory -----*/
575 list_delete (&head);
576 }
577
578
579 /*---------------------------------------------------------------------------*/
580 /*
581 Calculate the Mann-Whitney statistic and shift parameter for a single voxel.
582 */
583
process_voxel(int nvox,int m,int n,float * xarray,float * yarray,float * delta_hat,float * zvar)584 void process_voxel
585 (
586 int nvox, /* flag for voxel output */
587 int m, /* number of control subjects */
588 int n, /* number of treatment subjects */
589 float * xarray, /* array of control data (x data) */
590 float * yarray, /* array of treatment data (y data) */
591 float * delta_hat, /* estimated shift parameter */
592 float * zvar /* normalized Mann-Whitney statistic */
593 )
594
595 {
596 int i; /* array index */
597
598
599 /*----- check for voxel output -----*/
600 if (nvox > 0)
601 {
602 printf ("\nResults for voxel #%d : \n", nvox);
603
604 printf ("\n");
605 printf ("X data: \n");
606 for (i = 0; i < m; i++)
607 {
608 printf (" %6.1f", xarray[i]);
609 if (((i+1) % 10 == 0) && (i < m-1))
610 printf ("\n");
611 }
612
613 printf ("\n\n");
614 printf ("Y data: \n");
615 for (i = 0; i < n; i++)
616 {
617 printf (" %6.1f", yarray[i]);
618 if (((i+1) % 10 == 0) && (i < n-1))
619 printf ("\n");
620 }
621 printf ("\n");
622 }
623
624 /*----- calculate normalized Mann-Whitney statistic -----*/
625 calc_stat (nvox, m, n, xarray, yarray, zvar);
626
627
628 /*----- estimate shift parameter -----*/
629 calc_shift (nvox, m, n, xarray, yarray, delta_hat);
630
631 }
632
633
634 /*---------------------------------------------------------------------------*/
635 /*
636 Calculate the Wilcoxon-Mann-Whitney rank-sum statistics for all voxels
637 (by breaking the datasets into sub-volumes, if necessary).
638 */
639
calculate_results(NP_options * option_data,float * delta,float * zvar)640 void calculate_results
641 (
642 NP_options * option_data, /* user input options */
643 float * delta, /* estimated shift parameter */
644 float * zvar /* normalized Mann-Whitney statistic */
645 )
646
647 {
648 int i; /* dataset index */
649 int m; /* control sample size */
650 int n; /* treatment sample size */
651 int nxyz; /* number of voxels per dataset */
652 int num_datasets; /* total number of datasets */
653 int piece_size; /* number of voxels in dataset sub-volume */
654 int num_pieces; /* dataset is divided into this many pieces */
655 int piece; /* piece index */
656 int piece_len; /* number of voxels in current sub-volume */
657 int fim_offset; /* array offset to current sub-volume */
658 int ivox; /* index to voxels in current sub-volume */
659 int nvox; /* index of voxel within entire volume */
660 float delta_hat; /* estimate of shift parameter */
661 float z; /* normalized Mann-Whitney statistic */
662 float ** xfimar; /* array of pieces of X-datasets */
663 float ** yfimar; /* array of pieces of Y-datasets */
664 float * xarray; /* array of control data (X-data) */
665 float * yarray; /* array of treatment data (Y-data) */
666
667
668 /*----- initialize local variables -----*/
669 m = option_data->m;
670 n = option_data->n;
671 nxyz = option_data->nxyz;
672 num_datasets = m + n;
673
674
675 /*----- break problem into smaller pieces -----*/
676 piece_size = option_data->workmem * MEGA / (num_datasets * sizeof(float));
677 if (piece_size > nxyz) piece_size = nxyz;
678 num_pieces = (nxyz + piece_size - 1) / piece_size;
679 printf ("num_pieces = %d piece_size = %d \n", num_pieces, piece_size);
680
681 /*----- allocate memory space -----*/
682 xarray = (float *) malloc (sizeof(float) * m); MTEST(xarray);
683 yarray = (float *) malloc (sizeof(float) * n); MTEST(yarray);
684 xfimar = (float **) malloc (sizeof(float *) * m); MTEST(xfimar);
685 yfimar = (float **) malloc (sizeof(float *) * n); MTEST(yfimar);
686 for (i = 0; i < m; i++)
687 {
688 xfimar[i] = (float *) malloc(sizeof(float) * piece_size);
689 MTEST(xfimar[i]);
690 }
691 for (i = 0; i < n; i++)
692 {
693 yfimar[i] = (float *) malloc(sizeof(float) * piece_size);
694 MTEST(yfimar[i]);
695 }
696
697
698 /*----- loop over the pieces of the input datasets -----*/
699 nvox = 0;
700 for (piece = 0; piece < num_pieces; piece++)
701 {
702 printf ("piece = %d \n", piece);
703 fim_offset = piece * piece_size;
704 if (piece < num_pieces-1)
705 piece_len = piece_size;
706 else
707 piece_len = nxyz - fim_offset;
708
709 /*----- read in the X-data -----*/
710 for (i = 0; i < m; i++)
711 read_afni_data (option_data, option_data->xname[0][i],
712 piece_len, fim_offset, xfimar[i]);
713
714 /*----- read in the Y-data -----*/
715 for (i = 0; i < n; i++)
716 read_afni_data (option_data, option_data->xname[1][i],
717 piece_len, fim_offset, yfimar[i]);
718
719
720 /*----- loop over voxels in this piece -----*/
721 for (ivox = 0; ivox < piece_len; ivox++)
722 {
723 nvox += 1;
724
725 for (i = 0; i < m; i++)
726 xarray[i] = xfimar[i][ivox];
727 for (i = 0; i < n; i++)
728 yarray[i] = yfimar[i][ivox];
729
730
731 /*----- calculate results for this voxel -----*/
732 if (nvox == option_data->nvoxel)
733 process_voxel (nvox, m, n, xarray, yarray, &delta_hat, &z);
734 else
735 process_voxel (-1, m, n, xarray, yarray, &delta_hat, &z);
736
737
738 /*----- save results for this voxel -----*/
739 delta[ivox+fim_offset] = delta_hat;
740 zvar[ivox+fim_offset] = z;
741
742 }
743
744 } /* loop over pieces */
745
746
747 /*----- deallocate memory -----*/
748 free (xarray); xarray = NULL;
749 free (yarray); yarray = NULL;
750
751 for (i = 0; i < m; i++)
752 {
753 free (xfimar[i]); xfimar[i] = NULL;
754 }
755 free (xfimar); xfimar = NULL;
756
757 for (i = 0; i < n; i++)
758 {
759 free (yfimar[i]); yfimar[i] = NULL;
760 }
761 free (yfimar); yfimar = NULL;
762
763 }
764
765
766 /*---------------------------------------------------------------------------*/
767 /*
768 Generate the requested output.
769 */
770
output_results(int argc,char ** argv,NP_options * option_data,float * delta,float * zvar)771 void output_results
772 (
773 int argc, /* number of input arguments */
774 char ** argv, /* array of input arguments */
775 NP_options * option_data, /* user input options */
776 float * delta, /* estimated shift parameter */
777 float * zvar /* normalized Mann-Whitney rank-sum statistic */
778 )
779
780 {
781
782 /*----- write out afni fizt data file -----*/
783 write_afni_fizt (argc, argv, option_data, option_data->outfile,
784 delta, zvar);
785
786 }
787
788
789
790 /*---------------------------------------------------------------------------*/
791 /*
792 Routine to release memory and remove any remaining temporary data files.
793 */
794
terminate(NP_options ** option_data,float ** delta,float ** zvar)795 void terminate
796 (
797 NP_options ** option_data, /* user input options */
798 float ** delta, /* estimated shift parameter */
799 float ** zvar /* normalized Mann-Whitney rank-sum statistic */
800 )
801
802 {
803 int i, j; /* dataset indices */
804
805
806 /*----- deallocate memory -----*/
807 for (j = 0; j < (*option_data)->m; j++)
808 {
809 free ((*option_data)->xname[0][j]);
810 (*option_data)->xname[0][j] = NULL;
811 }
812 for (j = 0; j < (*option_data)->n; j++)
813 {
814 free ((*option_data)->xname[1][j]);
815 (*option_data)->xname[1][j] = NULL;
816 }
817 for (i = 0; i < 2; i++)
818 {
819 free ((*option_data)->xname[i]);
820 (*option_data)->xname[i] = NULL;
821 }
822 free ((*option_data)->xname);
823 (*option_data)->xname = NULL;
824
825 if ((*option_data)->outfile != NULL)
826 {
827 free ((*option_data)-> outfile);
828 (*option_data)->outfile = NULL;
829 }
830
831 free (*option_data); *option_data = NULL;
832
833 free (*delta); *delta = NULL;
834
835 free (*zvar); *zvar = NULL;
836 }
837
838
839 /*---------------------------------------------------------------------------*/
840 /*
841 Perform nonparametric Wilcoxon-Mann-Whitney rank-sum two sample test.
842 */
843
main(int argc,char ** argv)844 int main (int argc, char ** argv)
845 {
846 NP_options * option_data = NULL; /* user input options */
847 float * delta; /* estimated shift parameter */
848 float * zvar; /* normalized Mann-Whitney statistic */
849
850
851 /*----- Identify software -----*/
852 #if 0
853 printf ("\n\n");
854 printf ("Program: %s \n", PROGRAM_NAME);
855 printf ("Author: %s \n", PROGRAM_AUTHOR);
856 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
857 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
858 printf ("\n");
859 #endif
860
861 /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
862
863 PRINT_VERSION("3dMannWhitney") ; AUTHOR(PROGRAM_AUTHOR);
864 mainENTRY("3dMannWhitney main") ; machdep() ;
865
866 { int new_argc ; char ** new_argv ;
867 addto_args( argc , argv , &new_argc , &new_argv ) ;
868 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
869 }
870
871 /*----- program initialization -----*/
872 initialize (argc, argv, &option_data, &delta, &zvar);
873
874 /*----- calculate nonparameteric Mann-Whitney statistics -----*/
875 calculate_results (option_data, delta, zvar);
876
877 /*----- generate requested output -----*/
878 output_results (argc, argv, option_data, delta, zvar);
879
880 /*----- terminate program -----*/
881 terminate (&option_data, &delta, &zvar);
882
883 exit(0);
884 }
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902