1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1999-2003, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 /*---------------------------------------------------------------------------*/
8 /*
9 This sample program generates random stimulus functions.
10
11 File: RSFgen.c
12 Author: B. Douglas Ward
13 Date: 06 July 1999
14
15
16 Mod: Increased max. allowed number of input stimulus functions.
17 Date: 24 August 1999
18
19 Mod: Added option to create multiple column stimulus function files.
20 Date: 24 November 1999
21
22 Mod: Changed option label "-nstim" to "-num_stimts".
23 Date: 29 November 1999
24
25 Mod: Added option "-nblock" to specify block length for each stim fn.
26 Date: 15 December 1999
27
28 Mod: Added flag to expand array for block type design.
29 Date: 14 January 2000
30
31 Mod: Added -markov and -pzero options.
32 Date: 03 October 2000
33
34 Mod: Added -pseed option for permutation of stimulus labels.
35 Date: 27 April 2001
36
37 Mod: Changes to eliminate constraints on number of stimulus time series.
38 Date: 11 May 2001
39
40 Mod: Added call to AFNI_logger.
41 Date: 15 August 2001
42
43 Mod: Made -nblock option compatible with the Markov Chain options.
44 Date: 06 March 2002
45
46 Mod: Added -one_col option to write stimulus functions as a single
47 column of decimal integers.
48 Date: 18 April 2002
49
50 Mod: Added -quiet option.
51 Date: 30 December 2002
52
53 Mod: Added -table option, to generate random permutations of the rows
54 of an input column or table of numbers.
55 Date: 13 March 2003
56
57 */
58
59 /*---------------------------------------------------------------------------*/
60
61 #define PROGRAM_NAME "RSFgen" /* name of this program */
62 #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */
63 #define PROGRAM_INITIAL "06 July 1999" /* date of initial program release */
64 #define PROGRAM_LATEST "13 March 2003" /* date of latest program revision */
65
66 /*---------------------------------------------------------------------------*/
67
68 #include <math.h>
69 #include <stdlib.h>
70 #include <stdio.h>
71
72 #include "mrilib.h"
73 #include "matrix.h"
74
75 #include "randgen.c"
76
77 /*---------------------------------------------------------------------------*/
78 /*
79 Global variables.
80 */
81
82
83 int NT = 0; /* length of stimulus time series */
84 int nt = 0; /* length of simple time series (block length = 1) */
85 int num_stimts = 0; /* number of input stimuli (experimental conditions) */
86 int * num_reps = NULL; /* number of repetitions for each stimulus */
87 int * nblock = NULL; /* block length for each stimulus */
88 int expand = 0; /* flag to expand the array for block type design */
89 long seed = 1234567; /* random number seed */
90 long pseed = 0; /* stimulus permutation random number seed */
91 char * prefix = NULL; /* prefix for output .1D stimulus functions */
92 int one_file = 0; /* flag for place stim functions into a single file */
93 int one_col = 0; /* flag for write stim functions as a single column */
94 int markov = 0; /* flag for Markov process */
95 char * tpm_file = NULL; /* file name for input of transition prob. matrix */
96 float pzero = 0.0; /* zero (null) state probability */
97 int quiet = 0; /* flag for suppress screen output */
98 char * table_file = NULL; /* input time series file */
99
100
101 /*---------------------------------------------------------------------------*/
102 /*
103 Print error message and stop.
104 */
105
RSF_error(char * message)106 void RSF_error (char * message)
107 {
108 fprintf (stderr, "\n%s Error: %s \n", PROGRAM_NAME, message);
109 exit(1);
110 }
111
112
113 /*---------------------------------------------------------------------------*/
114
115 /** macro to test a malloc-ed pointer for validity **/
116
117 #define MTEST(ptr) \
118 if((ptr)==NULL) \
119 ( RSF_error ("Cannot allocate memory") )
120
121
122 /*---------------------------------------------------------------------------*/
123 /*
124 Identify software.
125 */
126
identify_software()127 void identify_software ()
128 {
129
130 #if 0
131 /*----- Identify software -----*/
132 printf ("\n\n");
133 printf ("Program: %s \n", PROGRAM_NAME);
134 printf ("Author: %s \n", PROGRAM_AUTHOR);
135 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
136 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
137 printf ("\n");
138 #endif
139
140 PRINT_VERSION("RSFgen") ; AUTHOR(PROGRAM_AUTHOR) ;
141 }
142
143
144 /*---------------------------------------------------------------------------*/
145 /*
146 Routine to display help menu for program RSFgen.
147 */
148
display_help_menu()149 void display_help_menu()
150 {
151 identify_software();
152
153 printf (
154 "Sample program to generate random stimulus functions. \n"
155 " \n"
156 "Usage: \n"
157 PROGRAM_NAME " \n"
158 "-nt n n = length of time series \n"
159 "-num_stimts p p = number of input stimuli (experimental conditions) \n"
160 "[-nblock i k] k = block length for stimulus i (1<=i<=p) \n"
161 " (default: k = 1) \n"
162 "[-seed s] s = random number seed \n"
163 "[-quiet] flag to suppress screen output \n"
164 "[-one_file] place stimulus functions into a single .1D file \n"
165 "[-one_col] write stimulus functions as a single column of decimal\n"
166 " integers (default: multiple columns of binary nos.) \n"
167 "[-prefix pname] pname = prefix for p output .1D stimulus functions \n"
168 " e.g., pname1.1D, pname2.1D, ..., pnamep.1D \n"
169 " \n"
170 "The following Random Permutation, Markov Chain, and Input Table options\n"
171 "are mutually exclusive. \n"
172 " \n"
173 "Random Permutation options: \n"
174 "-nreps i r r = number of repetitions for stimulus i (1<=i<=p) \n"
175 "[-pseed s] s = stim label permutation random number seed \n"
176 " p \n"
177 " Note: Require n >= Sum (r[i] * k[i]) \n"
178 " i=1 \n"
179 " \n"
180 "Markov Chain options: \n"
181 "-markov mfile mfile = file containing the transition prob. matrix \n"
182 "[-pzero z] probability of a zero (i.e., null) state \n"
183 " (default: z = 0) \n"
184 " \n"
185 "Input Table row permutation options: \n"
186 "[-table dfile] dfile = filename of column or table of numbers \n"
187 " Note: dfile may have a column selector attached \n"
188 " Note: With this option, all other input options, \n"
189 " except -seed and -prefix, are ignored \n"
190 " \n"
191 " \n"
192 "Warning: This program will overwrite pre-existing .1D files \n"
193 " \n"
194 );
195
196 exit(0);
197 }
198
199
200 /*---------------------------------------------------------------------------*/
201 /*
202 Routine to get user specified input options.
203 */
204
get_options(int argc,char ** argv)205 void get_options
206 (
207 int argc, /* number of input arguments */
208 char ** argv /* array of input arguments */
209 )
210
211 {
212 int nopt = 1; /* input option argument counter */
213 int ival; /* integer input */
214 float fval; /* float input */
215 long lval; /* long integer input */
216 char message[THD_MAX_NAME]; /* error message */
217 int i;
218
219
220 /*----- Does user request help menu? -----*/
221 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
222
223
224 /*----- add to program log -----*/
225 AFNI_logger (PROGRAM_NAME,argc,argv);
226
227
228 /*----- Main loop over input options -----*/
229 while (nopt < argc )
230 {
231
232 /*----- -nt n -----*/
233 if (strncmp(argv[nopt], "-nt", 3) == 0)
234 {
235 nopt++;
236 if (nopt >= argc) RSF_error ("need argument after -nt ");
237 sscanf (argv[nopt], "%d", &ival);
238 if (ival <= 0)
239 RSF_error ("illegal argument after -nt ");
240 NT = ival;
241 nopt++;
242 continue;
243 }
244
245
246 /*----- -num_stimts p -----*/
247 if (strncmp(argv[nopt], "-num_stimts", 11) == 0)
248 {
249 nopt++;
250 if (nopt >= argc) RSF_error ("need argument after -num_stimts ");
251 sscanf (argv[nopt], "%d", &ival);
252 if (ival <= 0)
253 RSF_error ("illegal argument after -num_stimts ");
254 num_stimts = ival;
255
256 /*----- Initialize repetition number array -----*/
257 num_reps = (int *) malloc (sizeof(int) * num_stimts);
258 MTEST (num_reps);
259 for (i = 0; i < num_stimts; i++)
260 num_reps[i] = 0;
261
262 /*----- Initialize block length array -----*/
263 nblock = (int *) malloc (sizeof(int) * num_stimts);
264 MTEST (nblock);
265 for (i = 0; i < num_stimts; i++)
266 nblock[i] = 1;
267
268 nopt++;
269 continue;
270 }
271
272
273 /*----- -nreps i r -----*/
274 if (strncmp(argv[nopt], "-nreps", 6) == 0)
275 {
276 nopt++;
277 if (nopt+1 >= argc) RSF_error ("need 2 arguments after -nreps ");
278 sscanf (argv[nopt], "%d", &ival);
279 if ((ival <= 0) || (ival > num_stimts))
280 RSF_error ("illegal i argument for -nreps i r ");
281 i = ival - 1;
282 nopt++;
283
284 sscanf (argv[nopt], "%d", &ival);
285 if (ival <= 0)
286 RSF_error ("illegal r argument for -nreps i r ");
287 num_reps[i] = ival;
288 nopt++;
289 continue;
290 }
291
292
293 /*----- -nblock i k -----*/
294 if (strncmp(argv[nopt], "-nblock", 7) == 0)
295 {
296 nopt++;
297 if (nopt+1 >= argc) RSF_error ("need 2 arguments after -nblock ");
298 sscanf (argv[nopt], "%d", &ival);
299 if ((ival <= 0) || (ival > num_stimts))
300 RSF_error ("illegal i argument for -nblock i k ");
301 i = ival - 1;
302 nopt++;
303
304 sscanf (argv[nopt], "%d", &ival);
305 if (ival <= 0)
306 RSF_error ("illegal k argument for -nblock i k ");
307 nblock[i] = ival;
308 if (ival > 1) expand = 1;
309 nopt++;
310 continue;
311 }
312
313
314 /*----- -seed s -----*/
315 if (strncmp(argv[nopt], "-seed", 5) == 0)
316 {
317 nopt++;
318 if (nopt >= argc) RSF_error ("need argument after -seed ");
319 sscanf (argv[nopt], "%ld", &lval);
320 if (lval <= 0)
321 RSF_error ("illegal argument after -seed ");
322 seed = lval;
323 nopt++;
324 continue;
325 }
326
327
328 /*----- -pseed s -----*/
329 if (strcmp(argv[nopt], "-pseed") == 0)
330 {
331 nopt++;
332 if (nopt >= argc) RSF_error ("need argument after -pseed ");
333 sscanf (argv[nopt], "%ld", &lval);
334 if (lval <= 0)
335 RSF_error ("illegal argument after -pseed ");
336 pseed = lval;
337 nopt++;
338 continue;
339 }
340
341
342 /*----- -quiet -----*/
343 if (strcmp(argv[nopt], "-quiet") == 0)
344 {
345 quiet = 1;
346 nopt++;
347 continue;
348 }
349
350
351 /*----- -one_file -----*/
352 if (strncmp(argv[nopt], "-one_file", 9) == 0)
353 {
354 one_file = 1;
355 nopt++;
356 continue;
357 }
358
359
360 /*----- -one_col -----*/
361 if (strcmp(argv[nopt], "-one_col") == 0)
362 {
363 one_col = 1;
364 nopt++;
365 continue;
366 }
367
368
369 /*----- -prefix pname -----*/
370 if (strncmp(argv[nopt], "-prefix", 7) == 0)
371 {
372 nopt++;
373 if (nopt >= argc) RSF_error ("need argument after -prefix ");
374 prefix = AFMALL(char, sizeof(char) * THD_MAX_NAME);
375 MTEST (prefix);
376 strcpy (prefix, argv[nopt]);
377 nopt++;
378 continue;
379 }
380
381
382 /*----- -markov mfile -----*/
383 if (strcmp(argv[nopt], "-markov") == 0)
384 {
385 markov = 1;
386 nopt++;
387 if (nopt >= argc) RSF_error ("need argument after -markov ");
388 tpm_file = AFMALL(char, sizeof(char) * THD_MAX_NAME);
389 MTEST (tpm_file);
390 strcpy (tpm_file, argv[nopt]);
391 nopt++;
392 continue;
393 }
394
395
396 /*----- -pzero z -----*/
397 if (strcmp(argv[nopt], "-pzero") == 0)
398 {
399 nopt++;
400 if (nopt >= argc) RSF_error ("need argument after -pzero ");
401 sscanf (argv[nopt], "%f", &fval);
402 if ((fval < 0.0) || (fval > 1.0))
403 RSF_error ("Require 0.0 <= pzero <= 1.0");
404 pzero = fval;
405 nopt++;
406 continue;
407 }
408
409
410 /*----- -table dfile -----*/
411 if (strcmp(argv[nopt], "-table") == 0)
412 {
413 nopt++;
414 if (nopt >= argc) RSF_error ("need argument after -table ");
415 table_file = AFMALL(char, sizeof(char) * THD_MAX_NAME);
416 MTEST (table_file);
417 strcpy (table_file, argv[nopt]);
418 nopt++;
419 continue;
420 }
421
422
423 /*----- Unknown command -----*/
424 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
425 RSF_error (message);
426
427 }
428 }
429
430
431 /*---------------------------------------------------------------------------*/
432 /*
433 Read input table.
434 */
435
read_table(char * table_file,MRI_IMAGE ** flim)436 void read_table
437 (
438 char * table_file, /* file containing input table */
439 MRI_IMAGE ** flim /* data structure containing input table */
440 )
441
442 {
443 int i;
444 char message[THD_MAX_NAME]; /* error message */
445
446 /*----- Read table -----*/
447 *flim = mri_read_1D(table_file);
448 if ((*flim) == NULL)
449 {
450 sprintf (message, "Unable to read table file: %s", table_file);
451 RSF_error (message);
452 }
453
454 /*----- Initialize control variables -----*/
455 NT = (*flim)->nx;
456 num_stimts = NT;
457 markov = 0;
458 pseed = 0;
459
460 /*----- Initialize repetition number array -----*/
461 num_reps = (int *) malloc (sizeof(int) * num_stimts);
462 MTEST (num_reps);
463 for (i = 0; i < num_stimts; i++)
464 num_reps[i] = 1;
465
466 /*----- Initialize block length array -----*/
467 nblock = (int *) malloc (sizeof(int) * num_stimts);
468 MTEST (nblock);
469 for (i = 0; i < num_stimts; i++)
470 nblock[i] = 1;
471
472 }
473
474
475 /*---------------------------------------------------------------------------*/
476 /*
477 Print input options.
478 */
479
print_options()480 void print_options ()
481
482 {
483 int i;
484
485 identify_software();
486
487 if (table_file != NULL)
488 {
489 printf ("table file = %s \n", table_file);
490 printf ("nt = %d \n", NT);
491 printf ("seed = %ld \n", seed);
492 printf ("output prefix = %s \n", prefix);
493 }
494 else
495 {
496 printf ("nt = %d \n", NT);
497 printf ("num_stimts = %d \n", num_stimts);
498 printf ("seed = %ld \n", seed);
499 if (pseed) printf ("pseed = %ld \n", pseed);
500 printf ("output prefix = %s \n", prefix);
501 if (markov)
502 {
503 printf ("TPM file = %s \n", tpm_file);
504 printf ("pzero = %f \n", pzero);
505 for (i = 0; i < num_stimts; i++)
506 printf ("nblock[%d] = %d \n", i+1, nblock[i]);
507 }
508 else
509 for (i = 0; i < num_stimts; i++)
510 printf ("nreps[%d] = %d nblock[%d] = %d \n",
511 i+1, num_reps[i], i+1, nblock[i]);
512 }
513 }
514
515
516 /*---------------------------------------------------------------------------*/
517 /*
518 Perform all program initialization.
519 */
520
initialize(int argc,char ** argv,int ** darray,int ** earray,MRI_IMAGE ** flim)521 void initialize
522 (
523 int argc, /* number of input arguments */
524 char ** argv, /* array of input arguments */
525 int ** darray, /* original design array (block length = 1) */
526 int ** earray, /* expanded design array (arbitrary block length) */
527 MRI_IMAGE ** flim /* data structure containing input table */
528 )
529
530 {
531 int i, total;
532
533
534 /*----- Get command line inputs -----*/
535 get_options (argc, argv);
536
537
538 /*----- Read input table -----*/
539 if (table_file != NULL) read_table (table_file, flim);
540
541
542 /*----- Print input options -----*/
543 if (! quiet) print_options ();
544
545
546 /*----- Check for valid inputs -----*/
547 if (NT == 0) RSF_error ("Must specify nt");
548 if (num_stimts == 0) RSF_error ("Must specify num_stimts");
549 total = 0;
550 nt = NT;
551
552 if (! markov)
553 {
554 for (i = 0; i < num_stimts; i++)
555 {
556 if (num_reps[i] == 0)
557 RSF_error ("Must specify nreps > 0 for each stimulus");
558 total += num_reps[i] * nblock[i];
559 nt -= num_reps[i] * (nblock[i] - 1);
560 }
561 if (total > NT) RSF_error ("Require nt >= Sum (r[i] * k[i]) ");
562 }
563
564
565 /*----- Allocate memory for experimental design -----*/
566 *darray = (int *) malloc (sizeof(int) * nt); MTEST (*darray);
567 *earray = (int *) malloc (sizeof(int) * NT); MTEST (*earray);
568
569
570 }
571
572
573 /*---------------------------------------------------------------------------*/
574 /*
575 Use Markov chain to create random stimulus functions.
576 */
577
markov_array(int * design)578 void markov_array (int * design)
579
580 {
581 int it, is, id, isprev;
582 float prob, cumprob;
583 matrix tpm;
584 char message[THD_MAX_NAME]; /* error message */
585
586
587 matrix_initialize (&tpm);
588
589
590 /*----- Read the transition probability matrix -----*/
591 matrix_file_read (tpm_file, num_stimts, num_stimts, &tpm, 1);
592 if (tpm.elts == NULL)
593 {
594 sprintf (message, "Unable to read Markov chain matrix from file: %s",
595 tpm_file);
596 RSF_error (message);
597 }
598 if (!quiet) matrix_sprint ("\nTPM matrix:", tpm);
599
600
601 /*----- Verify that the TPM has the correct form -----*/
602 for (is = 0; is < num_stimts; is++)
603 {
604 cumprob = 0.0;
605 for (it = 0; it < num_stimts; it++)
606 cumprob += tpm.elts[is][it];
607 if (cumprob < 0.9999)
608 {
609 sprintf (message, "Row %d of TPM sums to %f, which is < 1.0",
610 is, cumprob);
611 RSF_error (message);
612 }
613 if (cumprob > 1.0001)
614 {
615 sprintf (message, "Row %d of TPM sums to %f, which is > 1.0",
616 is, cumprob);
617 RSF_error (message);
618 }
619 }
620
621
622 /*----- Initialize the experimental design array -----*/
623 for (it = 0; it < NT; it++)
624 design[it] = 0;
625
626
627 /*----- Initialize random number seed -----*/
628 srand48 (seed);
629
630
631 /*----- Generate Markov process -----*/
632 isprev = (int) (rand_uniform(0.0,1.0)*num_stimts);
633 it = 0; id = 0;
634 while (it < NT)
635 {
636 if ((pzero > 0.0) && (rand_uniform(0.0,1.0) < pzero))
637 {
638 design[id] = 0;
639 id++; it++;
640 }
641 else
642 {
643 prob = rand_uniform(0.0,1.0);
644 cumprob = 0.0;
645 for (is = 0; is < num_stimts; is++)
646 {
647 cumprob += tpm.elts[isprev][is];
648 if (prob <= cumprob)
649 {
650 design[id] = is+1;
651 isprev = is;
652 id++; it += nblock[is];
653 break;
654 }
655 }
656 }
657 }
658
659 nt = id;
660
661
662 matrix_destroy (&tpm);
663
664 return;
665 }
666
667
668 /*---------------------------------------------------------------------------*/
669 /*
670 Fill the experimental design array.
671 */
672
fill_array(int * design)673 void fill_array (int * design)
674
675 {
676 int i, is, m;
677
678
679 for (i = 0; i < nt; i++)
680 design[i] = 0;
681
682 i = 0;
683 for (is = 0; is < num_stimts; is++)
684 {
685 for (m = 0; m < num_reps[is]; m++)
686 {
687 design[i] = is+1;
688 i++;
689 }
690 }
691
692
693 return;
694 }
695
696
697 /*---------------------------------------------------------------------------*/
698 /*
699 Permute the stimulus functions within the experimental design array.
700 */
701
permute_array(int * design)702 void permute_array (int * design)
703
704 {
705 int i, j, temp;
706 int is, nb;
707
708
709 /*----- Initialize random number seed -----*/
710 srand48 (pseed);
711
712
713 /*----- Determine total number of blocks -----*/
714 nb = 0;
715 for (is = 0; is < num_stimts; is++)
716 nb += num_reps[is];
717
718
719 /*----- Permute the blocks -----*/
720 for (i = 0; i < nb; i++)
721 {
722 j = rand_uniform(0.0,1.0) * nb;
723
724 /*----- Just in case -----*/
725 if (j < 0) j = 0;
726 if (j > nb-1) j = nb-1;
727
728 temp = design[i];
729 design[i] = design[j];
730 design[j] = temp;
731 }
732
733 return;
734 }
735
736
737 /*---------------------------------------------------------------------------*/
738 /*
739 Shuffle the experimental design array.
740 */
741
shuffle_array(int * design)742 void shuffle_array (int * design)
743
744 {
745 int i, j, temp;
746
747
748 /*----- Initialize random number seed -----*/
749 srand48 (seed);
750
751
752 for (i = 0; i < nt; i++)
753 {
754 j = rand_uniform(0.0,1.0) * nt;
755
756 /*----- Just in case -----*/
757 if (j < 0) j = 0;
758 if (j > nt-1) j = nt-1;
759
760 temp = design[i];
761 design[i] = design[j];
762 design[j] = temp;
763 }
764
765 return;
766 }
767
768
769 /*---------------------------------------------------------------------------*/
770 /*
771 Expand the experimental design array, to allow for block-type designs.
772 */
773
expand_array(int * darray,int * earray)774 void expand_array (int * darray, int * earray)
775
776 {
777 int i, j, k, m;
778
779 j = 0;
780 for (i = 0; i < nt, j < NT; i++)
781 {
782 m = darray[i];
783
784 if (m == 0)
785 {
786 earray[j] = 0;
787 j++;
788 }
789 else
790 {
791 for (k = 0; k < nblock[m-1]; k++)
792 {
793 earray[j] = m;
794 j++;
795 if (j >= NT) break;
796 }
797 }
798 }
799
800 return;
801 }
802
803
804 /*---------------------------------------------------------------------------*/
805 /*
806 Print array.
807 */
808
print_array(int * array,int n)809 void print_array (int * array, int n)
810
811 {
812 int i;
813
814 for (i = 0; i < n; i++)
815 {
816 printf (" %2d ", array[i]);
817 if ((i+1) % 20 == 0) printf ("\n");
818 }
819
820 printf ("\n");
821
822 return;
823 }
824
825
826 /*---------------------------------------------------------------------------*/
827 /*
828 Print labeled array.
829 */
830
sprint_array(char * str,int * array,int n)831 void sprint_array (char * str, int * array, int n)
832
833 {
834 int i;
835
836 if (!quiet)
837 {
838 printf ("%s \n", str);
839 print_array (array, n);
840 }
841
842 return;
843 }
844
845
846 /*---------------------------------------------------------------------------*/
847 /*
848 Write one time series array to specified file.
849 */
850
write_one_ts(char * filename,int * array)851 void write_one_ts (char * filename, int * array)
852 {
853 int i;
854 FILE * outfile = NULL;
855
856
857 outfile = fopen (filename, "w");
858
859
860 for (i = 0; i < NT; i++)
861 {
862 fprintf (outfile, "%d", array[i]);
863 fprintf (outfile, " \n");
864 }
865
866
867 fclose (outfile);
868 }
869
870
871 /*---------------------------------------------------------------------------*/
872 /*
873 Write multiple time series arrays to specified file.
874 */
875
write_many_ts(char * filename,int * design)876 void write_many_ts (char * filename, int * design)
877 {
878 int it, is;
879 FILE * outfile = NULL;
880
881
882 outfile = fopen (filename, "w");
883
884
885 for (it = 0; it < NT; it++)
886 {
887 if (one_col)
888 fprintf (outfile, " %d", design[it]);
889 else
890 for (is = 0; is < num_stimts; is++)
891 if (design[it] == is+1)
892 fprintf (outfile, " %d", 1);
893 else
894 fprintf (outfile, " %d", 0);
895 fprintf (outfile, " \n");
896 }
897
898
899 fclose (outfile);
900 }
901
902
903 /*---------------------------------------------------------------------------*/
904 /*
905 Write experimental design to stimulus function time series files.
906 */
907
write_results(char * prefix,int * design,int NT)908 void write_results (char * prefix, int * design, int NT)
909 {
910 char filename[THD_MAX_NAME]; /* output file name */
911 int * array; /* output binary sequence representing one
912 stimulus function. */
913 int is, i;
914
915
916 if (one_file || one_col)
917 {
918 sprintf (filename, "%s.1D", prefix);
919 if (!quiet) printf ("\nWriting file: %s\n", filename);
920 write_many_ts (filename, design);
921 }
922
923 else
924 {
925 /*----- Allocate memory for output array -----*/
926 array = (int *) malloc (sizeof(int) * NT);
927 MTEST (array);
928
929 for (is = 1; is <= num_stimts; is++)
930 {
931 sprintf (filename, "%s%d.1D", prefix, is);
932 if (!quiet) printf ("\nWriting file: %s\n", filename);
933 for (i = 0; i < NT; i++)
934 {
935 if (design[i] == is) array[i] = 1;
936 else array[i] = 0;
937 }
938 write_one_ts (filename, array);
939 }
940
941 /*----- Deallocate memory -----*/
942 free (array); array = NULL;
943 }
944
945 }
946
947
948 /*---------------------------------------------------------------------------*/
949 /*
950 Write permutation of input table.
951 */
952
write_table(char * prefix,int * design,MRI_IMAGE * flim)953 void write_table
954 (
955 char * prefix, /* prefix for output file */
956 int * design, /* permutation array */
957 MRI_IMAGE * flim /* data structure containing input table */
958 )
959
960 {
961 FILE * outfile = NULL; /* pointer to output file */
962 char filename[THD_MAX_NAME]; /* output file name */
963 int nx, ny; /* dimensions of input table */
964 int it, is, icol; /* row and column indices */
965 float * far; /* pointer to input table float array */
966
967
968 /*----- Assemble output file name -----*/
969 sprintf (filename, "%s.1D", prefix);
970 outfile = fopen (filename, "w");
971
972 /*----- Get dimensions of input table -----*/
973 nx = flim->nx ;
974 ny = flim->ny ;
975
976 /*----- Set pointer to input table float array -----*/
977 far = MRI_FLOAT_PTR (flim);
978
979 /*----- Loop over rows of table -----*/
980 for (it = 0; it <nx; it++)
981 {
982 /*----- Get permuted row index -----*/
983 is = design[it]-1;
984
985 /*----- Copy row from input table to output file -----*/
986 for (icol = 0; icol < ny; icol++)
987 {
988 fprintf (outfile, " %f", far[icol*nx+is]);
989 }
990 fprintf (outfile, " \n");
991 }
992
993 fclose (outfile);
994 }
995
996
997 /*---------------------------------------------------------------------------*/
998
main(int argc,char ** argv)999 int main
1000 (
1001 int argc, /* number of input arguments */
1002 char ** argv /* array of input arguments */
1003 )
1004
1005 {
1006 int * darray = NULL; /* design array (block length = 1) */
1007 int * earray = NULL; /* expanded array (arbitrary block length) */
1008 MRI_IMAGE * flim = NULL; /* data structure containing input table */
1009
1010
1011 machdep() ; mainENTRY("RSFgen") ;
1012
1013 /*----- Perform program initialization -----*/
1014 initialize (argc, argv, &darray, &earray, &flim);
1015
1016
1017 /*----- Use Markov chain to generate random stimulus functions -----*/
1018 if (markov)
1019 {
1020 markov_array (darray);
1021 sprint_array ("\nMarkov chain time series: ", darray, nt);
1022 }
1023
1024 /*----- Use random permutations to generate random stimulus functions -----*/
1025 else
1026 {
1027 /*----- Generate required number of repetitions of stim. fns. -----*/
1028 fill_array (darray);
1029 if (!quiet) sprint_array ("\nOriginal array: ", darray, nt);
1030
1031 /*----- Permute the stimulus functions -----*/
1032 if (pseed)
1033 {
1034 permute_array (darray);
1035 if (!quiet) sprint_array ("\nPermuted array: ", darray, nt);
1036 }
1037
1038 /*----- Randomize the order of the stimulus functions -----*/
1039 shuffle_array (darray);
1040 if (!quiet) sprint_array ("\nShuffled array: ", darray, nt);
1041
1042 }
1043
1044
1045 /*----- Expand the darray for block type designs -----*/
1046 expand_array (darray, earray);
1047 if (expand && (!quiet)) sprint_array ("\nExpanded array: ", earray, NT);
1048
1049
1050 /*----- Output results -----*/
1051 if (prefix != NULL)
1052 {
1053 if (flim == NULL)
1054 write_results (prefix, earray, NT);
1055 else
1056 write_table (prefix, earray, flim);
1057 }
1058
1059
1060 /*----- Deallocate memory -----*/
1061 if (darray != NULL) { free (darray); darray = NULL; }
1062 if (earray != NULL) { free (earray); earray = NULL; }
1063 if (flim != NULL) { mri_free(flim); flim = NULL; }
1064
1065 exit(0);
1066 }
1067
1068 /*---------------------------------------------------------------------------*/
1069
1070
1071
1072
1073