00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1999-2003, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   This sample program generates random stimulus functions.
00010 
00011   File:    RSFgen.c
00012   Author:  B. Douglas Ward
00013   Date:    06 July 1999
00014 
00015 
00016   Mod:     Increased max. allowed number of input stimulus functions.
00017   Date:    24 August 1999
00018 
00019   Mod:     Added option to create multiple column stimulus function files.
00020   Date:    24 November 1999
00021 
00022   Mod:     Changed option label "-nstim" to "-num_stimts".
00023   Date:    29 November 1999
00024 
00025   Mod:     Added option "-nblock" to specify block length for each stim fn.
00026   Date:    15 December 1999
00027 
00028   Mod:     Added flag to expand array for block type design.
00029   Date:    14 January 2000
00030 
00031   Mod:     Added -markov and -pzero options.
00032   Date:    03 October 2000
00033 
00034   Mod:     Added -pseed option for permutation of stimulus labels.
00035   Date:    27 April 2001
00036 
00037   Mod:     Changes to eliminate constraints on number of stimulus time series.
00038   Date:    11 May 2001
00039 
00040   Mod:     Added call to AFNI_logger.
00041   Date:    15 August 2001
00042 
00043   Mod:     Made -nblock option compatible with the Markov Chain options.
00044   Date:    06 March 2002
00045 
00046   Mod:     Added -one_col option to write stimulus functions as a single
00047            column of decimal integers.
00048   Date:    18 April 2002
00049 
00050   Mod:     Added -quiet option.
00051   Date:    30 December 2002
00052 
00053   Mod:     Added -table option, to generate random permutations of the rows
00054            of an input column or table of numbers.
00055   Date:    13 March 2003
00056 
00057 */
00058 
00059 /*---------------------------------------------------------------------------*/
00060 
00061 #define PROGRAM_NAME "RSFgen"                        /* name of this program */
00062 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
00063 #define PROGRAM_INITIAL "06 July 1999"    /* date of initial program release */
00064 #define PROGRAM_LATEST "13 March 2003"    /* date of latest program revision */
00065 
00066 /*---------------------------------------------------------------------------*/
00067 
00068 #include <math.h>
00069 #include <stdlib.h>
00070 #include <stdio.h>
00071 
00072 #include "mrilib.h"
00073 #include "matrix.h"
00074 
00075 #include "randgen.c"
00076 #include "matrix.c"
00077 
00078 
00079 /*---------------------------------------------------------------------------*/
00080 /*
00081   Global variables.
00082 */
00083 
00084 
00085 int NT = 0;             /* length of stimulus time series */
00086 int nt = 0;             /* length of simple time series (block length = 1) */
00087 int num_stimts = 0;     /* number of input stimuli (experimental conditions) */
00088 int * num_reps = NULL;  /* number of repetitions for each stimulus */
00089 int * nblock = NULL;    /* block length for each stimulus */
00090 int expand = 0;         /* flag to expand the array for block type design */
00091 long seed = 1234567;    /* random number seed */
00092 long pseed = 0;         /* stimulus permutation random number seed */
00093 char * prefix = NULL;   /* prefix for output .1D stimulus functions */
00094 int one_file = 0;       /* flag for place stim functions into a single file */
00095 int one_col = 0;        /* flag for write stim functions as a single column */
00096 int markov = 0;         /* flag for Markov process */
00097 char * tpm_file = NULL; /* file name for input of transition prob. matrix */
00098 float pzero = 0.0;      /* zero (null) state probability */
00099 int quiet = 0;          /* flag for suppress screen output */
00100 char * table_file = NULL;  /* input time series file */
00101 
00102 
00103 /*---------------------------------------------------------------------------*/
00104 /*
00105    Print error message and stop.
00106 */
00107 
00108 void RSF_error (char * message)
00109 {
00110   fprintf (stderr, "\n%s Error: %s \n", PROGRAM_NAME, message);
00111   exit(1);
00112 }
00113 
00114 
00115 /*---------------------------------------------------------------------------*/
00116 
00117 /** macro to test a malloc-ed pointer for validity **/
00118 
00119 #define MTEST(ptr) \
00120 if((ptr)==NULL) \
00121 ( RSF_error ("Cannot allocate memory") )
00122      
00123 
00124 /*---------------------------------------------------------------------------*/
00125 /*
00126   Identify software.
00127 */
00128 
00129 void identify_software ()
00130 {
00131   
00132   /*----- Identify software -----*/
00133   printf ("\n\n");
00134   printf ("Program:          %s \n", PROGRAM_NAME);
00135   printf ("Author:           %s \n", PROGRAM_AUTHOR); 
00136   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
00137   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
00138   printf ("\n");
00139 }
00140 
00141 
00142 /*---------------------------------------------------------------------------*/
00143 /*
00144   Routine to display help menu for program RSFgen.
00145 */
00146 
00147 void display_help_menu()
00148 {
00149   identify_software();
00150 
00151   printf (
00152     "Sample program to generate random stimulus functions.                  \n"
00153     "                                                                       \n"
00154     "Usage:                                                                 \n"
00155     PROGRAM_NAME "                                                          \n"
00156     "-nt n            n = length of time series                             \n"
00157     "-num_stimts p    p = number of input stimuli (experimental conditions) \n"
00158     "[-nblock i k]    k = block length for stimulus i  (1<=i<=p)            \n"
00159     "                     (default: k = 1)                                  \n"
00160     "[-seed s]        s = random number seed                                \n"
00161     "[-quiet]         flag to suppress screen output                        \n"
00162     "[-one_file]      place stimulus functions into a single .1D file       \n"
00163     "[-one_col]       write stimulus functions as a single column of decimal\n"
00164     "                   integers (default: multiple columns of binary nos.) \n"
00165     "[-prefix pname]  pname = prefix for p output .1D stimulus functions    \n"
00166     "                   e.g., pname1.1D, pname2.1D, ..., pnamep.1D          \n"
00167     "                                                                       \n"
00168     "The following Random Permutation, Markov Chain, and Input Table options\n"
00169     "are mutually exclusive.                                                \n"
00170     "                                                                       \n"
00171     "Random Permutation options:                                            \n"
00172     "-nreps i r       r = number of repetitions for stimulus i  (1<=i<=p)   \n"
00173     "[-pseed s]       s = stim label permutation random number seed         \n"
00174     "                                     p                                 \n"
00175     "                 Note: Require n >= Sum (r[i] * k[i])                  \n"
00176     "                                    i=1                                \n"
00177     "                                                                       \n"
00178     "Markov Chain options:                                                  \n"
00179     "-markov mfile    mfile = file containing the transition prob. matrix   \n"
00180     "[-pzero z]       probability of a zero (i.e., null) state              \n"
00181     "                     (default: z = 0)                                  \n"
00182     "                                                                       \n"
00183     "Input Table row permutation options:                                   \n"
00184     "[-table dfile]   dfile = filename of column or table of numbers        \n"
00185     "                 Note: dfile may have a column selector attached       \n"
00186     "                 Note: With this option, all other input options,      \n"
00187     "                       except -seed and -prefix, are ignored           \n"
00188     "                                                                       \n"
00189     "                                                                       \n"
00190     "Warning: This program will overwrite pre-existing .1D files            \n"
00191     "                                                                       \n"
00192     );
00193   
00194   exit(0);
00195 }
00196 
00197 
00198 /*---------------------------------------------------------------------------*/
00199 /*
00200   Routine to get user specified input options.
00201 */
00202 
00203 void get_options
00204 (
00205   int argc,                        /* number of input arguments */
00206   char ** argv                     /* array of input arguments */ 
00207 )
00208 
00209 {
00210   int nopt = 1;                     /* input option argument counter */
00211   int ival;                         /* integer input */
00212   float fval;                       /* float input */
00213   long lval;                        /* long integer input */
00214   char message[THD_MAX_NAME];       /* error message */
00215   int i;
00216 
00217 
00218   /*----- Does user request help menu? -----*/
00219   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();  
00220 
00221   
00222   /*----- add to program log -----*/
00223   AFNI_logger (PROGRAM_NAME,argc,argv); 
00224 
00225 
00226   /*----- Main loop over input options -----*/
00227   while (nopt < argc )
00228     {
00229      
00230       /*-----  -nt n  -----*/
00231       if (strncmp(argv[nopt], "-nt", 3) == 0)
00232         {
00233           nopt++;
00234           if (nopt >= argc)  RSF_error ("need argument after -nt ");
00235           sscanf (argv[nopt], "%d", &ival);
00236           if (ival <= 0)
00237             RSF_error ("illegal argument after -nt ");
00238           NT = ival;
00239           nopt++;
00240           continue;
00241         }
00242 
00243 
00244       /*-----  -num_stimts p  -----*/
00245       if (strncmp(argv[nopt], "-num_stimts", 11) == 0)
00246         {
00247           nopt++;
00248           if (nopt >= argc)  RSF_error ("need argument after -num_stimts ");
00249           sscanf (argv[nopt], "%d", &ival);
00250           if (ival <= 0)
00251             RSF_error ("illegal argument after -num_stimts ");
00252           num_stimts = ival;
00253 
00254           /*----- Initialize repetition number array -----*/
00255           num_reps = (int *) malloc (sizeof(int) * num_stimts);
00256           MTEST (num_reps);
00257           for (i = 0;  i < num_stimts;  i++)
00258             num_reps[i] = 0;
00259 
00260           /*----- Initialize block length array -----*/
00261           nblock = (int *) malloc (sizeof(int) * num_stimts);
00262           MTEST (nblock);
00263           for (i = 0;  i < num_stimts;  i++)
00264             nblock[i] = 1;
00265 
00266           nopt++;
00267           continue;
00268         }
00269 
00270 
00271       /*-----  -nreps i r  -----*/
00272       if (strncmp(argv[nopt], "-nreps", 6) == 0)
00273         {
00274           nopt++;
00275           if (nopt+1 >= argc)  RSF_error ("need 2 arguments after -nreps ");
00276           sscanf (argv[nopt], "%d", &ival);
00277           if ((ival <= 0) || (ival > num_stimts))
00278             RSF_error ("illegal i argument for -nreps i r ");
00279           i = ival - 1;
00280           nopt++;
00281 
00282           sscanf (argv[nopt], "%d", &ival);
00283           if (ival <= 0)
00284             RSF_error ("illegal r argument for -nreps i r ");
00285           num_reps[i] = ival;
00286           nopt++;
00287           continue;
00288         }
00289 
00290 
00291       /*-----  -nblock i k  -----*/
00292       if (strncmp(argv[nopt], "-nblock", 7) == 0)
00293         {
00294           nopt++;
00295           if (nopt+1 >= argc)  RSF_error ("need 2 arguments after -nblock ");
00296           sscanf (argv[nopt], "%d", &ival);
00297           if ((ival <= 0) || (ival > num_stimts))
00298             RSF_error ("illegal i argument for -nblock i k ");
00299           i = ival - 1;
00300           nopt++;
00301 
00302           sscanf (argv[nopt], "%d", &ival);
00303           if (ival <= 0)
00304             RSF_error ("illegal k argument for -nblock i k ");
00305           nblock[i] = ival;
00306           if (ival > 1)  expand = 1;
00307           nopt++;
00308           continue;
00309         }
00310 
00311 
00312       /*-----  -seed s  -----*/
00313       if (strncmp(argv[nopt], "-seed", 5) == 0)
00314         {
00315           nopt++;
00316           if (nopt >= argc)  RSF_error ("need argument after -seed ");
00317           sscanf (argv[nopt], "%ld", &lval);
00318           if (lval <= 0)
00319             RSF_error ("illegal argument after -seed ");
00320           seed = lval;
00321           nopt++;
00322           continue;
00323         }
00324 
00325 
00326       /*-----  -pseed s  -----*/
00327       if (strcmp(argv[nopt], "-pseed") == 0)
00328         {
00329           nopt++;
00330           if (nopt >= argc)  RSF_error ("need argument after -pseed ");
00331           sscanf (argv[nopt], "%ld", &lval);
00332           if (lval <= 0)
00333             RSF_error ("illegal argument after -pseed ");
00334           pseed = lval;
00335           nopt++;
00336           continue;
00337         }
00338 
00339 
00340       /*-----  -quiet  -----*/
00341       if (strcmp(argv[nopt], "-quiet") == 0)
00342         {
00343           quiet = 1;
00344           nopt++;
00345           continue;
00346         }
00347 
00348 
00349       /*-----  -one_file  -----*/
00350       if (strncmp(argv[nopt], "-one_file", 9) == 0)
00351         {
00352           one_file = 1;
00353           nopt++;
00354           continue;
00355         }
00356 
00357 
00358       /*-----  -one_col  -----*/
00359       if (strcmp(argv[nopt], "-one_col") == 0)
00360         {
00361           one_col = 1;
00362           nopt++;
00363           continue;
00364         }
00365 
00366 
00367       /*-----   -prefix pname   -----*/
00368       if (strncmp(argv[nopt], "-prefix", 7) == 0)
00369         {
00370           nopt++;
00371           if (nopt >= argc)  RSF_error ("need argument after -prefix ");
00372           prefix = AFMALL(char, sizeof(char) * THD_MAX_NAME);
00373           MTEST (prefix);
00374           strcpy (prefix, argv[nopt]);
00375           nopt++;
00376           continue;
00377         }
00378       
00379 
00380       /*-----   -markov mfile   -----*/
00381       if (strcmp(argv[nopt], "-markov") == 0)
00382         {
00383           markov = 1;
00384           nopt++;
00385           if (nopt >= argc)  RSF_error ("need argument after -markov ");
00386           tpm_file = AFMALL(char, sizeof(char) * THD_MAX_NAME);
00387           MTEST (tpm_file);
00388           strcpy (tpm_file, argv[nopt]);
00389           nopt++;
00390           continue;
00391         }
00392       
00393 
00394       /*-----   -pzero z   -----*/
00395       if (strcmp(argv[nopt], "-pzero") == 0)
00396         {
00397           nopt++;
00398           if (nopt >= argc)  RSF_error ("need argument after -pzero ");
00399           sscanf (argv[nopt], "%f", &fval);
00400           if ((fval < 0.0) || (fval > 1.0))  
00401             RSF_error ("Require  0.0 <= pzero <= 1.0");
00402           pzero = fval;
00403           nopt++;
00404           continue;
00405         }
00406       
00407 
00408       /*-----   -table dfile   -----*/
00409       if (strcmp(argv[nopt], "-table") == 0)
00410         {
00411           nopt++;
00412           if (nopt >= argc)  RSF_error ("need argument after -table ");
00413           table_file = AFMALL(char, sizeof(char) * THD_MAX_NAME);
00414           MTEST (table_file);
00415           strcpy (table_file, argv[nopt]);
00416           nopt++;
00417           continue;
00418         }
00419       
00420 
00421       /*----- Unknown command -----*/
00422       sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00423       RSF_error (message);
00424       
00425     }
00426 }
00427 
00428 
00429 /*---------------------------------------------------------------------------*/
00430 /*
00431   Read input table.
00432 */
00433 
00434 void read_table 
00435 (
00436   char * table_file,       /* file containing input table */
00437   MRI_IMAGE ** flim        /* data structure containing input table */
00438 )
00439 
00440 {
00441   int i;
00442   char message[THD_MAX_NAME];  /* error message */
00443 
00444   /*----- Read table -----*/
00445   *flim = mri_read_1D(table_file);
00446   if ((*flim) == NULL)  
00447     {
00448       sprintf (message,  "Unable to read table file: %s", table_file);
00449       RSF_error (message);
00450     }
00451 
00452   /*----- Initialize control variables -----*/
00453   NT = (*flim)->nx;
00454   num_stimts = NT;
00455   markov = 0;
00456   pseed = 0;
00457 
00458   /*----- Initialize repetition number array -----*/
00459   num_reps = (int *) malloc (sizeof(int) * num_stimts);
00460   MTEST (num_reps);
00461   for (i = 0;  i < num_stimts;  i++)
00462     num_reps[i] = 1;
00463   
00464   /*----- Initialize block length array -----*/
00465   nblock = (int *) malloc (sizeof(int) * num_stimts);
00466   MTEST (nblock);
00467   for (i = 0;  i < num_stimts;  i++)
00468     nblock[i] = 1;
00469   
00470 }
00471 
00472 
00473 /*---------------------------------------------------------------------------*/
00474 /*
00475   Print input options. 
00476 */
00477 
00478 void print_options ()
00479 
00480 {
00481   int i;
00482 
00483   identify_software();
00484   
00485   if (table_file != NULL)
00486     {
00487       printf ("table file    = %s \n", table_file);
00488       printf ("nt            = %d \n", NT);
00489       printf ("seed          = %ld \n", seed);
00490       printf ("output prefix = %s \n", prefix);
00491     }
00492   else
00493     {
00494       printf ("nt            = %d \n", NT);
00495       printf ("num_stimts    = %d \n", num_stimts);
00496       printf ("seed          = %ld \n", seed);
00497       if (pseed)  printf ("pseed         = %ld \n", pseed);
00498       printf ("output prefix = %s \n", prefix);
00499       if (markov)
00500         {
00501           printf ("TPM file      = %s \n", tpm_file);
00502           printf ("pzero         = %f \n", pzero);
00503           for (i = 0;  i < num_stimts;  i++)
00504             printf ("nblock[%d] = %d \n", i+1, nblock[i]);
00505         }
00506       else
00507         for (i = 0;  i < num_stimts;  i++)
00508           printf ("nreps[%d]  = %d    nblock[%d] = %d \n", 
00509                   i+1, num_reps[i], i+1, nblock[i]);
00510     }
00511 }
00512 
00513 
00514 /*---------------------------------------------------------------------------*/
00515 /*
00516   Perform all program initialization. 
00517 */
00518 
00519 void initialize 
00520 (  
00521   int argc,                /* number of input arguments */
00522   char ** argv,            /* array of input arguments */ 
00523   int ** darray,           /* original design array (block length = 1) */
00524   int ** earray,           /* expanded design array (arbitrary block length) */
00525   MRI_IMAGE ** flim        /* data structure containing input table */
00526 )
00527 
00528 { 
00529   int i, total;
00530 
00531 
00532   /*----- Get command line inputs -----*/
00533   get_options (argc, argv);
00534 
00535 
00536   /*----- Read input table -----*/
00537   if (table_file != NULL)  read_table (table_file, flim);
00538 
00539 
00540   /*----- Print input options -----*/
00541   if (! quiet)  print_options ();
00542 
00543 
00544   /*----- Check for valid inputs -----*/
00545   if (NT == 0)          RSF_error ("Must specify nt");
00546   if (num_stimts == 0)  RSF_error ("Must specify num_stimts");
00547   total = 0;
00548   nt = NT;
00549 
00550   if (! markov)
00551     {
00552       for (i = 0;  i < num_stimts;  i++)
00553         {
00554           if (num_reps[i] == 0)  
00555             RSF_error ("Must specify nreps > 0 for each stimulus");
00556           total += num_reps[i] * nblock[i];
00557           nt -= num_reps[i] * (nblock[i] - 1);
00558         }
00559       if (total > NT)  RSF_error ("Require  nt >= Sum (r[i] * k[i]) ");
00560     }
00561 
00562 
00563   /*----- Allocate memory for experimental design -----*/
00564   *darray = (int *) malloc (sizeof(int) * nt);   MTEST (*darray);
00565   *earray = (int *) malloc (sizeof(int) * NT);   MTEST (*earray);
00566 
00567 
00568 }
00569 
00570 
00571 /*---------------------------------------------------------------------------*/
00572 /*
00573   Use Markov chain to create random stimulus functions.
00574 */
00575 
00576 void markov_array (int * design)
00577 
00578 {
00579   int it, is, id, isprev;
00580   float prob, cumprob;
00581   matrix tpm;
00582   char message[THD_MAX_NAME];  /* error message */
00583 
00584 
00585   matrix_initialize (&tpm);
00586 
00587 
00588   /*----- Read the transition probability matrix -----*/
00589   matrix_file_read (tpm_file, num_stimts, num_stimts, &tpm, 1);
00590   if (tpm.elts == NULL)
00591     { 
00592       sprintf (message,  "Unable to read Markov chain matrix from file: %s", 
00593                tpm_file);
00594       RSF_error (message);
00595     }  
00596   if (!quiet)  matrix_sprint ("\nTPM matrix:", tpm);
00597 
00598 
00599   /*----- Verify that the TPM has the correct form -----*/
00600   for (is = 0;  is < num_stimts;  is++)
00601     {
00602       cumprob = 0.0;
00603       for (it = 0;  it < num_stimts;  it++)
00604         cumprob += tpm.elts[is][it];
00605       if (cumprob < 0.9999)  
00606         {
00607           sprintf (message, "Row %d of TPM sums to %f, which is < 1.0",
00608                    is, cumprob);
00609           RSF_error (message);
00610         }
00611       if (cumprob > 1.0001)  
00612         {
00613           sprintf (message, "Row %d of TPM sums to %f, which is > 1.0",
00614                    is, cumprob);
00615           RSF_error (message);
00616         }
00617     }
00618   
00619 
00620   /*----- Initialize the experimental design array -----*/
00621   for (it = 0;  it < NT;  it++)
00622     design[it] = 0;
00623 
00624 
00625   /*----- Initialize random number seed -----*/
00626   srand48 (seed);
00627 
00628 
00629   /*----- Generate Markov process -----*/
00630   isprev = (int) (rand_uniform(0.0,1.0)*num_stimts);
00631   it = 0;  id = 0;
00632   while (it < NT)
00633     {
00634       if ((pzero > 0.0) && (rand_uniform(0.0,1.0) < pzero))
00635         {
00636           design[id] = 0;
00637           id++;  it++;
00638         } 
00639       else
00640         {
00641           prob = rand_uniform(0.0,1.0);
00642           cumprob = 0.0;
00643           for (is = 0;  is < num_stimts;  is++)
00644             {
00645               cumprob += tpm.elts[isprev][is];
00646               if (prob <= cumprob)
00647                 {
00648                   design[id] = is+1;
00649                   isprev = is;
00650                   id++;  it += nblock[is];
00651                   break;
00652                 }
00653             }
00654         }
00655     }
00656 
00657   nt = id;
00658 
00659 
00660   matrix_destroy (&tpm);
00661 
00662   return;
00663 }
00664 
00665 
00666 /*---------------------------------------------------------------------------*/
00667 /*
00668   Fill the experimental design array.
00669 */
00670 
00671 void fill_array (int * design)
00672 
00673 {
00674   int i, is, m;
00675 
00676 
00677   for (i = 0;  i < nt;  i++)
00678     design[i] = 0;
00679 
00680   i = 0;
00681   for (is = 0;  is < num_stimts;  is++)
00682     {
00683       for (m = 0;  m < num_reps[is];  m++)
00684         {
00685           design[i] = is+1;
00686           i++;
00687         }
00688     }
00689 
00690 
00691   return;
00692 }
00693 
00694 
00695 /*---------------------------------------------------------------------------*/
00696 /*
00697   Permute the stimulus functions within the experimental design array.
00698 */
00699 
00700 void permute_array (int * design)
00701 
00702 {
00703   int i, j, temp;
00704   int is, nb;
00705 
00706 
00707   /*----- Initialize random number seed -----*/
00708   srand48 (pseed);
00709 
00710   
00711   /*----- Determine total number of blocks -----*/
00712   nb = 0;
00713   for (is = 0;  is < num_stimts;  is++)
00714     nb += num_reps[is];
00715 
00716 
00717   /*----- Permute the blocks -----*/
00718   for (i = 0;  i < nb;  i++)
00719     {
00720       j = rand_uniform(0.0,1.0) * nb;
00721 
00722       /*----- Just in case -----*/
00723       if (j < 0)  j = 0;
00724       if (j > nb-1)  j = nb-1;
00725 
00726       temp = design[i];
00727       design[i] = design[j];
00728       design[j] = temp;
00729     }
00730       
00731   return;
00732 }
00733 
00734 
00735 /*---------------------------------------------------------------------------*/
00736 /*
00737   Shuffle the experimental design array.
00738 */
00739 
00740 void shuffle_array (int * design)
00741 
00742 {
00743   int i, j, temp;
00744 
00745 
00746   /*----- Initialize random number seed -----*/
00747   srand48 (seed);
00748 
00749 
00750   for (i = 0;  i < nt;  i++)
00751     {
00752       j = rand_uniform(0.0,1.0) * nt;
00753 
00754       /*----- Just in case -----*/
00755       if (j < 0)  j = 0;
00756       if (j > nt-1)  j = nt-1;
00757 
00758       temp = design[i];
00759       design[i] = design[j];
00760       design[j] = temp;
00761     }
00762       
00763   return;
00764 }
00765 
00766 
00767 /*---------------------------------------------------------------------------*/
00768 /*
00769   Expand the experimental design array, to allow for block-type designs.
00770 */
00771 
00772 void expand_array (int * darray, int * earray)
00773 
00774 {
00775   int i, j, k, m;
00776 
00777   j = 0;
00778   for (i = 0;  i < nt, j < NT;  i++)
00779     {
00780       m = darray[i];
00781 
00782       if (m == 0)
00783         {
00784           earray[j] = 0;
00785           j++;
00786         }
00787       else
00788         {
00789           for (k = 0;  k < nblock[m-1];  k++)
00790             {
00791               earray[j] = m;
00792               j++;
00793               if (j >= NT)  break;
00794             }  
00795         }
00796     }
00797       
00798   return;
00799 }
00800 
00801 
00802 /*---------------------------------------------------------------------------*/
00803 /*
00804   Print array.
00805 */
00806 
00807 void print_array (int * array, int n)
00808 
00809 {
00810   int i;
00811 
00812   for (i = 0;  i < n;  i++)
00813     {
00814       printf (" %2d ", array[i]);
00815       if ((i+1) % 20 == 0)  printf ("\n");
00816     }
00817 
00818   printf ("\n");
00819 
00820   return;
00821 }
00822 
00823 
00824 /*---------------------------------------------------------------------------*/
00825 /*
00826   Print labeled array.
00827 */
00828 
00829 void sprint_array (char * str, int * array, int n)
00830 
00831 {
00832   int i;
00833 
00834   if (!quiet)
00835     {
00836       printf ("%s \n", str);
00837       print_array (array, n);
00838     }
00839 
00840   return;
00841 }
00842 
00843 
00844 /*---------------------------------------------------------------------------*/
00845 /*
00846   Write one time series array to specified file.
00847 */
00848 
00849 void write_one_ts (char * filename, int * array)
00850 {
00851   int i;
00852   FILE * outfile = NULL;
00853 
00854 
00855   outfile = fopen (filename, "w");
00856 
00857 
00858   for (i = 0;  i < NT;  i++)
00859     {
00860       fprintf (outfile, "%d", array[i]);
00861       fprintf (outfile, " \n");
00862     }
00863 
00864 
00865   fclose (outfile);
00866 }
00867 
00868 
00869 /*---------------------------------------------------------------------------*/
00870 /*
00871   Write multiple time series arrays to specified file.
00872 */
00873 
00874 void write_many_ts (char * filename, int * design)
00875 {
00876   int it, is;
00877   FILE * outfile = NULL;
00878 
00879 
00880   outfile = fopen (filename, "w");
00881 
00882 
00883   for (it = 0;  it < NT;  it++)
00884     {
00885       if (one_col)
00886         fprintf (outfile, "  %d", design[it]);
00887       else
00888         for (is = 0;  is < num_stimts;  is++)
00889           if (design[it] == is+1)
00890             fprintf (outfile, "  %d", 1);
00891           else
00892             fprintf (outfile, "  %d", 0);         
00893       fprintf (outfile, " \n");
00894     }
00895 
00896 
00897   fclose (outfile);
00898 }
00899 
00900 
00901 /*---------------------------------------------------------------------------*/
00902 /*
00903   Write experimental design to stimulus function time series files.
00904 */
00905 
00906 void write_results (char * prefix, int * design, int NT)
00907 {
00908   char filename[THD_MAX_NAME];    /* output file name */
00909   int * array;                    /* output binary sequence representing one
00910                                      stimulus function. */
00911   int is, i;
00912 
00913 
00914   if (one_file || one_col)
00915     {
00916       sprintf (filename, "%s.1D", prefix);
00917       if (!quiet)  printf ("\nWriting file: %s\n", filename);
00918       write_many_ts (filename, design);
00919     }
00920 
00921   else
00922     {
00923       /*----- Allocate memory for output array -----*/
00924       array = (int *) malloc (sizeof(int) * NT);
00925       MTEST (array);
00926 
00927       for (is = 1;  is <= num_stimts; is++)
00928         {
00929           sprintf (filename, "%s%d.1D", prefix, is);
00930           if (!quiet)  printf ("\nWriting file: %s\n", filename);
00931           for (i = 0;  i < NT;  i++)
00932             {
00933               if (design[i] == is)  array[i] = 1;
00934               else                  array[i] = 0;
00935             }
00936           write_one_ts (filename, array);
00937         }  
00938       
00939       /*----- Deallocate memory -----*/
00940       free (array);   array = NULL;
00941     }
00942 
00943 }
00944 
00945 
00946 /*---------------------------------------------------------------------------*/
00947 /*
00948   Write permutation of input table.
00949 */
00950 
00951 void write_table 
00952 (
00953   char * prefix,                  /* prefix for output file */
00954   int * design,                   /* permutation array */
00955   MRI_IMAGE * flim                /* data structure containing input table */
00956 )
00957 
00958 {
00959   FILE * outfile = NULL;          /* pointer to output file */
00960   char filename[THD_MAX_NAME];    /* output file name */
00961   int nx, ny;                     /* dimensions of input table */
00962   int it, is, icol;               /* row and column indices */
00963   float * far;                    /* pointer to input table float array */ 
00964 
00965 
00966   /*----- Assemble output file name -----*/
00967   sprintf (filename, "%s.1D", prefix);
00968   outfile = fopen (filename, "w");
00969 
00970   /*----- Get dimensions of input table -----*/
00971   nx = flim->nx ;
00972   ny = flim->ny ;
00973 
00974   /*----- Set pointer to input table float array -----*/
00975   far = MRI_FLOAT_PTR (flim);
00976 
00977   /*----- Loop over rows of table -----*/
00978   for (it = 0;  it <nx;  it++)
00979     {
00980       /*----- Get permuted row index -----*/
00981       is = design[it]-1;
00982 
00983       /*----- Copy row from input table to output file -----*/
00984       for (icol = 0;  icol < ny;  icol++)
00985         {
00986           fprintf (outfile, "  %f", far[icol*nx+is]);
00987         }
00988       fprintf (outfile, " \n");
00989     }
00990 
00991   fclose (outfile);
00992 }
00993 
00994 
00995 /*---------------------------------------------------------------------------*/
00996 
00997 int main 
00998 (
00999   int argc,                /* number of input arguments */
01000   char ** argv             /* array of input arguments */ 
01001 )
01002 
01003 {
01004   int * darray = NULL;     /* design array (block length = 1) */
01005   int * earray = NULL;     /* expanded array (arbitrary block length) */
01006   MRI_IMAGE * flim = NULL; /* data structure containing input table */
01007 
01008   
01009   /*----- Perform program initialization -----*/
01010   initialize (argc, argv, &darray, &earray, &flim);
01011 
01012 
01013   /*----- Use Markov chain to generate random stimulus functions -----*/
01014   if (markov)
01015     {
01016       markov_array (darray);
01017       sprint_array ("\nMarkov chain time series: ", darray, nt);
01018     }
01019 
01020   /*----- Use random permutations to generate random stimulus functions -----*/
01021   else
01022     {
01023       /*----- Generate required number of repetitions of stim. fns. -----*/
01024       fill_array (darray);
01025       if (!quiet)  sprint_array ("\nOriginal array: ", darray, nt);
01026       
01027       /*----- Permute the stimulus functions -----*/
01028       if (pseed)  
01029         {
01030           permute_array (darray);
01031           if (!quiet)  sprint_array ("\nPermuted array: ", darray, nt);
01032         }
01033 
01034       /*----- Randomize the order of the stimulus functions -----*/
01035       shuffle_array (darray);
01036       if (!quiet)  sprint_array ("\nShuffled array: ", darray, nt);
01037 
01038     }
01039 
01040   
01041   /*----- Expand the darray for block type designs -----*/
01042   expand_array (darray, earray);
01043   if (expand && (!quiet))  sprint_array ("\nExpanded array: ", earray, NT);
01044   
01045 
01046   /*----- Output results -----*/
01047   if (prefix != NULL)  
01048     {
01049       if (flim == NULL)
01050         write_results (prefix, earray, NT);
01051       else
01052         write_table (prefix, earray, flim);
01053     }
01054 
01055 
01056   /*----- Deallocate memory -----*/
01057   if (darray != NULL)  { free (darray);   darray = NULL; }
01058   if (earray != NULL)  { free (earray);   earray = NULL; }
01059   if (flim   != NULL)  { mri_free(flim);  flim = NULL; }
01060   
01061   exit(0);
01062 }
01063 
01064 /*---------------------------------------------------------------------------*/
01065 
01066 
01067 
01068