00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2001, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*
00008   This program performs the nonparametric Wilcoxon signed-rank test
00009   for paired comparisons of two samples. The output consists of an AFNI 'fizt'
00010   dataset; the first sub-brick contains an estimate of the treatment effect, 
00011   the second sub-brick contains the normalized Wilcoxon signed-rank statistic.
00012 
00013   File:    3dWilcoxon.c
00014   Author:  B. Douglas Ward
00015   Date:    23 July 1997
00016 
00017   Mod:     Added changes for incorporating History notes.
00018   Date:    08 September 1999
00019 
00020   Mod:     Replaced dataset input code with calls to THD_open_dataset,
00021            to allow operator selection of individual sub-bricks for input.
00022   Date:    03 December 1999
00023 
00024   Mod:     Added call to AFNI_logger.
00025   Date:    15 August 2001
00026 
00027   Mod:     Modified routine write_afni_fizt of NPstats.c so that all output 
00028            subbricks will now have the scaled short integer format.
00029   Date:    14 March 2002
00030 
00031   Mod:     Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
00032   Date:    02 December 2002
00033 */
00034 
00035 /*---------------------------------------------------------------------------*/
00036 
00037 #define PROGRAM_NAME "3dWilcoxon"                    /* name of this program */
00038 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
00039 #define PROGRAM_INITIAL "23 July 1997"    /* date of initial program release */
00040 #define PROGRAM_LATEST  "02 Dec  2002"    /* date of latest program revision */
00041 
00042 /*---------------------------------------------------------------------------*/
00043 
00044 #include <stdio.h>
00045 #include <math.h>
00046 #include "mrilib.h"
00047 
00048 #define MAX_OBSERVATIONS 100     /* max. number of observations per cell */
00049 #define MAX_NAME_LENGTH THD_MAX_NAME   /* max. string length for file names */ 
00050 #define MEGA  1048576            /* one megabyte */
00051 
00052 /*---------------------------------------------------------------------------*/
00053 
00054 typedef struct NP_options
00055 { 
00056   int   datum;                  /* data type for "intensity" data subbrick */
00057   char  session[MAX_NAME_LENGTH];     /* name of output directory */
00058 
00059   
00060   int   nvoxel;                 /* number of voxel for special output */
00061 
00062   int   m;                      /* number of X observations */
00063   int   n;                      /* number of Y observations */
00064 
00065   char  *** xname;              /* names of the input data files */
00066   char  * first_dataset;        /* name of the first data set */
00067    
00068   int   nx, ny, nz;             /* data set dimensions */
00069   int   nxyz;                   /* number of voxels per image */
00070 
00071   int workmem;                  /* working memory */
00072 
00073   char  * outfile;              /* name of output file  */
00074 
00075 
00076 } NP_options;
00077 
00078 
00079 #include "NPstats.c"
00080 
00081 
00082 /*---------------------------------------------------------------------------*/
00083 /*
00084    Routine to display 3dWilcoxon help menu.
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    Routine to initialize input options.
00127 */
00128 
00129 void initialize_options (NP_options * option_data)
00130 {
00131   int i;          /* index */
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   /*----- allocate memory for storing data file names -----*/
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    Routine to get user specified Wilcoxon signed-rank test options.
00165 */
00166 
00167 void get_options (int argc, char ** argv, NP_options * option_data)
00168 {
00169   int nopt = 1;                  /* input option argument counter */
00170   int ival;                      /* integer input */
00171   int nijk;                      /* count of data files */     
00172   float fval;                    /* float input */
00173   THD_3dim_dataset * dset=NULL;             /* test whether data set exists */
00174   char message[MAX_NAME_LENGTH];            /* error message */
00175 
00176 
00177   /*----- does user request help menu? -----*/
00178   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();
00179 
00180   
00181   /*----- add to program log -----*/
00182   AFNI_logger (PROGRAM_NAME,argc,argv); 
00183 
00184 
00185   /*----- initialize the input options -----*/
00186   initialize_options (option_data);
00187 
00188   
00189   /*----- main loop over input options -----*/
00190   while (nopt < argc)
00191     {
00192       
00193       
00194       /*-----   -datum type   -----*/
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 ;  /* go to next arg */
00209       }
00210       
00211       
00212       /*-----   -session dirname    -----*/
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       /*-----   -voxel num  -----*/
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       /*-----   -workmem megabytes  -----*/
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       /*-----   -dset level filename   -----*/
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           /*--- check whether input files exist ---*/
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           /*--- check number of selected sub-bricks ---*/
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       /*-----   -out filename   -----*/
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       /*----- unknown command -----*/
00309       NP_error ("unrecognized command line option ");
00310     }
00311 
00312 }
00313 
00314 
00315 /*---------------------------------------------------------------------------*/
00316 /*
00317   Routine to check for valid inputs.
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   Routine to perform all Wilcoxon initialization.
00338 */
00339 
00340 void initialize 
00341 (
00342   int argc,                    /* number of input arguments */
00343   char ** argv,                /* array of input arguments */ 
00344   NP_options ** option_data,   /* user input options */
00345   float ** delta,              /* estimated shift parameter */
00346   float **zvar                 /* normalized Wilcoxon statistic */
00347 )
00348 
00349 {
00350   
00351   
00352   /*----- allocate memory space for input data -----*/   
00353   *option_data = (NP_options *) malloc(sizeof(NP_options));
00354   if (*option_data == NULL)
00355     NP_error ("memory allocation error");
00356   
00357   /*----- get command line inputs -----*/
00358   get_options(argc, argv, *option_data);
00359   
00360   /*----- use first data set to get data set dimensions -----*/
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   /*----- check for valid inputs -----*/
00369   check_for_valid_inputs (*option_data);
00370     
00371   /*----- check whether output files already exist -----*/
00372   check_one_output_file (*option_data, (*option_data)->outfile);
00373 
00374   /*----- allocate memory -----*/
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   Calculate the normalized Wilcoxon signed-rank statistic.
00389 */
00390 
00391 void calc_stat 
00392 (
00393   int nvox,                          /* flag for voxel output */
00394   int n,                             /* sample size */
00395   float * xarray,                    /* array of control data (x data) */
00396   float * yarray,                    /* array of treatment data (y data) */
00397   float * zvar                       /* normalized Wilcoxon statistic */
00398 )
00399 
00400 {
00401   const float EPSILON = 1.0e-10;     /* minimum variance limit */
00402   int i;                      /* array indices */
00403   node * head = NULL;         /* points to head of list */
00404   node * ptr = NULL;          /* points to current position in list */
00405   float wp;                   /* signed-rank statistic */
00406   float ewp;                  /* expected value of wp */
00407   float varwp;                /* variance of wp */
00408   int d;                      /* count of number of ties */
00409   float diff;                 /* difference in pair of observations */
00410 
00411 
00412   /*----- enter and sort absolute values of differences -----*/
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   /*----- calculate sum of ranks of positive differences -----*/
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   /*----- calculate expected value of Wilcoxon signed-rank statistic -----*/
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   /*----- calculate variance of Wilcoxon signed-rank statisitc -----*/
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   /*----- calculate normalized Wilcoxon signed-rank statistic -----*/
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   /*----- deallocate memory -----*/
00488   list_delete (&head);
00489 }
00490 
00491 
00492 /*---------------------------------------------------------------------------*/
00493 /*
00494   Calculate the estimated shift parameter, using the median of the Walsh
00495   averages.
00496 */
00497 
00498 void calc_shift 
00499 (
00500   int nvox,                          /* flag for voxel output */
00501   int n,                             /* sample size */
00502   float * xarray,                    /* array of control data (x data) */
00503   float * yarray,                    /* array of treatment data (y data) */
00504   float * delta_hat                  /* median of differences */
00505 )
00506 
00507 {
00508   int i, j;                           /* array indices */
00509   int mn;                             /* number of Walsh averages */
00510   node * head = NULL;                 /* points to head of list */
00511   float d1, d2;                       /* two paired differences */
00512   int count;                          /* list print counter */
00513 
00514 
00515   /*----- set up ordered array of Walsh averages -----*/
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   /*----- if output requested, write the ordered Walsh averages -----*/
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   /*----- find median of Walsh averages -----*/
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   /*----- deallocate memory -----*/
00549   list_delete (&head);
00550 }
00551 
00552 
00553 /*---------------------------------------------------------------------------*/
00554 /*
00555   Calculate the Wilcoxon statistic and shift parameter for a single voxel.
00556 */
00557 
00558 void process_voxel
00559 (
00560   int nvox,                          /* flag for voxel output */
00561   int n,                             /* sample size */
00562   float * xarray,                    /* array of control data (x data) */
00563   float * yarray,                    /* array of treatment data (y data) */
00564   float * delta_hat,                 /* estimated shift parameter */
00565   float * zvar                       /* normalized Wilcoxon statistic */
00566 )
00567 
00568 {
00569   int i;                             /* array index */
00570 
00571   /*----- check for voxel output  -----*/
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   /*----- calculate normalized Wilcoxon statistic -----*/
00597   calc_stat (nvox, n, xarray, yarray, zvar);
00598 
00599 
00600   /*----- estimate shift parameter -----*/
00601   calc_shift (nvox, n, xarray, yarray, delta_hat);
00602 
00603 }
00604 
00605 
00606 /*---------------------------------------------------------------------------*/
00607 /*
00608   Calculate the Wilcoxon signed-rank statistics for all voxels  (by breaking 
00609   the datasets into sub-volumes, if necessary).
00610 */
00611 
00612 void calculate_results 
00613 (
00614   NP_options * option_data,    /* user input options */
00615   float * delta,               /* estimated shift parameter */
00616   float * zvar                 /* normalized Wilcoxon sign-rank statistic */
00617 )
00618 
00619 {
00620   int i;                       /* dataset index */
00621   int n;                       /* sample size */
00622   int nxyz;                    /* number of voxels per dataset */
00623   int num_datasets;            /* total number of datasets */
00624   int piece_size;              /* number of voxels in dataset sub-volume */
00625   int num_pieces;              /* dataset is divided into this many pieces */
00626   int piece;                   /* piece index */
00627   int piece_len;               /* number of voxels in current sub-volume */
00628   int fim_offset;              /* array offset to current sub-volume */
00629   int ivox;                    /* index to voxels in current sub-volume */
00630   int nvox;                    /* index of voxel within entire volume */
00631   float delta_hat;             /* estimate of shift parameter */
00632   float z;                     /* normalized Wilcoxon rank sum statistic */
00633   float ** xfimar;             /* array of pieces of X-datasets */
00634   float ** yfimar;             /* array of pieces of Y-datasets */
00635   float * xarray;              /* array of control data (X-data) */
00636   float * yarray;              /* array of treatment data (Y-data) */
00637 
00638 
00639   /*----- initialize local variables -----*/
00640   n = option_data->n;
00641   nxyz = option_data->nxyz;
00642   num_datasets = 2 * n;
00643 
00644 
00645   /*----- break problem into smaller pieces -----*/
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   /*----- allocate memory space -----*/
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   /*----- loop over the pieces of the input datasets -----*/
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       /*----- read in the X-data -----*/
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       /*----- read in the Y-data -----*/
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       /*----- loop over voxels in this piece -----*/
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           /*----- calculate results for this voxel -----*/
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           /*----- save results for this voxel -----*/
00710           delta[ivox+fim_offset] = delta_hat;
00711           zvar[ivox+fim_offset] = z;
00712   
00713         } 
00714           
00715     } /* loop over pieces */
00716 
00717 
00718   /*----- deallocate memory -----*/
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   Generate the requested output.
00740 */
00741 
00742 void output_results 
00743 (
00744   int argc,                    /* number of input arguments */
00745   char ** argv,                /* array of input arguments */ 
00746   NP_options * option_data,    /* user input options */
00747   float * delta,               /* estimated shift parameter */
00748   float * zvar                 /* normalized Wilcoxon sign-rank statistic */
00749 )
00750 
00751 {
00752 
00753   /*----- write out afni fizt data file -----*/
00754   write_afni_fizt (argc, argv, option_data, option_data->outfile, 
00755                    delta, zvar);
00756 
00757 }
00758 
00759 
00760 
00761 /*---------------------------------------------------------------------------*/
00762 /*
00763    Routine to release memory and remove any remaining temporary data files.
00764 */
00765 
00766 void terminate (NP_options ** option_data, float ** delta, float **zvar)
00767 {
00768    int i, j;                       /* dataset indices */
00769 
00770 
00771    /*----- deallocate memory -----*/
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    Perform nonparametric Wilcoxon signed-rank paired sample test.
00807 */
00808  
00809 int main (int argc, char ** argv)
00810 {
00811   NP_options * option_data = NULL;   /* user input options */
00812   float * delta;                     /* estimated shift parameter */
00813   float * zvar;                      /* normalized Wilcoxon statistic */
00814   
00815   
00816   /*----- Identify software -----*/
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    /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
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   /*----- program initialization -----*/
00833   initialize (argc, argv, &option_data, &delta, &zvar);
00834   
00835   /*----- calculate nonparameteric Wilcoxon signed-rank statistics -----*/
00836   calculate_results (option_data, delta, zvar);
00837   
00838   /*----- generate requested output -----*/
00839   output_results (argc, argv, option_data, delta, zvar);
00840   
00841   /*----- terminate program -----*/
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