00511 {
00512   int i, j;                          /* array indices */
00513   int mn;                            /* number of Walsh differences */
00514   node * head = NULL;                /* points to head of list */
00515   int count;                         /* list print counter */
00516 
00517 
00518   /*----- enter and sort array of differences -----*/
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   /*----- if output requested, write the array of ordered differences -----*/
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   /*----- find median of differences -----*/
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   /*----- deallocate memory -----*/
00547   list_delete (&head);
00548 }

00402 {
00403   const float EPSILON = 1.0e-10;      /* minimum variance limit */
00404   int i, j;                   /* array indices */
00405   node * head = NULL;         /* points to head of list */        
00406   node * ptr = NULL;          /* points to current position in list */
00407   int NN;                     /* total number of sample points */
00408   float wy;                   /* rank sum statistic */
00409   float ewy;                  /* expected value of wy */
00410   float varwy;                /* variance of wy */
00411   int d;                      /* count of number of ties */ 
00412   float corr;                 /* correction to variance to account for ties */
00413   float rank;                 /* rank of data point */
00414 
00415 
00416   /*----- enter and sort x-array -----*/
00417   for (i = 0;  i < m;  i++)
00418     node_addvalue (&head, xarray[i]);
00419 
00420 
00421   /*----- enter and sort y-array -----*/
00422   for (j = 0;  j < n;  j++)
00423     node_addvalue (&head, yarray[j]); 
00424 
00425 
00426   /*----- if display voxel, write the ranks of the input data -----*/
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   /*----- calculate rank sum of y-values -----*/
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   /*----- calculate expected value of rank sum statistic -----*/
00461   NN = m + n;
00462   ewy = 0.5 * n * (NN+1);
00463   if (nvox > 0)  printf ("E(Wy) = %f \n", ewy);
00464 
00465   /*----- calculate variance of rank sum statisitc -----*/
00466   varwy = m * n * (NN+1) / 12.0;
00467 
00468   /*----- correction to variance due to ties in data -----*/
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   /*----- calculate normalized Wilcoxon-Mann-Whitney rank-sum statistic -----*/
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   /*----- deallocate memory -----*/
00491   list_delete (&head);
00492 }

00620 {
00621   int i;                       /* dataset index */
00622   int m;                       /* control sample size */  
00623   int n;                       /* treatment sample size */
00624   int nxyz;                    /* number of voxels per dataset */
00625   int num_datasets;            /* total number of datasets */
00626   int piece_size;              /* number of voxels in dataset sub-volume */
00627   int num_pieces;              /* dataset is divided into this many pieces */
00628   int piece;                   /* piece index */
00629   int piece_len;               /* number of voxels in current sub-volume */
00630   int fim_offset;              /* array offset to current sub-volume */
00631   int ivox;                    /* index to voxels in current sub-volume */
00632   int nvox;                    /* index of voxel within entire volume */
00633   float delta_hat;             /* estimate of shift parameter */
00634   float z;                     /* normalized Mann-Whitney statistic */
00635   float ** xfimar;             /* array of pieces of X-datasets */
00636   float ** yfimar;             /* array of pieces of Y-datasets */
00637   float * xarray;              /* array of control data (X-data) */
00638   float * yarray;              /* array of treatment data (Y-data) */
00639 
00640 
00641   /*----- initialize local variables -----*/
00642   m = option_data->m;
00643   n = option_data->n;
00644   nxyz = option_data->nxyz;
00645   num_datasets = m + n;
00646 
00647 
00648   /*----- break problem into smaller pieces -----*/
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   /*----- allocate memory space -----*/
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   /*----- loop over the pieces of the input datasets -----*/
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       /*----- read in the X-data -----*/
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       /*----- read in the Y-data -----*/
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       /*----- loop over voxels in this piece -----*/
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           /*----- calculate results for this voxel -----*/
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           /*----- save results for this voxel -----*/
00712           delta[ivox+fim_offset] = delta_hat;
00713           zvar[ivox+fim_offset] = z;
00714   
00715         } 
00716           
00717     } /* loop over pieces */
00718 
00719 
00720   /*----- deallocate memory -----*/
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 }

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 }

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 }

00169 {
00170   int nopt = 1;                  /* input option argument counter */
00171   int ival;                      /* integer input */
00172   int nijk;                      /* count of data files */     
00173   float fval;                    /* float input */
00174   THD_3dim_dataset * dset=NULL;             /* test whether data set exists */
00175   char message[MAX_NAME_LENGTH];            /* error message */
00176 
00177 
00178   /*----- does user request help menu? -----*/
00179   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();
00180 
00181   
00182   /*----- add to program log -----*/
00183   AFNI_logger (PROGRAM_NAME,argc,argv); 
00184 
00185 
00186   /*----- initialize the input options -----*/
00187   initialize_options (option_data);
00188 
00189   
00190   /*----- main loop over input options -----*/
00191   while (nopt < argc)
00192     {
00193       
00194       
00195       /*-----   -datum type   -----*/
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 ;  /* go to next arg */
00210       }
00211       
00212       
00213       /*-----   -session dirname    -----*/
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       /*-----   -voxel num  -----*/
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       /*-----   -workmem megabytes  -----*/
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       /*-----   -dset level filename   -----*/
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           /*--- check whether input files exist ---*/
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           /*--- check number of selected sub-bricks ---*/
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       /*-----   -out filename   -----*/
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       /*----- unknown command -----*/
00310       NP_error ("unrecognized command line option ");
00311     }
00312 
00313 }

00350 {
00351   
00352   
00353   /*----- allocate memory space for input data -----*/   
00354   *option_data = (NP_options *) malloc(sizeof(NP_options));
00355   if (*option_data == NULL)
00356     NP_error ("memory allocation error");
00357   
00358   /*----- get command line inputs -----*/
00359   get_options(argc, argv, *option_data);
00360   
00361   /*----- use first data set to get data set dimensions -----*/
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   /*----- check for valid inputs -----*/
00370   check_for_valid_inputs (*option_data);
00371     
00372   /*----- check whether output files already exist -----*/
00373   check_one_output_file (*option_data, (*option_data)->outfile);
00374 
00375   /*----- allocate memory -----*/
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 }

00131 {
00132   int i;          /* index */
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   /*----- allocate memory for storing data file names -----*/
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 }

00818 {
00819   NP_options * option_data = NULL;   /* user input options */
00820   float * delta;                     /* estimated shift parameter */
00821   float * zvar;                      /* normalized Mann-Whitney statistic */
00822  
00823 
00824   /*----- Identify software -----*/
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    /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
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   /*----- program initialization -----*/
00841   initialize (argc, argv, &option_data, &delta, &zvar);
00842   
00843   /*----- calculate nonparameteric Mann-Whitney statistics -----*/
00844   calculate_results (option_data, delta, zvar);
00845   
00846   /*----- generate requested output -----*/
00847   output_results (argc, argv, option_data, delta, zvar);
00848 
00849   /*----- terminate program -----*/
00850   terminate (&option_data, &delta, &zvar);
00851 
00852   exit(0);
00853 }

00753 {
00754 
00755   /*----- write out afni fizt data file -----*/
00756   write_afni_fizt (argc, argv, option_data, option_data->outfile, 
00757                    delta, zvar);
00758 
00759 }

00567 {
00568   int i;                             /* array index */
00569 
00570 
00571   /*----- check for voxel output  -----*/
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   /*----- calculate normalized Mann-Whitney statistic -----*/
00598   calc_stat (nvox, m, n, xarray, yarray, zvar);
00599 
00600 
00601   /*----- estimate shift parameter -----*/
00602   calc_shift (nvox, m, n, xarray, yarray, delta_hat);
00603 
00604 }

00775 {
00776    int i, j;                       /* dataset indices */
00777 
00778 
00779    /*----- deallocate memory -----*/
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 }