00418 {
00419   const float EPSILON = 1.0e-10;      /* protection from roundoff error */
00420   int i, j;                   /* array indices */
00421   node * head = NULL;         /* points to head of list */        
00422   node * ptr = NULL;          /* points to current position in list */
00423   int NN;                     /* total number of sample points */
00424   float rsum;                 /* sum of squares of ranks */
00425   int d;                      /* count of number of ties */ 
00426   float corr;                 /* correction to account for ties */
00427   float rank;                 /* rank of data point */
00428   float ranksum;              /* sum of ranks for ith treatment */
00429   float qnum;                 /* numerator of Friedman statistic */
00430   float qden;                 /* denominator of Friedman statistic */
00431   float best_rank;            /* best average rank for a treatment */
00432   float ** rank_array;        /* array to store ranks for all observations */
00433 
00434 
00435 
00436   /*----- allocate memory for storing ranks -----*/
00437   rank_array = (float **) malloc (sizeof(float *) * s);
00438   for (i = 0;  i < s;  i++)
00439     rank_array[i] = (float *) malloc (sizeof(float) * n);
00440 
00441 
00442   /*----- loop over blocks -----*/
00443   corr = 0.0;
00444   for (j = 0;  j < n;  j++)
00445     {
00446 
00447       /*----- enter and sort data for each treatment within block j -----*/
00448       for (i = 0;  i < s;  i++)
00449         node_addvalue (&head, xarray[i][j]);
00450 
00451 
00452       /*----- store the ranks for each treatment within block j -----*/
00453       for (i = 0;  i < s;  i++)
00454         rank_array[i][j] = node_get_rank (head, xarray[i][j]);
00455 
00456       
00457       /*----- calculate the ties correction factor -----*/
00458       ptr = head;
00459       while (ptr != NULL)
00460         {
00461           d = ptr->d;
00462           corr += d*d*d - d;
00463           ptr = ptr->next;
00464         }
00465 
00466       list_delete (&head);
00467 
00468     } /* j loop */
00469 
00470 
00471   /*----- if display voxel, write the ranks of the input data -----*/
00472   if (nvox > 0)
00473     {
00474       printf ("\n");
00475       for (i = 0;  i < s;  i++)
00476         {
00477           printf ("Y%d ranks: ", i+1);
00478           for (j = 0;  j < n;  j++)
00479             {
00480               rank = rank_array[i][j];
00481               printf (" %6.1f", rank);
00482               if (((j+1) % 8 == 0) && (j < n-1))
00483                 printf ("\n          ");
00484             }
00485           printf ("\n");
00486           if (n > 8)  printf ("\n");
00487         }
00488       printf ("\n");
00489       for (i = 0;  i < s;  i++)
00490         {
00491           printf ("Y%d: ", i+1);
00492           ranksum = 0.0;
00493           for (j = 0;  j < n;  j++)
00494             {
00495               rank = rank_array[i][j];
00496               ranksum += rank;
00497             }
00498           printf ("   Rank sum = %6.1f    Rank average = %6.1f \n", 
00499                   ranksum, ranksum/n);
00500         }
00501       printf ("\n");
00502     }
00503 
00504 
00505   /*----- calculate the sum of the ranks -----*/
00506   rsum = 0.0;
00507   *best = 0.0;
00508   best_rank = (s + 1.0) / 2.0 + EPSILON;
00509   for (i = 0;  i < s;  i++)
00510     {
00511       ranksum = 0.0;
00512       for (j = 0;  j < n;  j++)
00513         ranksum += rank_array[i][j];
00514       rsum += ranksum * ranksum;
00515 
00516       if (ranksum/n > best_rank)
00517         {
00518           *best = (float) (i+1);
00519           best_rank = ranksum / n;
00520         }
00521     }
00522 
00523 
00524   /*----- numerator of Friedman statistic -----*/
00525   qnum = (12.0/(n*s*(s+1)))*rsum - 3.0*n*(s+1);
00526 
00527   /*----- denominator of Friedman statistic -----*/
00528   qden = 1.0 - (corr / (n*s*(s*s-1)));
00529 
00530   /*----- calculate Friedman statistic -----*/
00531   if (qden < EPSILON)
00532     *qstat = 0.0;
00533   else
00534     *qstat = qnum / qden;
00535   if (nvox > 0)  printf ("Q = %f \n", *qstat);
00536 
00537 
00538   /*----- deallocate memory -----*/
00539   for (i = 0;  i < s;  i++)
00540     {
00541       free (rank_array[i]);
00542       rank_array[i] = NULL;
00543     }
00544   free (rank_array);
00545   rank_array = NULL;
00546   
00547 }

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

00332 {
00333   int i, n;
00334   char message[MAX_NAME_LENGTH];            /* error message */
00335 
00336 
00337   n = option_data->n[0];
00338   if (n < 1)
00339     NP_error ("Sample size is too small");
00340   
00341   for (i = 1;  i < option_data->s;  i++)
00342     if (option_data->n[i] != n)
00343       NP_error ("Must have equal sample sizes for all treatments");
00344      
00345 
00346   if (option_data->nvoxel > option_data->nxyz)
00347     NP_error ("argument of -voxel is too large");
00348 
00349 }

00087 {
00088   printf 
00089     (
00090      "This program performs nonparametric Friedman test for               \n"
00091      "randomized complete block design experiments.                     \n\n"
00092      "Usage:                                                              \n"
00093      "3dFriedman                                                          \n"
00094      "-levels s                      s = number of treatments             \n"
00095      "-dset 1 filename               data set for treatment #1            \n"
00096      " . . .                           . . .                              \n"
00097      "-dset 1 filename               data set for treatment #1            \n"
00098      " . . .                           . . .                              \n"
00099      "-dset s filename               data set for treatment #s            \n"
00100      " . . .                           . . .                              \n"
00101      "-dset s filename               data set for treatment #s            \n"
00102      "                                                                    \n"
00103      "[-workmem mega]                number of megabytes of RAM to use    \n"
00104      "                                 for statistical workspace          \n"
00105      "[-voxel num]                   screen output for voxel # num        \n"
00106      "-out prefixname                Friedman statistics are written      \n"
00107      "                                 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 }

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

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

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->s = 0;
00141   
00142   for (i = 0;  i < MAX_TREATMENTS;  i++)
00143     option_data->n[i] = 0;
00144 
00145   option_data->workmem = 12;
00146  
00147   /*----- allocate memory for storing data file names -----*/
00148   option_data->xname = (char ***) malloc (sizeof(char **) * MAX_TREATMENTS);
00149   for (i = 0;  i < MAX_TREATMENTS;  i++)
00150     option_data->xname[i]
00151       = (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
00152 
00153   option_data->first_dataset = NULL;
00154   
00155   option_data->nx = 0;
00156   option_data->ny = 0;
00157   option_data->nz = 0;
00158   option_data->nxyz = 0;
00159 
00160   option_data->outfile = NULL;
00161 
00162 }

00804 {
00805   NP_options * option_data = NULL;    /* user input options */
00806   float * best;                       /* index of best treatment */
00807   float * qstat;                      /* Friedman statistic */
00808  
00809   
00810   /*----- Identify software -----*/
00811   printf ("\n\n");
00812   printf ("Program: %s \n", PROGRAM_NAME);
00813   printf ("Author:  %s \n", PROGRAM_AUTHOR); 
00814   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
00815   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
00816   printf ("\n");
00817 
00818 
00819    /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
00820 
00821    machdep() ; 
00822    { int new_argc ; char ** new_argv ;
00823      addto_args( argc , argv , &new_argc , &new_argv ) ;
00824      if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00825    }
00826 
00827 
00828   /*----- program initialization -----*/
00829   initialize (argc, argv, &option_data, &best, &qstat);
00830   
00831   /*----- calculate nonparameteric Friedman statistics -----*/
00832   calculate_results (option_data, best, qstat);
00833   
00834   /*----- generate requested output -----*/
00835   output_results (argc, argv, option_data, best, qstat);
00836 
00837   /*----- terminate program -----*/
00838   terminate (&option_data, &best, &qstat);
00839 
00840 }

00742 {
00743 
00744   /*----- write out afni fict data file -----*/
00745   write_afni_fict (argc, argv, option_data, option_data->outfile, 
00746                    best, qstat, option_data->s - 1);
00747 
00748 }

00565 {
00566   int i;                             /* treatment index */
00567   int j;                             /* array index */
00568 
00569 
00570   /*----- check for voxel output  -----*/
00571   if (nvox > 0)
00572     {
00573       printf ("\nResults for voxel #%d : \n\n", nvox);
00574 
00575       for (i = 0;  i < s;  i++)
00576         {
00577           printf ("Y%d data:  ", i+1);
00578           for (j = 0;  j < n;  j++)
00579             {
00580             printf (" %6.1f", xarray[i][j]);
00581             if (((j+1) % 8 == 0) && (j < n-1))
00582               printf ("\n          ");
00583             }
00584           printf ("\n");
00585           if (n > 8)  printf ("\n");
00586         }
00587       if (n <= 8)  printf ("\n");
00588     }
00589 
00590 
00591   /*----- calculate Friedman statistic -----*/
00592   calc_stat (nvox, s, n, xarray, best, qstat);
00593 
00594 }

00764 {
00765    int i, j;                       /* dataset indices */
00766 
00767 
00768    /*----- deallocate memory -----*/
00769    for (i = 0;  i < (*option_data)->s;  i++)
00770      for (j = 0;  j < (*option_data)->n[i];  j++)
00771        {
00772          free ((*option_data)->xname[i][j]);
00773          (*option_data)->xname[i][j] = NULL;
00774        }
00775    for (i = 0;  i < (*option_data)->s;  i++)
00776      {
00777        free ((*option_data)->xname[i]);
00778        (*option_data)->xname[i] = NULL;
00779      }
00780    free ((*option_data)->xname);
00781    (*option_data)->xname = NULL;
00782 
00783    if ((*option_data)->outfile != NULL)
00784    {
00785       free ((*option_data)-> outfile);
00786       (*option_data)->outfile = NULL;
00787    }
00788 
00789    free (*option_data);   *option_data = NULL;
00790 
00791    free (*best);          *best = NULL;
00792 
00793    free (*qstat);         *qstat = NULL;
00794 }