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