00420 {
00421   const float EPSILON = 1.0e-10;      /* protection from roundoff error */
00422   int i, j;                   /* array indices */
00423   node * head = NULL;         /* points to head of list */        
00424   node * ptr = NULL;          /* points to current position in list */
00425   int NN;                     /* total number of sample points */
00426   float rsum;                 /* sum of squares of ranks */
00427   int d;                      /* count of number of ties */ 
00428   float corr;                 /* correction to account for ties */
00429   float rank;                 /* rank of data point */
00430   float ranksum;              /* sum of ranks for ith treatment */
00431   float knum;                 /* numerator of Kruskal-Wallis statistic */
00432   float kden;                 /* denominator of Kruskal-Wallis statistic */
00433   float best_rank;            /* best average rank for a treatment */
00434 
00435 
00436   /*----- count total number of datasets -----*/
00437   NN = 0;
00438   for (i = 0;  i < s;  i++)
00439     NN += n[i];
00440 
00441 
00442   /*----- enter data arrays -----*/
00443   for (i = 0;  i < s;  i++)
00444     for (j = 0;  j < n[i];  j++)
00445       node_addvalue (&head, xarray[i][j]);
00446 
00447 
00448   /*----- if display voxel, write the ranks of the input data -----*/
00449   if (nvox > 0)
00450     {
00451       printf ("\n");
00452       for (i = 0;  i < s;  i++)
00453         {
00454           printf ("Y%d ranks: ", i+1);
00455           for (j = 0;  j < n[i];  j++)
00456             {
00457               rank = node_get_rank (head, xarray[i][j]);
00458               printf (" %6.1f", rank);
00459               if (((j+1) % 8 == 0) && (j < n[i]-1))
00460                 printf ("\n          ");                  
00461             }
00462           printf ("\n");
00463         }
00464       printf ("\n");
00465       for (i = 0;  i < s;  i++)
00466         {
00467           printf ("Y%d: ", i+1);
00468           ranksum = 0.0;
00469           for (j = 0;  j < n[i];  j++)
00470             {
00471               rank = node_get_rank (head, xarray[i][j]);
00472               ranksum += rank;
00473             }
00474           printf ("   Rank sum = %6.1f    Rank average = %6.1f \n", 
00475                   ranksum, ranksum/n[i]);
00476         }
00477       printf ("\n");
00478     }
00479 
00480 
00481   /*----- calculate the sum of the rank averages -----*/
00482   rsum = 0.0;
00483   *best = 0.0;
00484   best_rank = (NN + 1.0) / 2.0 + EPSILON;
00485   for (i = 0;  i < s;  i++)
00486     {
00487       ranksum = 0.0;
00488       for (j = 0;  j < n[i];  j++)
00489         ranksum += node_get_rank (head, xarray[i][j]);
00490       rsum += ranksum * ranksum / n[i];
00491 
00492       if (ranksum / n[i] > best_rank)
00493         {
00494           *best = (float) (i+1);
00495           best_rank = ranksum / n[i];
00496         }
00497     }
00498 
00499 
00500   /*----- numerator of Kruskal-Wallis statistic -----*/
00501   knum = (12.0/(NN*(NN+1)))*rsum - 3.0*(NN+1);
00502 
00503 
00504   /*----- denominator of Kruskal-Wallis statistic -----*/
00505   corr = 0.0;
00506   ptr = head;
00507   while (ptr != NULL)
00508     {
00509       d = ptr->d;
00510       corr += d*d*d - d;
00511       ptr = ptr->next;
00512     }
00513   kden = 1.0 - (corr / (NN*NN*NN-NN));
00514 
00515 
00516   /*----- calculate Kruskal-Wallis statistic -----*/
00517   if (kden < EPSILON)
00518     *kstat = 0.0;
00519   else
00520     *kstat = knum / kden;
00521   if (nvox > 0)  printf ("K = %f \n", *kstat);
00522 
00523 
00524   /*----- deallocate memory -----*/
00525   list_delete (&head);
00526 }

00589 {
00590   int i, j, m;                 /* array indices */
00591   int s;                       /* number of treatments */  
00592   int * n;                     /* number of observations per treatment */
00593   int nxyz;                    /* number of voxels per dataset */
00594   int num_datasets;            /* total number of datasets */
00595   int piece_size;              /* number of voxels in dataset sub-volume */
00596   int num_pieces;              /* dataset is divided into this many pieces */
00597   int piece;                   /* piece index */
00598   int piece_len;               /* number of voxels in current sub-volume */
00599   int fim_offset;              /* array offset to current sub-volume */
00600   int ivox;                    /* index to voxels in current sub-volume */
00601   int nvox;                    /* index of voxel within entire volume */
00602   float b;                     /* index of best treatment */
00603   float k;                     /* Kruskal-Wallis statistic */
00604   float ** xfimar;             /* array of sub-volumes of datasets */
00605   float ** xarray;             /* array of data arrays */
00606 
00607 
00608   /*----- initialize local variables -----*/
00609   s = option_data->s;
00610   nxyz = option_data->nxyz;
00611   num_datasets = 0;
00612   n = (int *) malloc (sizeof(int) * s);  MTEST(n);
00613   for (i = 0;  i < s;  i++)
00614     {
00615       n[i] = option_data->n[i];
00616       num_datasets += n[i];
00617     }
00618 
00619 
00620   /*----- break problem into smaller pieces -----*/
00621   piece_size = option_data->workmem * MEGA / (num_datasets * sizeof(float));
00622   if (piece_size > nxyz)  piece_size = nxyz;
00623   num_pieces = (nxyz + piece_size - 1) / piece_size;
00624   printf ("num_pieces = %d    piece_size = %d \n", num_pieces, piece_size);    
00625 
00626   
00627   /*----- allocate memory space -----*/
00628   xarray = (float **) malloc (sizeof(float *) * s);  MTEST(xarray);
00629   for (i = 0;  i < s;  i++)
00630     {
00631       xarray[i] = (float *) malloc (sizeof(float) * option_data->n[i]);
00632       MTEST(xarray[i]);
00633     }
00634 
00635   xfimar = (float **) malloc (sizeof(float *) * num_datasets);  MTEST(xfimar);
00636   for (i = 0;  i < num_datasets;  i++)
00637     {
00638       xfimar[i] = (float *) malloc (sizeof(float) * piece_size);  
00639       MTEST(xfimar[i]);
00640     }
00641 
00642 
00643   /*----- loop over the pieces of the input datasets -----*/
00644   nvox = 0;
00645   for (piece = 0;  piece < num_pieces;  piece++)
00646     {
00647       printf ("piece = %d \n", piece);
00648       fim_offset = piece * piece_size;
00649       if (piece < num_pieces-1)
00650         piece_len = piece_size;
00651       else
00652         piece_len = nxyz - fim_offset;
00653 
00654 
00655       /*----- read in sub-volume of data from each dataset -----*/
00656       m = 0;
00657       for (i = 0;  i < s;  i++)
00658         for (j = 0;  j < option_data->n[i];  j++)
00659           {
00660             read_afni_data (option_data, option_data->xname[i][j],
00661                             piece_len, fim_offset, xfimar[m]);
00662             m++;
00663           }
00664 
00665 
00666       /*----- loop over voxels in this piece -----*/
00667       for (ivox = 0;  ivox < piece_len;  ivox++)
00668         {
00669           nvox += 1;
00670 
00671           m = 0;
00672           for (i = 0;  i < s;  i++)
00673             for (j = 0;  j < option_data->n[i];  j++)
00674               {
00675                 xarray[i][j] = xfimar[m][ivox];
00676                 m++;
00677               }
00678 
00679 
00680           /*----- calculate results for this voxel -----*/
00681           if (nvox == option_data->nvoxel)
00682             process_voxel (nvox, s, n, xarray, &b, &k);
00683           else
00684             process_voxel (-1, s, n, xarray, &b, &k);
00685     
00686 
00687           /*----- save results for this voxel -----*/
00688           best[ivox+fim_offset] = b;
00689           kstat[ivox+fim_offset] = k;
00690   
00691         } 
00692           
00693     }  /* loop over pieces */
00694 
00695 
00696   /*----- deallocate memory -----*/
00697   free (n);   n = NULL;
00698   
00699   for (i = 0;  i < s;  i++)
00700     {
00701       free (xarray[i]);   xarray[i] = NULL;
00702     }
00703   free (xarray);   xarray = NULL;
00704 
00705   for (i = 0;  i < num_datasets;  i++)
00706     {
00707       free (xfimar[i]);   xfimar[i] = NULL;
00708     }
00709   free (xfimar);   xfimar = NULL;
00710 }

00335 {
00336   int i;
00337   char message[MAX_NAME_LENGTH];            /* error message */
00338 
00339 
00340   for (i = 0;  i < option_data->s;  i++)
00341     if (option_data->n[i] < 1)
00342       {
00343         sprintf(message,"Too few data sets for treatment level %d \n", 
00344                 i+1);
00345         NP_error (message);
00346       } 
00347 
00348   if (option_data->nvoxel > option_data->nxyz)
00349     NP_error ("argument of -voxel is too large");
00350 
00351 }

00089 {
00090   printf 
00091     (
00092      "This program performs nonparametric Kruskal-Wallis test for         \n"
00093      "comparison of multiple treatments.                                \n\n"
00094      "Usage:                                                              \n"
00095      "3dKruskalWallis                                                     \n"
00096      "-levels s                      s = number of treatments             \n"
00097      "-dset 1 filename               data set for treatment #1            \n"
00098      " . . .                           . . .                              \n"
00099      "-dset 1 filename               data set for treatment #1            \n"
00100      " . . .                           . . .                              \n"
00101      "-dset s filename               data set for treatment #s            \n"
00102      " . . .                           . . .                              \n"
00103      "-dset s filename               data set for treatment #s            \n"
00104      "                                                                    \n"
00105      "[-workmem mega]                number of megabytes of RAM to use    \n"
00106      "                                 for statistical workspace          \n"
00107      "[-voxel num]                   screen output for voxel # num        \n"
00108      "-out prefixnamename            Kruskal-Wallis statistics are written\n"
00109      "                                 to file prefixname                 \n"
00110      "\n");
00111   
00112   
00113   printf
00114     (
00115      "\n"
00116      "N.B.: For this program, the user must specify 1 and only 1 sub-brick  \n"
00117      "      with each -dset command. That is, if an input dataset contains  \n"
00118      "      more than 1 sub-brick, a sub-brick selector must be used, e.g.: \n"
00119      "      -dset 2 'fred+orig[3]'                                          \n"
00120      );
00121 
00122   printf("\n" MASTER_SHORTHELP_STRING ) ;
00123   
00124   exit(0);
00125 }

00174 {
00175   int nopt = 1;                  /* input option argument counter */
00176   int ival;                      /* integer input */
00177   int nijk;                      /* count of data files */     
00178   float fval;                    /* float input */
00179   THD_3dim_dataset * dset=NULL;             /* test whether data set exists */
00180   char message[MAX_NAME_LENGTH];            /* error message */
00181 
00182 
00183   /*----- does user request help menu? -----*/
00184   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();
00185 
00186   
00187   /*----- add to program log -----*/
00188   AFNI_logger (PROGRAM_NAME,argc,argv); 
00189 
00190 
00191   /*----- initialize the input options -----*/
00192   initialize_options (option_data);
00193 
00194   
00195   /*----- main loop over input options -----*/
00196   while (nopt < argc)
00197     {
00198       
00199       
00200       /*-----   -datum type   -----*/
00201       if( strncmp(argv[nopt],"-datum",6) == 0 ){
00202         if( ++nopt >= argc ) NP_error("need an argument after -datum!") ;
00203         
00204         if( strcmp(argv[nopt],"short") == 0 ){
00205           option_data->datum = MRI_short ;
00206         } else if( strcmp(argv[nopt],"float") == 0 ){
00207           option_data->datum = MRI_float ;
00208         } else {
00209           char buf[256] ;
00210           sprintf(buf,
00211                   "-datum of type '%s' is not supported in 3dKruskalWallis.",
00212                   argv[nopt] ) ;
00213           NP_error(buf) ;
00214         }
00215         nopt++ ; continue ;  /* go to next arg */
00216       }
00217       
00218       
00219       /*-----   -session dirname    -----*/
00220       if( strncmp(argv[nopt],"-session",6) == 0 ){
00221         nopt++ ;
00222         if( nopt >= argc ) NP_error("need argument after -session!") ;
00223         strcpy(option_data->session , argv[nopt++]) ;
00224         continue ;
00225       }
00226       
00227       
00228       /*-----   -voxel num  -----*/
00229       if (strncmp(argv[nopt], "-voxel", 6) == 0)
00230         {
00231           nopt++;
00232           if (nopt >= argc)  NP_error ("need argument after -voxel ");
00233           sscanf (argv[nopt], "%d", &ival);
00234           if (ival <= 0)
00235             NP_error ("illegal argument after -voxel ");
00236           option_data->nvoxel = ival;
00237           nopt++;
00238           continue;
00239         }
00240       
00241       
00242       /*-----   -workmem megabytes  -----*/
00243 
00244       if( strncmp(argv[nopt],"-workmem",6) == 0 ){
00245          nopt++ ;
00246          if( nopt >= argc ) NP_error ("need argument after -workmem!") ;
00247          sscanf (argv[nopt], "%d", &ival);
00248          if( ival <= 0 ) NP_error ("illegal argument after -workmem!") ;
00249          option_data->workmem = ival ;
00250          nopt++ ; continue ;
00251       }
00252 
00253 
00254       /*-----   -levels s  -----*/
00255       if (strncmp(argv[nopt], "-levels", 7) == 0)
00256         {
00257           nopt++;
00258           if (nopt >= argc)  NP_error ("need argument after -levels ");
00259           sscanf (argv[nopt], "%d", &ival);
00260           if ((ival <= 0) || (ival > MAX_TREATMENTS))
00261             NP_error ("illegal argument after -levels ");
00262           option_data->s = ival;
00263           nopt++;
00264           continue;
00265         }
00266       
00267       
00268       /*-----   -dset level filename   -----*/
00269       if (strncmp(argv[nopt], "-dset", 5) == 0)
00270         {
00271           nopt++;
00272           if (nopt+1 >= argc)  NP_error ("need 2 arguments after -dset ");
00273           sscanf (argv[nopt], "%d", &ival);
00274           if ((ival <= 0) || (ival > option_data->s))
00275             NP_error ("illegal argument after -dset ");
00276           
00277           option_data->n[ival-1] += 1;
00278 
00279           if (option_data->n[ival-1] > MAX_OBSERVATIONS)
00280             NP_error ("too many data files");
00281           nijk = option_data->n[ival-1];
00282           
00283           /*--- check whether input files exist ---*/
00284           nopt++;
00285           dset = THD_open_dataset( argv[nopt] ) ;
00286           if( ! ISVALID_3DIM_DATASET(dset) )
00287             {
00288              sprintf(message,"Unable to open dataset file %s\n", argv[nopt]);
00289              NP_error (message);
00290             }
00291 
00292           /*--- check number of selected sub-bricks ---*/
00293           if (DSET_NVALS(dset) != 1)
00294             {
00295               sprintf(message,"Must specify exactly 1 sub-brick for file %s\n",
00296                       argv[nopt]);
00297               NP_error (message);
00298             }
00299 
00300           THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00301           
00302           option_data->xname[ival-1][nijk-1] 
00303             =  malloc (sizeof(char) * MAX_NAME_LENGTH);
00304           strcpy (option_data->xname[ival-1][nijk-1], argv[nopt]);
00305           nopt++;
00306           continue;
00307         }
00308       
00309       
00310       /*-----   -out filename   -----*/
00311       if (strncmp(argv[nopt], "-out", 4) == 0)
00312         {
00313           nopt++;
00314           if (nopt >= argc)  NP_error ("need argument after -out ");
00315           option_data->outfile = malloc (sizeof(char) * MAX_NAME_LENGTH);
00316           strcpy (option_data->outfile, argv[nopt]);
00317           nopt++;
00318           continue;
00319         }
00320             
00321       
00322       /*----- unknown command -----*/
00323       NP_error ("unrecognized command line option ");
00324     }
00325 
00326 }

00368 {
00369   
00370   
00371   /*----- allocate memory space for input data -----*/   
00372   *option_data = (NP_options *) malloc(sizeof(NP_options));
00373   if (*option_data == NULL)
00374     NP_error ("memory allocation error");
00375   
00376   /*----- get command line inputs -----*/
00377   get_options(argc, argv, *option_data);
00378   
00379   /*----- use first data set to get data set dimensions -----*/
00380   (*option_data)->first_dataset = (*option_data)->xname[0][0];
00381   get_dimensions (*option_data);
00382   printf ("Data set dimensions:  nx = %d  ny = %d  nz = %d  nxyz = %d \n",
00383           (*option_data)->nx, (*option_data)->ny,
00384           (*option_data)->nz, (*option_data)->nxyz);
00385   
00386 
00387   /*----- check for valid inputs -----*/
00388   check_for_valid_inputs (*option_data);
00389     
00390   /*----- check whether output files already exist -----*/
00391   check_one_output_file (*option_data, (*option_data)->outfile);
00392 
00393   /*----- allocate memory -----*/
00394   *best = (float *) malloc(sizeof(float) * (*option_data)->nxyz);
00395   if (*best == NULL)
00396     NP_error ("memory allocation error");
00397   *kstat = (float *) malloc(sizeof(float) * (*option_data)->nxyz);
00398   if (*kstat == NULL)
00399     NP_error ("memory allocation error");
00400  
00401   
00402 }

00134 {
00135   int i;          /* index */
00136   
00137   option_data->datum = ILLEGAL_TYPE;
00138   strcpy (option_data->session, "./");
00139  
00140 
00141   option_data->nvoxel = -1;
00142   
00143   option_data->s = 0;
00144   
00145   for (i = 0;  i < MAX_TREATMENTS;  i++)
00146     option_data->n[i] = 0;
00147 
00148   option_data->workmem = 12;
00149  
00150   /*----- allocate memory for storing data file names -----*/
00151   option_data->xname = (char ***) malloc (sizeof(char **) * MAX_TREATMENTS);
00152   for (i = 0;  i < MAX_TREATMENTS;  i++)
00153     option_data->xname[i]
00154       = (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
00155 
00156   option_data->first_dataset = NULL;
00157   
00158   option_data->nx = 0;
00159   option_data->ny = 0;
00160   option_data->nz = 0;
00161   option_data->nxyz = 0;
00162 
00163   option_data->outfile = NULL;
00164 
00165 }

00789 {
00790   NP_options * option_data = NULL;    /* user input options */
00791   float * best;                       /* index of best treatment */
00792   float * kstat;                      /* Kruskal-Wallis statistic */
00793 
00794   
00795   /*----- Identify software -----*/
00796   printf ("\n\n");
00797   printf ("Program: %s \n", PROGRAM_NAME);
00798   printf ("Author:  %s \n", PROGRAM_AUTHOR); 
00799   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
00800   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
00801   printf ("\n");
00802 
00803    /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
00804 
00805    machdep() ; 
00806    { int new_argc ; char ** new_argv ;
00807      addto_args( argc , argv , &new_argc , &new_argv ) ;
00808      if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00809    }
00810 
00811 
00812   /*----- program initialization -----*/
00813   initialize (argc, argv, &option_data, &best, &kstat);
00814   
00815   /*----- calculate nonparameteric Kruskal-Wallis statistics -----*/
00816   calculate_results (option_data, best, kstat);
00817   
00818   /*----- generate requested output -----*/
00819   output_results (argc, argv, option_data, best, kstat);
00820 
00821   /*----- terminate program -----*/
00822   terminate (&option_data, &best, &kstat);
00823 
00824   exit(0);
00825 }

00727 {
00728 
00729   /*----- write out afni fict data file -----*/
00730   write_afni_fict (argc, argv, option_data, option_data->outfile, 
00731                    best, kstat, option_data->s - 1);
00732 
00733 }

00544 {
00545   int i;                             /* treatment index */
00546   int j;                             /* array index */
00547 
00548 
00549   /*----- check for voxel output  -----*/
00550   if (nvox > 0)
00551     {
00552       printf ("\nResults for voxel #%d : \n\n", nvox);
00553 
00554       for (i = 0;  i < s;  i++)
00555         {
00556           printf ("Y%d data:  ", i+1);
00557           for (j = 0;  j < n[i];  j++)
00558             {
00559               printf (" %6.1f", xarray[i][j]);
00560               if (((j+1) % 8 == 0) && (j < n[i]-1))
00561                 printf ("\n          ");
00562             }
00563           printf ("\n");
00564           if (n[i] > 8)  printf ("\n");
00565         }
00566       if (n[s] <= 8)  printf ("\n");
00567     }
00568 
00569 
00570   /*----- calculate Kruskal-Wallis statistic -----*/
00571   calc_stat (nvox, s, n, xarray, best, kstat);
00572 
00573 }

00749 {
00750    int i, j;                       /* dataset indices */
00751 
00752 
00753    /*----- deallocate memory -----*/
00754    for (i = 0;  i < (*option_data)->s;  i++)
00755      for (j = 0;  j < (*option_data)->n[i];  j++)
00756        {
00757          free ((*option_data)->xname[i][j]);
00758          (*option_data)->xname[i][j] = NULL;
00759        }
00760    for (i = 0;  i < (*option_data)->s;  i++)
00761      {
00762        free ((*option_data)->xname[i]);
00763        (*option_data)->xname[i] = NULL;
00764      }
00765    free ((*option_data)->xname);
00766    (*option_data)->xname = NULL;
00767 
00768    if ((*option_data)->outfile != NULL)
00769    {
00770       free ((*option_data)-> outfile);
00771       (*option_data)->outfile = NULL;
00772    }
00773 
00774    free (*option_data);   *option_data = NULL;
00775 
00776    free (*best);          *best = NULL;
00777 
00778    free (*kstat);         *kstat = NULL;
00779 }