00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #define PROGRAM_NAME "3dWilcoxon"
00038 #define PROGRAM_AUTHOR "B. Douglas Ward"
00039 #define PROGRAM_INITIAL "23 July 1997"
00040 #define PROGRAM_LATEST "02 Dec 2002"
00041
00042
00043
00044 #include <stdio.h>
00045 #include <math.h>
00046 #include "mrilib.h"
00047
00048 #define MAX_OBSERVATIONS 100
00049 #define MAX_NAME_LENGTH THD_MAX_NAME
00050 #define MEGA 1048576
00051
00052
00053
00054 typedef struct NP_options
00055 {
00056 int datum;
00057 char session[MAX_NAME_LENGTH];
00058
00059
00060 int nvoxel;
00061
00062 int m;
00063 int n;
00064
00065 char *** xname;
00066 char * first_dataset;
00067
00068 int nx, ny, nz;
00069 int nxyz;
00070
00071 int workmem;
00072
00073 char * outfile;
00074
00075
00076 } NP_options;
00077
00078
00079 #include "NPstats.c"
00080
00081
00082
00083
00084
00085
00086
00087 void display_help_menu()
00088 {
00089 printf
00090 (
00091 "This program performs the nonparametric Wilcoxon signed-rank test \n"
00092 "for paired comparisons of two samples. \n\n"
00093 "Usage: \n"
00094 "3dWilcoxon \n"
00095 "-dset 1 filename data set for X observations \n"
00096 " . . . . . . \n"
00097 "-dset 1 filename data set for X observations \n"
00098 "-dset 2 filename data set for Y observations \n"
00099 " . . . . . . \n"
00100 "-dset 2 filename data set for Y observations \n"
00101 " \n"
00102 "[-workmem mega] number of megabytes of RAM to use \n"
00103 " for statistical workspace \n"
00104 "[-voxel num] screen output for voxel # num \n"
00105 "-out prefixname estimated population delta and \n"
00106 " Wilcoxon signed-rank statistics are\n"
00107 " written to file prefixname \n"
00108 "\n");
00109
00110 printf
00111 (
00112 "\n"
00113 "N.B.: For this program, the user must specify 1 and only 1 sub-brick \n"
00114 " with each -dset command. That is, if an input dataset contains \n"
00115 " more than 1 sub-brick, a sub-brick selector must be used, e.g.: \n"
00116 " -dset 2 'fred+orig[3]' \n"
00117 );
00118
00119 printf("\n" MASTER_SHORTHELP_STRING ) ;
00120 exit(0);
00121 }
00122
00123
00124
00125
00126
00127
00128
00129 void initialize_options (NP_options * option_data)
00130 {
00131 int i;
00132
00133 option_data->datum = ILLEGAL_TYPE;
00134 strcpy (option_data->session, "./");
00135
00136
00137 option_data->nvoxel = -1;
00138
00139 option_data->m = 0;
00140 option_data->n = 0;
00141
00142 option_data->workmem = 12;
00143
00144
00145 option_data->xname = (char ***) malloc (sizeof(char **) * 2);
00146 for (i = 0; i < 2; i++)
00147 option_data->xname[i]
00148 = (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
00149
00150 option_data->first_dataset = NULL;
00151
00152 option_data->nx = 0;
00153 option_data->ny = 0;
00154 option_data->nz = 0;
00155 option_data->nxyz = 0;
00156
00157 option_data->outfile = NULL;
00158
00159 }
00160
00161
00162
00163
00164
00165
00166
00167 void get_options (int argc, char ** argv, NP_options * option_data)
00168 {
00169 int nopt = 1;
00170 int ival;
00171 int nijk;
00172 float fval;
00173 THD_3dim_dataset * dset=NULL;
00174 char message[MAX_NAME_LENGTH];
00175
00176
00177
00178 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
00179
00180
00181
00182 AFNI_logger (PROGRAM_NAME,argc,argv);
00183
00184
00185
00186 initialize_options (option_data);
00187
00188
00189
00190 while (nopt < argc)
00191 {
00192
00193
00194
00195 if( strncmp(argv[nopt],"-datum",6) == 0 ){
00196 if( ++nopt >= argc ) NP_error("need an argument after -datum!") ;
00197
00198 if( strcmp(argv[nopt],"short") == 0 ){
00199 option_data->datum = MRI_short ;
00200 } else if( strcmp(argv[nopt],"float") == 0 ){
00201 option_data->datum = MRI_float ;
00202 } else {
00203 char buf[256] ;
00204 sprintf(buf,"-datum of type '%s' is not supported in 3dMannWhitney.",
00205 argv[nopt] ) ;
00206 NP_error(buf) ;
00207 }
00208 nopt++ ; continue ;
00209 }
00210
00211
00212
00213 if( strncmp(argv[nopt],"-session",6) == 0 ){
00214 nopt++ ;
00215 if( nopt >= argc ) NP_error("need argument after -session!") ;
00216 strcpy(option_data->session , argv[nopt++]) ;
00217 continue ;
00218 }
00219
00220
00221
00222 if (strncmp(argv[nopt], "-voxel", 6) == 0)
00223 {
00224 nopt++;
00225 if (nopt >= argc) NP_error ("need argument after -voxel ");
00226 sscanf (argv[nopt], "%d", &ival);
00227 if (ival <= 0)
00228 NP_error ("illegal argument after -voxel ");
00229 option_data->nvoxel = ival;
00230 nopt++;
00231 continue;
00232 }
00233
00234
00235
00236
00237 if( strncmp(argv[nopt],"-workmem",6) == 0 ){
00238 nopt++ ;
00239 if( nopt >= argc ) NP_error ("need argument after -workmem!") ;
00240 sscanf (argv[nopt], "%d", &ival);
00241 if( ival <= 0 ) NP_error ("illegal argument after -workmem!") ;
00242 option_data->workmem = ival ;
00243 nopt++ ; continue ;
00244 }
00245
00246
00247
00248 if (strncmp(argv[nopt], "-dset", 5) == 0)
00249 {
00250 nopt++;
00251 if (nopt+1 >= argc) NP_error ("need 2 arguments after -dset ");
00252 sscanf (argv[nopt], "%d", &ival);
00253 if ((ival <= 0) || (ival > 2))
00254 NP_error ("illegal argument after -dset ");
00255
00256 if (ival == 1)
00257 {
00258 option_data->m += 1;
00259 nijk = option_data->m;
00260 }
00261 else
00262 {
00263 option_data->n += 1;
00264 nijk = option_data->n;
00265 }
00266 if (nijk > MAX_OBSERVATIONS)
00267 NP_error ("too many data files");
00268
00269
00270 nopt++;
00271 dset = THD_open_dataset( argv[nopt] ) ;
00272 if( ! ISVALID_3DIM_DATASET(dset) )
00273 {
00274 sprintf(message,"Unable to open dataset file %s\n", argv[nopt]);
00275 NP_error (message);
00276 }
00277
00278
00279 if (DSET_NVALS(dset) != 1)
00280 {
00281 sprintf(message,"Must specify exactly 1 sub-brick for file %s\n",
00282 argv[nopt]);
00283 NP_error (message);
00284 }
00285
00286 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00287
00288 option_data->xname[ival-1][nijk-1]
00289 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00290 strcpy (option_data->xname[ival-1][nijk-1], argv[nopt]);
00291 nopt++;
00292 continue;
00293 }
00294
00295
00296
00297 if (strncmp(argv[nopt], "-out", 4) == 0)
00298 {
00299 nopt++;
00300 if (nopt >= argc) NP_error ("need argument after -out ");
00301 option_data->outfile = malloc (sizeof(char) * MAX_NAME_LENGTH);
00302 strcpy (option_data->outfile, argv[nopt]);
00303 nopt++;
00304 continue;
00305 }
00306
00307
00308
00309 NP_error ("unrecognized command line option ");
00310 }
00311
00312 }
00313
00314
00315
00316
00317
00318
00319
00320 void check_for_valid_inputs (NP_options * option_data)
00321 {
00322
00323 if (option_data->m != option_data->n)
00324 NP_error ("Must have equal sample sizes for paired comparisons!");
00325
00326 if (option_data->n < 1)
00327 NP_error ("too few data sets ");
00328
00329 if (option_data->nvoxel > option_data->nxyz)
00330 NP_error ("argument of -voxel is too large");
00331
00332 }
00333
00334
00335
00336
00337
00338
00339
00340 void initialize
00341 (
00342 int argc,
00343 char ** argv,
00344 NP_options ** option_data,
00345 float ** delta,
00346 float **zvar
00347 )
00348
00349 {
00350
00351
00352
00353 *option_data = (NP_options *) malloc(sizeof(NP_options));
00354 if (*option_data == NULL)
00355 NP_error ("memory allocation error");
00356
00357
00358 get_options(argc, argv, *option_data);
00359
00360
00361 (*option_data)->first_dataset = (*option_data)->xname[0][0];
00362 get_dimensions (*option_data);
00363 printf ("Data set dimensions: nx = %d ny = %d nz = %d nxyz = %d \n",
00364 (*option_data)->nx, (*option_data)->ny,
00365 (*option_data)->nz, (*option_data)->nxyz);
00366
00367
00368
00369 check_for_valid_inputs (*option_data);
00370
00371
00372 check_one_output_file (*option_data, (*option_data)->outfile);
00373
00374
00375 *delta = (float *) malloc(sizeof(float) * (*option_data)->nxyz);
00376 if (*delta == NULL)
00377 NP_error ("memory allocation error");
00378 *zvar = (float *) malloc(sizeof(float) * (*option_data)->nxyz);
00379 if (*zvar == NULL)
00380 NP_error ("memory allocation error");
00381
00382
00383 }
00384
00385
00386
00387
00388
00389
00390
00391 void calc_stat
00392 (
00393 int nvox,
00394 int n,
00395 float * xarray,
00396 float * yarray,
00397 float * zvar
00398 )
00399
00400 {
00401 const float EPSILON = 1.0e-10;
00402 int i;
00403 node * head = NULL;
00404 node * ptr = NULL;
00405 float wp;
00406 float ewp;
00407 float varwp;
00408 int d;
00409 float diff;
00410
00411
00412
00413 if (nvox > 0) printf ("\n\nY - X: \n");
00414 for (i = 0; i < n; i++)
00415 {
00416 diff = yarray[i] - xarray[i];
00417 if (nvox > 0)
00418 {
00419 printf (" %6.1f", diff);
00420 if (((i+1) % 10 == 0) && (i < n-1))
00421 printf ("\n");
00422 }
00423 node_addvalue (&head, fabs(diff));
00424 }
00425
00426
00427
00428 if (nvox > 0) printf ("\n\nSigned Ranks: \n");
00429 wp = 0.0;
00430 for (i = 0; i < n; i++)
00431 {
00432 diff = yarray[i] - xarray[i];
00433 if (diff > 0.0)
00434 wp += node_get_rank (head, diff);
00435 if (nvox > 0)
00436 {
00437 if (diff > 0.0)
00438 printf (" %6.1f", node_get_rank(head, diff));
00439 else if (diff < 0.0)
00440 printf (" %6.1f", -node_get_rank(head, -diff));
00441 else
00442 printf (" %6.1f", 0.0);
00443 if (((i+1) % 10 == 0) && (i < n-1))
00444 printf ("\n");
00445 }
00446 }
00447 if (nvox > 0) printf ("\n\nW+ = %f \n", wp);
00448
00449
00450
00451 ewp = n * (n+1) / 4.0;
00452 ptr = head;
00453 if (ptr->fval == 0.0)
00454 {
00455 d = ptr->d;
00456 ewp -= d * (d+1) / 4.0;
00457 }
00458 if (nvox > 0) printf ("E(W+) = %f \n", ewp);
00459
00460
00461
00462 varwp = n * (n+1) * (2*n+1) / 24.0;
00463 ptr = head;
00464 if (ptr->fval == 0.0)
00465 {
00466 d = ptr->d;
00467 varwp -= d * (d+1) * (2*d+1) / 24.0;
00468 ptr = ptr->next;
00469 }
00470 while (ptr != NULL)
00471 {
00472 d = ptr->d;
00473 varwp -= d * (d-1) * (d+1) / 48.0;
00474 ptr = ptr->next;
00475 }
00476 if (nvox > 0) printf ("Var(W+) = %f \n", varwp);
00477
00478
00479
00480 if (varwp < EPSILON)
00481 *zvar = 0.0;
00482 else
00483 *zvar = (wp - ewp) / sqrt(varwp);
00484 if (nvox > 0) printf ("Z = %f \n", *zvar);
00485
00486
00487
00488 list_delete (&head);
00489 }
00490
00491
00492
00493
00494
00495
00496
00497
00498 void calc_shift
00499 (
00500 int nvox,
00501 int n,
00502 float * xarray,
00503 float * yarray,
00504 float * delta_hat
00505 )
00506
00507 {
00508 int i, j;
00509 int mn;
00510 node * head = NULL;
00511 float d1, d2;
00512 int count;
00513
00514
00515
00516 for (i = 0; i < n; i++)
00517 {
00518 d1 = yarray[i] - xarray[i];
00519 for (j = i; j < n; j++)
00520 {
00521 d2 = yarray[j] - xarray[j];
00522 node_addvalue (&head, (d1 + d2) / 2.0);
00523 }
00524 }
00525
00526
00527 if (nvox > 0)
00528 {
00529 printf ("\n");
00530 printf ("Ordered Walsh averages: \n");
00531 count = 0;
00532 list_print (head, &count);
00533 printf ("\n");
00534 }
00535
00536
00537
00538 mn = n*(n+1)/2;
00539 *delta_hat = node_get_median (head, mn);
00540
00541 if (nvox > 0)
00542 {
00543 printf ("\n");
00544 printf ("Delta hat = %f \n\n", *delta_hat);
00545 }
00546
00547
00548
00549 list_delete (&head);
00550 }
00551
00552
00553
00554
00555
00556
00557
00558 void process_voxel
00559 (
00560 int nvox,
00561 int n,
00562 float * xarray,
00563 float * yarray,
00564 float * delta_hat,
00565 float * zvar
00566 )
00567
00568 {
00569 int i;
00570
00571
00572 if (nvox > 0)
00573 {
00574 printf ("\n\nResults for voxel #%d : \n", nvox);
00575
00576 printf ("\n");
00577 printf ("X data: \n");
00578 for (i = 0; i < n; i++)
00579 {
00580 printf (" %6.1f", xarray[i]);
00581 if (((i+1) % 10 == 0) && (i < n-1))
00582 printf ("\n");
00583 }
00584
00585 printf ("\n\n");
00586 printf ("Y data: \n");
00587 for (i = 0; i < n; i++)
00588 {
00589 printf (" %6.1f", yarray[i]);
00590 if (((i+1) % 10 == 0) && (i < n-1))
00591 printf ("\n");
00592 }
00593 }
00594
00595
00596
00597 calc_stat (nvox, n, xarray, yarray, zvar);
00598
00599
00600
00601 calc_shift (nvox, n, xarray, yarray, delta_hat);
00602
00603 }
00604
00605
00606
00607
00608
00609
00610
00611
00612 void calculate_results
00613 (
00614 NP_options * option_data,
00615 float * delta,
00616 float * zvar
00617 )
00618
00619 {
00620 int i;
00621 int n;
00622 int nxyz;
00623 int num_datasets;
00624 int piece_size;
00625 int num_pieces;
00626 int piece;
00627 int piece_len;
00628 int fim_offset;
00629 int ivox;
00630 int nvox;
00631 float delta_hat;
00632 float z;
00633 float ** xfimar;
00634 float ** yfimar;
00635 float * xarray;
00636 float * yarray;
00637
00638
00639
00640 n = option_data->n;
00641 nxyz = option_data->nxyz;
00642 num_datasets = 2 * n;
00643
00644
00645
00646 piece_size = option_data->workmem * MEGA / (num_datasets * sizeof(float));
00647 if (piece_size > nxyz) piece_size = nxyz;
00648 num_pieces = (nxyz + piece_size - 1) / piece_size;
00649 printf ("num_pieces = %d piece_size = %d \n", num_pieces, piece_size);
00650
00651
00652
00653 xarray = (float *) malloc (sizeof(float) * n); MTEST(xarray);
00654 yarray = (float *) malloc (sizeof(float) * n); MTEST(yarray);
00655 xfimar = (float **) malloc (sizeof(float *) * n); MTEST(xfimar);
00656 yfimar = (float **) malloc (sizeof(float *) * n); MTEST(yfimar);
00657 for (i = 0; i < n; i++)
00658 {
00659 xfimar[i] = (float *) malloc(sizeof(float) * piece_size);
00660 MTEST(xfimar[i]);
00661 }
00662 for (i = 0; i < n; i++)
00663 {
00664 yfimar[i] = (float *) malloc(sizeof(float) * piece_size);
00665 MTEST(yfimar[i]);
00666 }
00667
00668
00669
00670 nvox = 0;
00671 for (piece = 0; piece < num_pieces; piece++)
00672 {
00673 printf ("piece = %d \n", piece);
00674 fim_offset = piece * piece_size;
00675 if (piece < num_pieces-1)
00676 piece_len = piece_size;
00677 else
00678 piece_len = nxyz - fim_offset;
00679
00680
00681 for (i = 0; i < n; i++)
00682 read_afni_data (option_data, option_data->xname[0][i],
00683 piece_len, fim_offset, xfimar[i]);
00684
00685
00686 for (i = 0; i < n; i++)
00687 read_afni_data (option_data, option_data->xname[1][i],
00688 piece_len, fim_offset, yfimar[i]);
00689
00690
00691
00692 for (ivox = 0; ivox < piece_len; ivox++)
00693 {
00694 nvox += 1;
00695
00696 for (i = 0; i < n; i++)
00697 xarray[i] = xfimar[i][ivox];
00698 for (i = 0; i < n; i++)
00699 yarray[i] = yfimar[i][ivox];
00700
00701
00702
00703 if (nvox == option_data->nvoxel)
00704 process_voxel (nvox, n, xarray, yarray, &delta_hat, &z);
00705 else
00706 process_voxel (-1, n, xarray, yarray, &delta_hat, &z);
00707
00708
00709
00710 delta[ivox+fim_offset] = delta_hat;
00711 zvar[ivox+fim_offset] = z;
00712
00713 }
00714
00715 }
00716
00717
00718
00719 free (xarray); xarray = NULL;
00720 free (yarray); yarray = NULL;
00721
00722 for (i = 0; i < n; i++)
00723 {
00724 free (xfimar[i]); xfimar[i] = NULL;
00725 }
00726 free (xfimar); xfimar = NULL;
00727
00728 for (i = 0; i < n; i++)
00729 {
00730 free (yfimar[i]); yfimar[i] = NULL;
00731 }
00732 free (yfimar); yfimar = NULL;
00733
00734 }
00735
00736
00737
00738
00739
00740
00741
00742 void output_results
00743 (
00744 int argc,
00745 char ** argv,
00746 NP_options * option_data,
00747 float * delta,
00748 float * zvar
00749 )
00750
00751 {
00752
00753
00754 write_afni_fizt (argc, argv, option_data, option_data->outfile,
00755 delta, zvar);
00756
00757 }
00758
00759
00760
00761
00762
00763
00764
00765
00766 void terminate (NP_options ** option_data, float ** delta, float **zvar)
00767 {
00768 int i, j;
00769
00770
00771
00772 for (j = 0; j < (*option_data)->n; j++)
00773 {
00774 free ((*option_data)->xname[0][j]);
00775 (*option_data)->xname[0][j] = NULL;
00776 }
00777 for (j = 0; j < (*option_data)->n; j++)
00778 {
00779 free ((*option_data)->xname[1][j]);
00780 (*option_data)->xname[1][j] = NULL;
00781 }
00782 for (i = 0; i < 2; i++)
00783 {
00784 free ((*option_data)->xname[i]);
00785 (*option_data)->xname[i] = NULL;
00786 }
00787 free ((*option_data)->xname);
00788 (*option_data)->xname = NULL;
00789
00790 if ((*option_data)->outfile != NULL)
00791 {
00792 free ((*option_data)-> outfile);
00793 (*option_data)->outfile = NULL;
00794 }
00795
00796 free (*option_data); *option_data = NULL;
00797
00798 free (*delta); *delta = NULL;
00799
00800 free (*zvar); *zvar = NULL;
00801 }
00802
00803
00804
00805
00806
00807
00808
00809 int main (int argc, char ** argv)
00810 {
00811 NP_options * option_data = NULL;
00812 float * delta;
00813 float * zvar;
00814
00815
00816
00817 printf ("\n\n");
00818 printf ("Program: %s \n", PROGRAM_NAME);
00819 printf ("Author: %s \n", PROGRAM_AUTHOR);
00820 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
00821 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
00822 printf ("\n");
00823
00824
00825
00826 machdep() ;
00827 { int new_argc ; char ** new_argv ;
00828 addto_args( argc , argv , &new_argc , &new_argv ) ;
00829 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00830 }
00831
00832
00833 initialize (argc, argv, &option_data, &delta, &zvar);
00834
00835
00836 calculate_results (option_data, delta, zvar);
00837
00838
00839 output_results (argc, argv, option_data, delta, zvar);
00840
00841
00842 terminate (&option_data, &delta, &zvar);
00843
00844 exit(0);
00845 }
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862