1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2001, 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 program performs multiple linear regression analysis across AFNI
10 3d datasets.
11
12 File: 3dRegAna.c
13 Author: B. Douglas Ward
14 Date: 10 October 1997
15
16
17 Mod: Changes to allow output of AFNI "bucket" type dataset.
18 Date: 17 December 1997
19
20 Mod: Converted some matrix data structures to vector data structures.
21 Date: 18 December 1997
22
23 Mod: Modified calculation of required disk space to account for the
24 'bucket' dataset.
25 Date: 08 January 1998
26
27 Mod: Some of the regression analysis software has been moved to
28 the include file RegAna.c.
29 Date: 20 August 1998
30
31 Mod: Restructured matrix calculations to improve execution speed.
32 Date: 16 December 1998
33
34 Mod: Add use of the argument list extenstion routine addto_args
35 to allow the last switch '-@' to get further command line
36 arguments from stdin
37 Date: 22 February 1999 -- RWCox
38
39 Mod: Added changes for incorporating History notes.
40 Date: 09 September 1999
41
42 Mod: Replaced dataset input code with calls to THD_open_dataset,
43 to allow operator selection of individual sub-bricks for input.
44 Date: 03 December 1999
45
46 Mod: Added call to AFNI_logger.
47 Date: 15 August 2001
48
49 Mod: Modified routine write_afni_data so that all output
50 subbricks will now have the scaled short integer format.
51 Date: 14 March 2002
52
53 Mod: Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
54 Date: 02 December 2002
55
56 Mod: Check for proper ':' usage in -model parameters.
57 Date: 24 August 2005 [rickr]
58
59 Mod: Changed default workmem to 750 MB (from 12).
60 Date: 08 December 2005 [rickr]
61 */
62
63 /*---------------------------------------------------------------------------*/
64
65 #define PROGRAM_NAME "3dRegAna" /* name of this program */
66 #define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */
67 #define PROGRAM_INITIAL "10 Oct 1997" /* date of initial program release */
68 #define PROGRAM_LATEST "24 Aug 2005" /* date of latest program revision */
69
70 /*---------------------------------------------------------------------------*/
71
72 #if 0
73 #define SUFFIX ".3dregana" /* suffix for temporary files */
74 #else
75 char SUFFIX[1024] = ".3dregana." ;
76 #endif
77
78 #include <stdio.h>
79 #include <math.h>
80 #include "mrilib.h"
81 #include "matrix.h"
82
83
84 #include "RegAna.c"
85
86 /*---------------------------------------------------------------------------*/
87
88 /*** HP-UX ***/
89 #ifdef HP
90 # define DF "bdf ."
91 #endif
92
93 /*** DEC OSF1 ***/
94 #ifdef OSF1
95 # define DF "df -k ."
96 #endif
97
98 /*** SGI IRIX ***/
99 #ifdef SGI
100 # define DF "df -k ."
101 #endif
102
103 /*** SunOS or Solaris ***/
104 #if defined(SOLARIS) || defined(SUN)
105 # define DF "df -k"
106 #endif
107
108 /*** IBM RS6000 ***/
109 #ifdef RS6000
110 #endif
111
112 /*** Linux 1.2.x ***/
113 #ifdef LINUX
114 # define DF "df -k ."
115 #endif
116
117 /*** other? ***/
118 #ifndef DF
119 # define DF "df"
120 #endif
121
122
123 /*---------------------------------------------------------------------------*/
124
125 #define MAX_XVARS 101 /* max. number of independent variables */
126 #define MAX_OBSERVATIONS 1000 /* max. number of input datasets */
127 #define MAX_NAME_LENGTH THD_MAX_NAME /* max. string length for file names */
128 #define MEGA 1048576 /* one megabyte */
129
130 static char * commandline = NULL ; /* command line for history notes */
131
132
133 typedef struct model
134 {
135 int p; /* number of parameters in full model */
136 int flist[MAX_XVARS]; /* list of parameters in full model */
137 int q; /* number of parameters in reduced model */
138 int rlist[MAX_XVARS]; /* list of parameters in reduced model */
139 } model;
140
141
142 typedef struct RA_options
143 {
144 int datum; /* data type for "intensity" data subbrick */
145 char session[MAX_NAME_LENGTH]; /* name of output directory */
146
147 char ** yname; /* names of the input data files */
148 char * first_dataset; /* name of the first data set */
149
150 int nx, ny, nz; /* data set dimensions */
151 int nxyz; /* number of voxels per image */
152
153 int diskspace; /* if true, print out how much disk space
154 is required for program execution */
155 int workmem; /* working memory */
156 int piece_size; /* number of voxels in dataset piece */
157 int num_pieces; /* dataset is divided into this many pieces */
158
159 float rms_min; /* minimum rms error to reject constant model */
160 float fdisp; /* minimum F-value for screen output */
161
162 int * levels; /* indices for repeat observations */
163 int * counts; /* counts of repeat observations */
164 int c; /* number of unique rows in ind. var. matrix */
165 float flofmax; /* max. allowed F due to lack of fit */
166
167 int numf; /* number of F-stats volumes */
168 int numr; /* number of R^2 volumes */
169 int numt; /* number of t-stats volumes */
170 int numc; /* number of reg. coef. volumes */
171
172 char ** fcoef_filename; /* file names for reg. coefs. and F-stats */
173 char ** rcoef_filename; /* file names for reg. coefs. and R^2 */
174 char ** tcoef_filename; /* file names for reg. coefs. and t-stats */
175
176 int numfiles; /* number of 2 sub-brick output files */
177
178 char * bucket_filename; /* file name for bucket dataset */
179 int numbricks; /* number of sub-bricks in bucket dataset */
180 int * brick_type; /* indicates type of sub-brick */
181 int * brick_coef; /* regression coefficient number for sub-brick*/
182 char ** brick_label; /* character string label for sub-brick */
183
184 } RA_options;
185
186
187 /*---------------------------------------------------------------------------*/
188 /*
189 Print error message and stop.
190 */
191
RA_error(char * message)192 void RA_error (char * message)
193 {
194 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
195 exit(1);
196 }
197
198
199 /*---------------------------------------------------------------------------*/
200 /*
201 Routine to display 3dRegAna help menu.
202 */
203
display_help_menu()204 void display_help_menu()
205 {
206 printf
207 (
208 "This program performs multiple linear regression analysis. \n\n"
209 "Usage: \n"
210 "3dRegAna \n"
211 "-rows n number of input datasets \n"
212 "-cols m number of X variables \n"
213 "-xydata X11 X12 ... X1m filename X variables and Y observations \n"
214 " . . \n"
215 " . . \n"
216 " . . \n"
217 "-xydata Xn1 Xn2 ... Xnm filename X variables and Y observations \n"
218 " \n"
219 "-model i1 ... iq : j1 ... jr definition of linear regression model; \n"
220 " reduced model: \n"
221 " Y = f(Xj1,...,Xjr) \n"
222 " full model: \n"
223 " Y = f(Xj1,...,Xjr,Xi1,...,Xiq) \n"
224 " \n"
225 "[-diskspace] print out disk space required for program execution\n"
226 "[-workmem mega] number of megabytes of RAM to use for statistical \n"
227 " workspace (default = 750 (was 12)) \n"
228 "[-rmsmin r] r = minimum rms error to reject constant model \n"
229 "[-fdisp fval] display (to screen) results for those voxels \n"
230 " whose F-statistic is > fval \n"
231 " \n"
232 "[-flof alpha] alpha = minimum p value for F due to lack of fit \n"
233 " \n"
234 " \n"
235 "The following commands generate individual AFNI 2 sub-brick datasets: \n"
236 " \n"
237 "[-fcoef k prefixname] estimate of kth regression coefficient \n"
238 " along with F-test for the regression \n"
239 " is written to AFNI `fift' dataset \n"
240 "[-rcoef k prefixname] estimate of kth regression coefficient \n"
241 " along with coef. of mult. deter. R^2 \n"
242 " is written to AFNI `fith' dataset \n"
243 "[-tcoef k prefixname] estimate of kth regression coefficient \n"
244 " along with t-test for the coefficient \n"
245 " is written to AFNI `fitt' dataset \n"
246 " \n"
247 " \n"
248 "The following commands generate one AFNI 'bucket' type dataset: \n"
249 " \n"
250 "[-bucket n prefixname] create one AFNI 'bucket' dataset having \n"
251 " n sub-bricks; n=0 creates default output;\n"
252 " output 'bucket' is written to prefixname \n"
253 "The mth sub-brick will contain: \n"
254 "[-brick m coef k label] kth parameter regression coefficient \n"
255 "[-brick m fstat label] F-stat for significance of regression \n"
256 "[-brick m rstat label] coefficient of multiple determination R^2 \n"
257 "[-brick m tstat k label] t-stat for kth regression coefficient \n"
258 "\n"
259 "[-datum DATUM] write the output in DATUM format. \n"
260 " Choose from short (default) or float.\n"
261 "\n" );
262
263 printf
264 (
265 "\n"
266 "N.B.: For this program, the user must specify 1 and only 1 sub-brick \n"
267 " with each -xydata command. That is, if an input dataset contains\n"
268 " more than 1 sub-brick, a sub-brick selector must be used, e.g.: \n"
269 " -xydata 2.17 4.59 7.18 'fred+orig[3]' \n"
270 );
271
272 printf("\n" MASTER_SHORTHELP_STRING ) ;
273
274 PRINT_COMPILE_DATE ; exit(0);
275 }
276
277
278 /*---------------------------------------------------------------------------*/
279
280 /** macro to open a dataset and make it ready for processing **/
281
282 #define DOPEN(ds,name) \
283 do{ int pv ; (ds) = THD_open_dataset((name)) ; \
284 if( !ISVALID_3DIM_DATASET((ds)) ){ \
285 fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \
286 if( (ds)->daxes->nxx!=nx || (ds)->daxes->nyy!=ny || \
287 (ds)->daxes->nzz!=nz ){ \
288 fprintf(stderr,"*** Axes mismatch: %s\n",(name)) ; exit(1) ; } \
289 if( DSET_NVALS((ds)) != 1 ){ \
290 fprintf(stderr,"*** Must specify 1 sub-brick: %s\n",(name));exit(1);}\
291 DSET_load((ds)) ; \
292 pv = DSET_PRINCIPAL_VALUE((ds)) ; \
293 if( DSET_ARRAY((ds),pv) == NULL ){ \
294 fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); } \
295 if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){ \
296 fprintf(stderr,"*** Can't use complex data: %s\n",(name)); exit(1); \
297 } \
298 break ; } while (0)
299
300
301 /*---------------------------------------------------------------------------*/
302
303 /** macro to return pointer to correct location in brick for current processing **/
304
305 #define SUB_POINTER(ds,vv,ind,ptr) \
306 do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){ \
307 default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1); \
308 case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ; \
309 (ptr) = (void *)( fim + (ind) ) ; \
310 } break ; \
311 case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ; \
312 (ptr) = (void *)( fim + (ind) ) ; \
313 } break ; \
314 case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ; \
315 (ptr) = (void *)( fim + (ind) ) ; \
316 } break ; } break ; } while(0)
317
318
319 /*---------------------------------------------------------------------------*/
320 /*
321 Get the dimensions of the 3d AFNI data sets.
322 */
323
get_dimensions(RA_options * option_data)324 void get_dimensions
325 (
326 RA_options * option_data /* user input options */
327 )
328
329 {
330
331 THD_3dim_dataset * dset=NULL;
332
333 /*----- read first dataset to get dimensions, etc. -----*/
334
335 dset = THD_open_dataset( option_data->first_dataset ) ;
336 CHECK_OPEN_ERROR(dset,option_data->first_dataset) ;
337
338 /*----- data set dimensions in voxels -----*/
339 option_data->nx = dset->daxes->nxx ;
340 option_data->ny = dset->daxes->nyy ;
341 option_data->nz = dset->daxes->nzz ;
342 option_data->nxyz = option_data->nx * option_data->ny * option_data->nz ;
343
344 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
345
346 }
347
348
349 /*---------------------------------------------------------------------------*/
350 /*
351 Read one AFNI data set from file 'filename'.
352 The data is converted to floating point (in ffim).
353 */
354
read_afni_data(RA_options * option_data,char * filename,int piece_len,int fim_offset,float * ffim)355 void read_afni_data
356 (
357 RA_options * option_data, /* user input options */
358 char * filename, /* input AFNI dataset file name */
359 int piece_len, /* number of voxels in current piece */
360 int fim_offset, /* array offset to current piece */
361 float * ffim /* input data in floating point format */
362 )
363
364 {
365 int iv; /* index number of intensity sub-brick */
366 THD_3dim_dataset * dset=NULL; /* data set pointer */
367 void * vfim = NULL; /* image data pointer */
368 int nx, ny, nz, nxyz; /* data set dimensions in voxels */
369
370 nx = option_data->nx;
371 ny = option_data->ny;
372 nz = option_data->nz;
373 nxyz = option_data->nxyz;
374
375
376 /*----- read in the data -----*/
377 DOPEN (dset, filename) ;
378 iv = DSET_PRINCIPAL_VALUE(dset) ;
379
380 /*----- convert it to floats (in ffim) -----*/
381 SUB_POINTER (dset, iv, fim_offset, vfim) ;
382 EDIT_coerce_scale_type (piece_len, DSET_BRICK_FACTOR(dset,iv),
383 DSET_BRICK_TYPE(dset,iv), vfim, /* input */
384 MRI_float ,ffim ) ; /* output */
385
386 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
387 }
388
389
390 /*---------------------------------------------------------------------------*/
391 /*
392 Routine to check whether one piece file already exists.
393 */
394
check_piece(char * filename)395 void check_piece
396 (
397 char * filename /* name of piece file */
398 )
399
400 {
401 FILE * far = NULL; /* pointer to temporary file */
402 char sfilename[MAX_NAME_LENGTH]; /* name of temporary file */
403 char message[MAX_NAME_LENGTH+128]; /* error message: fname+message */
404
405
406 /*----- see if piece file already exists -----*/
407 strcpy (sfilename, filename);
408 strcat (sfilename, SUFFIX);
409 far = fopen (sfilename, "r");
410 if (far != NULL)
411 {
412 sprintf (message, "temporary file %s already exists. ", sfilename);
413 RA_error (message);
414 }
415
416 }
417
418
419 /*---------------------------------------------------------------------------*/
420 /*
421 Routine to read a piece of floating point volume.
422 */
423
read_piece(char * filename,float * fin,int size)424 void read_piece
425 (
426 char * filename, /* root name of piece file */
427 float * fin, /* input floating point piece */
428 int size /* number of voxels in piece */
429 )
430
431 {
432 char sfilename[MAX_NAME_LENGTH]; /* piece file name */
433 char message[MAX_NAME_LENGTH+128]; /* error message: fname + message */
434 FILE * far = NULL; /* floating point input file */
435 int count; /* number of data items read from disk */
436
437
438 /*----- input file name -----*/
439 strcpy (sfilename, filename);
440 strcat (sfilename, SUFFIX);
441
442 /*----- open temporary data file for input -----*/
443 far = fopen (sfilename, "r");
444 if (far == NULL)
445 {
446 sprintf (message, "Unable to open temporary file %s", sfilename);
447 RA_error (message);
448 }
449
450 /*----- read 3d data file -----*/
451 count = fread (fin, sizeof(float), size, far);
452 fclose (far);
453
454 /*----- error in reading file? -----*/
455 if (count != size)
456 {
457 sprintf (message, "Error in reading temporary file %s", sfilename);
458 RA_error (message);
459 }
460 }
461
462
463 /*---------------------------------------------------------------------------*/
464 /*
465 Routine to write a piece of floating point volume to a (temporary)
466 disk file.
467 */
468
write_piece(char * filename,float * fout,int size)469 void write_piece
470 (
471 char * filename, /* root name of output piece file */
472 float * fout, /* floating point output piece */
473 int size /* number of voxels in piece */
474 )
475
476 {
477 char sfilename[MAX_NAME_LENGTH]; /* piece file name */
478 char message[MAX_NAME_LENGTH+128]; /* error message */
479 FILE * far = NULL; /* floating point output file */
480 int count; /* number of data items written to disk */
481
482
483 /*----- output file name -----*/
484 strcpy (sfilename, filename);
485 strcat (sfilename, SUFFIX);
486
487 /*----- first, see if file already exists -----*/
488 far = fopen (sfilename, "r");
489 if (far != NULL)
490 {
491 sprintf (message, "Temporary file %s already exists. ", sfilename);
492 RA_error (message);
493 }
494
495 /*----- open temporary data file for output -----*/
496 far = fopen (sfilename, "w");
497 if (far == NULL)
498 {
499 sprintf (message, "Unable to write temporary file %s ", sfilename);
500 RA_error (message);
501 }
502
503 /*----- write 3d data set -----*/
504 count = fwrite (fout, sizeof(float), size, far);
505 fclose (far);
506
507 /*----- error in writing file? -----*/
508 if (count != size)
509 {
510 sprintf (message, "Error in writing temporary file %s ", sfilename);
511 RA_error (message);
512 }
513
514 }
515
516
517 /*---------------------------------------------------------------------------*/
518 /*
519 Routine to delete a file containing a piece of floating point volume.
520 */
521
delete_piece(char * filename)522 void delete_piece
523 (
524 char * filename /* root name of piece file to be deleted */
525 )
526
527 {
528 char sfilename[MAX_NAME_LENGTH]; /* file name */
529
530 /*----- build file name -----*/
531 strcpy (sfilename, filename);
532 strcat (sfilename, SUFFIX);
533
534 /*----- delete 3d data file -----*/
535 remove (sfilename);
536
537 }
538
539
540 /*---------------------------------------------------------------------------*/
541 /*
542 Allocate memory space for multiple pieces.
543 */
544
allocate_pieces(matrix xdata,model * regmodel,RA_options * option_data,float *** yfimar,float ** freg_piece,float ** rsqr_piece,float *** coef_piece,float *** tcoef_piece)545 void allocate_pieces
546 (
547 matrix xdata, /* independent variable matrix */
548 model * regmodel, /* linear regression model */
549 RA_options * option_data, /* user input options */
550
551 float *** yfimar, /* array of pieces of Y-datasets */
552 float ** freg_piece, /* piece F-statistic for the full regression model */
553 float ** rsqr_piece, /* piece coeff. of multiple determination R^2 */
554 float *** coef_piece, /* piece regression parameters */
555 float *** tcoef_piece /* piece t-statistics for regression parameters */
556 )
557
558 {
559 int n; /* number of datasets */
560 int p; /* number of parameters in the full model */
561 int * flist = NULL; /* list of parameters in the full model */
562 int i; /* dataset index */
563 int ip; /* parameter index */
564 int ix; /* x-variable index */
565 int piece_size; /* number of voxels in datasset piece */
566 int ib; /* sub-brick index */
567 int nbricks; /* number of sub-bricks in bucket dataset */
568
569
570 /*----- initialize local variables -----*/
571 n = xdata.rows;
572 p = regmodel->p;
573 flist = regmodel->flist;
574 piece_size = option_data->piece_size;
575 nbricks = option_data->numbricks;
576
577
578 /*----- allocate memory space for input data -----*/
579 *yfimar = (float **) malloc (sizeof(float *) * n);
580 MTEST (*yfimar);
581 for (i = 0; i < n; i++)
582 {
583 (*yfimar)[i] = (float *) malloc(sizeof(float) * piece_size);
584 MTEST ((*yfimar)[i]);
585 }
586
587
588 /*----- allocate memory space for F-statistics data -----*/
589 for (ip = 0; ip < p; ip++)
590 {
591 ix = flist[ip];
592
593 if (option_data->fcoef_filename[ix] != NULL)
594 {
595 if (*freg_piece == NULL)
596 {
597 *freg_piece = (float *) malloc (sizeof(float) * piece_size);
598 MTEST (*freg_piece);
599 }
600 }
601 }
602
603 for (ib = 0; ib < nbricks; ib++)
604 {
605 if (option_data->brick_type[ib] == FUNC_FT_TYPE)
606 {
607 if (*freg_piece == NULL)
608 {
609 *freg_piece = (float *) malloc (sizeof(float) * piece_size);
610 MTEST (*freg_piece);
611 }
612 }
613 }
614
615
616 /*----- allocate memory space for coef. of mult. determination R^2 -----*/
617 for (ip = 0; ip < p; ip++)
618 {
619 ix = flist[ip];
620
621 if (option_data->rcoef_filename[ix] != NULL)
622 {
623 if (*rsqr_piece == NULL)
624 {
625 *rsqr_piece = (float *) malloc (sizeof(float) * piece_size);
626 MTEST (*rsqr_piece);
627 }
628 }
629 }
630
631 for (ib = 0; ib < nbricks; ib++)
632 {
633 if (option_data->brick_type[ib] == FUNC_THR_TYPE)
634 {
635 if (*rsqr_piece == NULL)
636 {
637 *rsqr_piece = (float *) malloc (sizeof(float) * piece_size);
638 MTEST (*rsqr_piece);
639 }
640 }
641 }
642
643
644 /*----- allocate memory space for parameters and t-stats -----*/
645 *coef_piece = (float **) malloc (sizeof(float *) * MAX_XVARS);
646 MTEST (*coef_piece);
647
648 *tcoef_piece = (float **) malloc (sizeof(float *) * MAX_XVARS);
649 MTEST (*tcoef_piece);
650
651 for (ix = 0; ix < MAX_XVARS; ix++)
652 {
653 (*coef_piece)[ix] = NULL;
654 (*tcoef_piece)[ix] = NULL;
655 }
656
657 for (ip = 0; ip < p; ip++)
658 {
659 ix = flist[ip];
660
661 if ( (option_data->fcoef_filename[ix] != NULL)
662 || (option_data->rcoef_filename[ix] != NULL)
663 || (option_data->tcoef_filename[ix] != NULL))
664 {
665 (*coef_piece)[ix] =
666 (float *) malloc (sizeof(float) * piece_size);
667 MTEST ((*coef_piece)[ix]);
668 }
669
670 if (option_data->tcoef_filename[ix] != NULL)
671 {
672 (*tcoef_piece)[ix] =
673 (float *) malloc (sizeof(float) * piece_size);
674 MTEST ((*tcoef_piece)[ix]);
675 }
676 }
677
678 for (ib = 0; ib < nbricks; ib++)
679 {
680 if (option_data->brick_type[ib] == FUNC_FIM_TYPE)
681 {
682 ix = option_data->brick_coef[ib];
683 if ((*coef_piece)[ix] == NULL)
684 {
685 (*coef_piece)[ix] = (float *) malloc (sizeof(float)*piece_size);
686 MTEST ((*coef_piece)[ix]);
687 }
688 }
689 }
690
691 for (ib = 0; ib < nbricks; ib++)
692 {
693 if (option_data->brick_type[ib] == FUNC_TT_TYPE)
694 {
695 ix = option_data->brick_coef[ib];
696 if ((*tcoef_piece)[ix] == NULL)
697 {
698 (*tcoef_piece)[ix] = (float *) malloc (sizeof(float)*piece_size);
699 MTEST ((*tcoef_piece)[ix]);
700 }
701 }
702 }
703
704 }
705
706
707 /*---------------------------------------------------------------------------*/
708 /*
709 Save multiple piece data into temporary files.
710 */
711
save_pieces(int piece,int piece_len,float * freg_piece,float * rsqr_piece,float ** coef_piece,float ** tcoef_piece)712 void save_pieces
713 (
714 int piece, /* current piece index */
715 int piece_len, /* number of voxels in current piece */
716
717 float * freg_piece, /* piece F-statistic for the full regression model */
718 float * rsqr_piece, /* piece coeff. of multiple determination R^2 */
719 float ** coef_piece, /* piece regression parameters */
720 float ** tcoef_piece /* piece t-statistics for regression parameters */
721 )
722
723 {
724 int ip; /* parameter index */
725 char filename[MAX_NAME_LENGTH]; /* name for temporary data file */
726
727
728 /*----- save piece containing F-statistics -----*/
729 if (freg_piece != NULL)
730 {
731 sprintf (filename, "freg.p%d", piece);
732 write_piece (filename, freg_piece, piece_len);
733 }
734
735
736 /*----- save piece containing R^2 -----*/
737 if (rsqr_piece != NULL)
738 {
739 sprintf (filename, "rsqr.p%d", piece);
740 write_piece (filename, rsqr_piece, piece_len);
741 }
742
743
744 /*----- save pieces containing regression coefficients -----*/
745 if (coef_piece != NULL)
746 {
747 for (ip = 0; ip < MAX_XVARS; ip++)
748 if (coef_piece[ip] != NULL)
749 {
750 sprintf (filename, "coef.%d.p%d", ip, piece);
751 write_piece (filename, coef_piece[ip], piece_len);
752 }
753 }
754
755
756 /*----- save pieces containing t-statistics -----*/
757 if (tcoef_piece != NULL)
758 {
759 for (ip = 0; ip < MAX_XVARS; ip++)
760 if (tcoef_piece[ip] != NULL)
761 {
762 sprintf (filename, "tcoef.%d.p%d", ip, piece);
763 write_piece (filename, tcoef_piece[ip], piece_len);
764 }
765 }
766
767 }
768
769
770 /*---------------------------------------------------------------------------*/
771 /*
772 Deallocate memory space for multiple pieces.
773 */
774
deallocate_pieces(int n,float *** yfimar,float ** freg_piece,float ** rsqr_piece,float *** coef_piece,float *** tcoef_piece)775 void deallocate_pieces
776 (
777 int n, /* number of datasets */
778 float *** yfimar, /* array of pieces of Y-datasets */
779 float ** freg_piece, /* piece F-statistic for the full regression model */
780 float ** rsqr_piece, /* piece coeff. of multiple determination R^2 */
781 float *** coef_piece, /* piece regression parameters */
782 float *** tcoef_piece /* piece t-statistics for regression parameters */
783 )
784
785 {
786 int i; /* dataset index */
787 int ip; /* parameter index */
788
789
790 /*----- deallocate memory space for input data -----*/
791 if (*yfimar != NULL)
792 {
793 for (i = 0; i < n; i++)
794 {
795 free ((*yfimar)[i]); (*yfimar)[i] = NULL;
796 }
797 free (*yfimar);
798 *yfimar = NULL;
799 }
800
801 /*----- deallocate memory space for F-statistics data -----*/
802 if (*freg_piece != NULL)
803 {
804 free (*freg_piece);
805 *freg_piece = NULL;
806 }
807
808 /*----- deallocate memory space for coef. of mult. determination R^2 -----*/
809 if (*rsqr_piece != NULL)
810 {
811 free (*rsqr_piece);
812 *rsqr_piece = NULL;
813 }
814
815 /*----- deallocate memory space for regression coefficients -----*/
816 if (*coef_piece != NULL)
817 {
818 for (ip = 0; ip < MAX_XVARS; ip++)
819 if ((*coef_piece)[ip] != NULL)
820 {
821 free ((*coef_piece)[ip]);
822 (*coef_piece)[ip] = NULL;
823 }
824 free (*coef_piece);
825 *coef_piece = NULL;
826 }
827
828 /*----- deallocate memory space for t-statistics -----*/
829 if (*tcoef_piece != NULL)
830 {
831 for (ip = 0; ip < MAX_XVARS; ip++)
832 if ((*tcoef_piece)[ip] != NULL)
833 {
834 free ((*tcoef_piece)[ip]);
835 (*tcoef_piece)[ip] = NULL;
836 }
837 free (*tcoef_piece);
838 *tcoef_piece = NULL;
839 }
840
841 }
842
843
844 /*---------------------------------------------------------------------------*/
845 /*
846 Check whether any of the piece files corresponding to a volume already exist.
847 */
848
check_volume(char * filename,int num_pieces)849 void check_volume
850 (
851 char * filename, /* root file name */
852 int num_pieces /* dataset is divided into this many pieces */
853 )
854
855 {
856 int piece; /* piece index */
857 char sfilename[MAX_NAME_LENGTH]; /* name for temporary data file */
858
859
860 /*----- loop over the temporary data file pieces -----*/
861 for (piece = 0; piece < num_pieces; piece++)
862 {
863 /*----- check for piece data file -----*/
864 sprintf (sfilename, "%s.p%d", filename, piece);
865 check_piece (sfilename);
866
867 } /* loop over pieces */
868
869 }
870
871
872 /*---------------------------------------------------------------------------*/
873 /*
874 Read 3d volume consisting of multiple piece files.
875 */
876
read_volume(char * filename,float * volume,int nxyz,int piece_size,int num_pieces)877 void read_volume
878 (
879 char * filename, /* root file name */
880 float * volume, /* floating point volume data */
881 int nxyz, /* number of voxels per volume */
882 int piece_size, /* number of voxels in dataset piece */
883 int num_pieces /* dataset is divided into this many pieces */
884 )
885
886 {
887 int piece; /* piece index */
888 int piece_len; /* number of voxels in current piece */
889 int fim_offset; /* array offset to current piece */
890 char sfilename[MAX_NAME_LENGTH]; /* name for temporary data file */
891
892
893 /*----- loop over the temporary data file pieces -----*/
894 for (piece = 0; piece < num_pieces; piece++)
895 {
896
897 /*----- current offset into volume -----*/
898 fim_offset = piece * piece_size;
899
900 /*----- size of current piece -----*/
901 if (piece < num_pieces-1)
902 piece_len = piece_size;
903 else
904 piece_len = nxyz - fim_offset;
905
906 /*----- read in the piece data -----*/
907 sprintf (sfilename, "%s.p%d", filename, piece);
908 read_piece (sfilename, volume + fim_offset, piece_len);
909
910 } /* loop over pieces */
911
912 }
913
914
915 /*---------------------------------------------------------------------------*/
916 /*
917 Delete 3d volume consisting of multiple piece files.
918 */
919
delete_volume(char * filename,int nxyz,int piece_size,int num_pieces)920 void delete_volume
921 (
922 char * filename, /* root file name */
923 int nxyz, /* number of voxels per volume */
924 int piece_size, /* number of voxels in dataset piece */
925 int num_pieces /* dataset is divided into this many pieces */
926 )
927
928 {
929 int piece; /* piece index */
930 char sfilename[MAX_NAME_LENGTH]; /* name for temporary data file */
931
932
933 /*----- loop over the temporary data file pieces -----*/
934 for (piece = 0; piece < num_pieces; piece++)
935 {
936
937 /*----- delete the piece data file -----*/
938 sprintf (sfilename, "%s.p%d", filename, piece);
939 delete_piece (sfilename);
940
941 } /* loop over pieces */
942
943 }
944
945
946 /*---------------------------------------------------------------------------*/
947 /*
948 Routine to initialize input options.
949 */
950
initialize_options(model * regmodel,RA_options * option_data)951 void initialize_options
952 (
953 model * regmodel, /* regression model */
954 RA_options * option_data /* user input options */
955 )
956
957 {
958 int ip; /* index */
959
960
961 regmodel->p = 0;
962 regmodel->q = 0;
963
964 option_data->datum = ILLEGAL_TYPE;
965 strcpy (option_data->session, "./");
966
967 option_data->yname = NULL;
968 option_data->first_dataset = NULL;
969
970 option_data->nx = 0;
971 option_data->ny = 0;
972 option_data->nz = 0;
973 option_data->nxyz = 0;
974
975 option_data->diskspace = 0;
976 option_data->workmem = 750; /* from 12 - 8 Dec 2005 [rickr] */
977 option_data->piece_size = 0;
978 option_data->num_pieces = 0;
979
980 option_data->rms_min = 0.0;
981 option_data->fdisp = -1.0;
982
983 option_data->levels = NULL;
984 option_data->counts = NULL;
985 option_data->c = 0;
986 option_data->flofmax = -1;
987
988 option_data->numf = 0;
989 option_data->numr = 0;
990 option_data->numc = 0;
991 option_data->numt = 0;
992
993 option_data->fcoef_filename = NULL;
994 option_data->rcoef_filename = NULL;
995 option_data->tcoef_filename = NULL;
996
997 option_data->numfiles = 0;
998
999 /*----- allocate memory for storing data file names -----*/
1000 option_data->yname
1001 = (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
1002 for (ip = 0; ip < MAX_OBSERVATIONS; ip++)
1003 option_data->yname[ip] = NULL;
1004
1005
1006 /*----- allocate memory space and initialize pointers for filenames -----*/
1007 option_data->fcoef_filename =
1008 (char **) malloc (sizeof(char *) * MAX_XVARS);
1009 if (option_data->fcoef_filename == NULL)
1010 RA_error ("Unable to allocate memory for fcoef_filename");
1011
1012 option_data->rcoef_filename =
1013 (char **) malloc (sizeof(char *) * MAX_XVARS);
1014 if (option_data->rcoef_filename == NULL)
1015 RA_error ("Unable to allocate memory for rcoef_filename");
1016
1017 option_data->tcoef_filename =
1018 (char **) malloc (sizeof(char *) * MAX_XVARS);
1019 if (option_data->tcoef_filename == NULL)
1020 RA_error ("Unable to allocate memory for tcoef_filename");
1021
1022 for (ip = 0; ip < MAX_XVARS; ip++)
1023 {
1024 option_data->fcoef_filename[ip] = NULL;
1025 option_data->rcoef_filename[ip] = NULL;
1026 option_data->tcoef_filename[ip] = NULL;
1027 }
1028
1029
1030 /*----- initialize bucket dataset options -----*/
1031 option_data->bucket_filename = NULL;
1032 option_data->numbricks = -1;
1033 option_data->brick_type = NULL;
1034 option_data->brick_coef = NULL;
1035 option_data->brick_label = NULL;
1036
1037 }
1038
1039
1040 /*---------------------------------------------------------------------------*/
1041 /*
1042 Routine to sort model variable index numbers. (This is for convenience only.)
1043 */
1044
sort_model_indices(model * regmodel)1045 void sort_model_indices
1046 (
1047 model * regmodel /* regression model */
1048 )
1049
1050 {
1051 int i, j, temp; /* model variable index numbers */
1052
1053
1054 /*----- sort full model indices into ascending order -----*/
1055 for (i = 0; i < regmodel->p - 1; i++)
1056 for (j = i+1; j < regmodel->p; j++)
1057 if (regmodel->flist[i] > regmodel->flist[j])
1058 {
1059 temp = regmodel->flist[i];
1060 regmodel->flist[i] = regmodel->flist[j];
1061 regmodel->flist[j] = temp;
1062 }
1063 else if (regmodel->flist[i] == regmodel->flist[j])
1064 RA_error ("Duplicate variable indices in model definition");
1065
1066
1067 /*----- sort reduced model indices into ascending order -----*/
1068 for (i = 0; i < regmodel->q - 1; i++)
1069 for (j = i+1; j < regmodel->q; j++)
1070 if (regmodel->rlist[i] > regmodel->rlist[j])
1071 {
1072 temp = regmodel->rlist[i];
1073 regmodel->rlist[i] = regmodel->rlist[j];
1074 regmodel->rlist[j] = temp;
1075 }
1076
1077 }
1078
1079
1080 /*---------------------------------------------------------------------------*/
1081 /*
1082 Routine to get user specified regression analysis inputs.
1083 */
1084
get_inputs(int argc,char ** argv,matrix * xdata,model * regmodel,RA_options * option_data)1085 void get_inputs
1086 (
1087 int argc, /* number of input arguments */
1088 char ** argv, /* array of input arguments */
1089 matrix * xdata, /* independent variable matrix */
1090 model * regmodel, /* linear regression model */
1091 RA_options * option_data /* user input options */
1092 )
1093
1094 #define MAX_BRICKS 100
1095 {
1096 int nopt = 1; /* input option argument counter */
1097 int ival, index; /* integer input */
1098 float fval; /* float input */
1099 int rows=0, cols=0; /* number rows and columns for x matrix */
1100 int irows=0, jcols; /* data point counters */
1101 THD_3dim_dataset * dset=NULL; /* test whether data set exists */
1102 char message[MAX_NAME_LENGTH]; /* error message */
1103 char label[MAX_NAME_LENGTH]; /* sub-brick label */
1104 int ibrick; /* sub-brick index */
1105 int nbricks; /* number of sub-bricks in the bucket */
1106 int p; /* number of parameters */
1107 int ip, ix; /* parameter indices */
1108
1109
1110 /*----- does user request help menu? -----*/
1111 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
1112
1113
1114 /*----- add to program log -----*/
1115 AFNI_logger (PROGRAM_NAME,argc,argv);
1116
1117
1118 /*----- initialize the input options -----*/
1119 initialize_options (regmodel, option_data);
1120
1121
1122 /*----- main loop over input options -----*/
1123 while (nopt < argc && argv[nopt][0] == '-')
1124 {
1125
1126 /*----- -datum type -----*/
1127 if( strncmp(argv[nopt],"-datum",6) == 0 ){
1128 if( ++nopt >= argc ) RA_error("need an argument after -datum!") ;
1129
1130 if( strcmp(argv[nopt],"short") == 0 ){
1131 option_data->datum = MRI_short ;
1132 } else if( strcmp(argv[nopt],"float") == 0 ){
1133 option_data->datum = MRI_float ;
1134 } else {
1135 char buf[256] ;
1136 sprintf(buf,"-datum of type '%s' is not supported in 3dRegAna.",
1137 argv[nopt] ) ;
1138 RA_error(buf) ;
1139 }
1140 nopt++ ; continue ; /* go to next arg */
1141 }
1142
1143
1144 /*----- -session dirname -----*/
1145 if( strncmp(argv[nopt],"-session",8) == 0 ){
1146 nopt++ ;
1147 if( nopt >= argc ) RA_error("need argument after -session!") ;
1148 strcpy(option_data->session , argv[nopt++]) ;
1149 continue ;
1150 }
1151
1152
1153 /*----- -diskspace -----*/
1154 if( strncmp(argv[nopt],"-diskspace",10) == 0 )
1155 {
1156 option_data->diskspace = 1;
1157 nopt++ ; continue ; /* go to next arg */
1158 }
1159
1160
1161 /*----- -workmem megabytes -----*/
1162
1163 if( strncmp(argv[nopt],"-workmem",6) == 0 ){
1164 nopt++ ;
1165 if( nopt >= argc ) RA_error ("need argument after -workmem!") ;
1166 sscanf (argv[nopt], "%d", &ival);
1167 if( ival <= 0 ) RA_error ("illegal argument after -workmem!") ;
1168 option_data->workmem = ival ;
1169 nopt++ ; continue ;
1170 }
1171
1172
1173 /*----- -rmsmin r -----*/
1174 if (strncmp(argv[nopt], "-rmsmin", 7) == 0)
1175 {
1176 nopt++;
1177 if (nopt >= argc) RA_error ("need argument after -rmsmin ");
1178 sscanf (argv[nopt], "%f", &fval);
1179 if (fval < 0.0)
1180 RA_error ("illegal argument after -rmsmin ");
1181 option_data->rms_min = fval;
1182 nopt++;
1183 continue;
1184 }
1185
1186
1187 /*----- -flof alpha -----*/
1188 if (strncmp(argv[nopt], "-flof", 6) == 0)
1189 {
1190 nopt++;
1191 if (nopt >= argc) RA_error ("need argument after -flof ");
1192 sscanf (argv[nopt], "%f", &fval);
1193 if ((fval < 0.0) || (fval > 1.0))
1194 RA_error ("illegal argument after -flof ");
1195 option_data->flofmax = fval;
1196 nopt++;
1197 continue;
1198 }
1199
1200
1201 /*----- -fdisp fval -----*/
1202 if (strncmp(argv[nopt], "-fdisp", 6) == 0)
1203 {
1204 nopt++;
1205 if (nopt >= argc) RA_error ("need argument after -fdisp ");
1206 sscanf (argv[nopt], "%f", &fval);
1207 option_data->fdisp = fval;
1208 nopt++;
1209 continue;
1210 }
1211
1212
1213 /*----- -rows n -----*/
1214 if (strncmp(argv[nopt], "-rows", 5) == 0)
1215 {
1216 nopt++;
1217 if (nopt >= argc) RA_error ("need argument after -rows ");
1218 sscanf (argv[nopt], "%d", &ival);
1219 if ((ival <= 0) || (ival > MAX_OBSERVATIONS))
1220 RA_error ("illegal argument after -rows ");
1221 rows = ival;
1222 nopt++;
1223 continue;
1224 }
1225
1226
1227 /*----- -cols m -----*/
1228 if (strncmp(argv[nopt], "-cols", 5) == 0)
1229 {
1230 nopt++;
1231 if (nopt >= argc) RA_error ("need argument after -cols ");
1232 sscanf (argv[nopt], "%d", &ival);
1233 if ((ival <= 0) || (ival > MAX_XVARS))
1234 RA_error ("illegal argument after -cols ");
1235 cols = ival;
1236 nopt++;
1237
1238 matrix_create (rows, cols+1, xdata);
1239 irows = -1;
1240 continue;
1241 }
1242
1243
1244 /*----- -xydata x1 ... xm y -----*/
1245 if (strncmp(argv[nopt], "-xydata", 7) == 0)
1246 {
1247 nopt++;
1248 if (nopt + cols >= argc)
1249 RA_error ("need cols+1 arguments after -xydata ");
1250
1251 irows++;
1252 if (irows > rows-1) RA_error ("too many data files");
1253
1254 xdata->elts[irows][0] = 1.0;
1255 for (jcols = 1; jcols <= cols; jcols++)
1256 {
1257 sscanf (argv[nopt], "%f", &fval);
1258 xdata->elts[irows][jcols] = fval;
1259 nopt++;
1260 }
1261
1262 /*--- check whether input files exist ---*/
1263 dset = THD_open_dataset( argv[nopt] ) ;
1264 CHECK_OPEN_ERROR(dset,argv[nopt]) ;
1265 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
1266
1267 option_data->yname[irows]
1268 = malloc (sizeof(char) * MAX_NAME_LENGTH);
1269 strcpy (option_data->yname[irows], argv[nopt]);
1270 nopt++;
1271 continue;
1272 }
1273
1274
1275 /*----- -model -----*/
1276 if (strncmp(argv[nopt], "-model", 6) == 0)
1277 {
1278 nopt++;
1279 if (nopt >= argc) RA_error ("need arguments after -model ");
1280
1281 while ((nopt < argc)
1282 && (strncmp(argv[nopt], "-", 1) != 0)
1283 && (strncmp(argv[nopt], ":", 1) != 0))
1284 {
1285 sscanf (argv[nopt], "%d", &ival);
1286 if ((ival <= 0) || (ival > cols))
1287 {
1288 fprintf(stderr,"bad Xi integer, '%s'\n\n", argv[nopt]);
1289 RA_error ("illegal argument after -model ");
1290 }
1291 regmodel->flist[regmodel->p] = ival;
1292
1293 regmodel->p++;
1294 nopt++;
1295 }
1296
1297 if (strncmp(argv[nopt], ":", 1) == 0)
1298 {
1299 if (isdigit(argv[nopt][1])) /* e.g. the user gives :3 ... */
1300 {
1301 fprintf(stderr,
1302 "please separate ':' from -model argument in '%s'\n\n",
1303 argv[nopt]);
1304 RA_error ("illegal argument after -model ");
1305 }
1306 nopt++;
1307
1308 while ((nopt < argc)
1309 && (strncmp(argv[nopt], "-", 1) != 0))
1310 {
1311 sscanf (argv[nopt], "%d", &ival);
1312 if ((ival < 0) || (ival > cols))
1313 {
1314 fprintf(stderr,"bad Xj integer, '%s'\n\n", argv[nopt]);
1315 RA_error ("illegal argument after -model ");
1316 }
1317 regmodel->flist[regmodel->p] = ival;
1318 regmodel->rlist[regmodel->q] = ival;
1319 regmodel->p++;
1320 regmodel->q++;
1321 nopt++;
1322 }
1323 }
1324 /* incorrect usage, complain and fail 24 Aug 2005 [rickr] */
1325 else
1326 {
1327 fprintf(stderr,"missing ':' between sets of -model parameters\n"
1328 " (is ':' attached to a number?)\n\n");
1329 RA_error ("bad -model usage ");
1330 }
1331
1332 /*----- sort model variable indices into ascending order -----*/
1333 sort_model_indices (regmodel);
1334
1335 continue;
1336 }
1337
1338
1339 /*----- -fcoef k filename -----*/
1340 if (strncmp(argv[nopt], "-fcoef", 6) == 0)
1341 {
1342 nopt++;
1343 if (nopt+1 >= argc) RA_error ("need 2 arguments after -fcoef ");
1344 sscanf (argv[nopt], "%d", &ival);
1345 if ((ival < 0) || (ival > cols))
1346 RA_error ("illegal argument after -fcoef ");
1347 index = ival;
1348 nopt++;
1349
1350 option_data->fcoef_filename[index] =
1351 malloc (sizeof(char) * MAX_NAME_LENGTH);
1352 if (option_data->fcoef_filename[index] == NULL)
1353 RA_error ("Unable to allocate memory for fcoef_filename");
1354 strcpy (option_data->fcoef_filename[index], argv[nopt]);
1355
1356 nopt++;
1357 continue;
1358 }
1359
1360
1361 /*----- -rcoef k filename -----*/
1362 if (strncmp(argv[nopt], "-rcoef", 6) == 0)
1363 {
1364 nopt++;
1365 if (nopt+1 >= argc) RA_error ("need 2 arguments after -rcoef ");
1366 sscanf (argv[nopt], "%d", &ival);
1367 if ((ival < 0) || (ival > cols))
1368 RA_error ("illegal argument after -rcoef ");
1369 index = ival;
1370 nopt++;
1371
1372 option_data->rcoef_filename[index] =
1373 malloc (sizeof(char) * MAX_NAME_LENGTH);
1374 if (option_data->rcoef_filename[index] == NULL)
1375 RA_error ("Unable to allocate memory for rcoef_filename");
1376 strcpy (option_data->rcoef_filename[index], argv[nopt]);
1377
1378 nopt++;
1379 continue;
1380 }
1381
1382
1383 /*----- -tcoef k filename -----*/
1384 if (strncmp(argv[nopt], "-tcoef", 6) == 0)
1385 {
1386 nopt++;
1387 if (nopt+1 >= argc) RA_error ("need 2 arguments after -tcoef ");
1388 sscanf (argv[nopt], "%d", &ival);
1389 if ((ival < 0) || (ival > cols))
1390 RA_error ("illegal argument after -tcoef ");
1391 index = ival;
1392 nopt++;
1393
1394 option_data->tcoef_filename[index] =
1395 malloc (sizeof(char) * MAX_NAME_LENGTH);
1396 if (option_data->tcoef_filename[index] == NULL)
1397 RA_error ("Unable to allocate memory for tcoef_filename");
1398 strcpy (option_data->tcoef_filename[index], argv[nopt]);
1399
1400 nopt++;
1401 continue;
1402 }
1403
1404
1405 /*----- -bucket n prefixname -----*/
1406 if (strncmp(argv[nopt], "-bucket", 7) == 0)
1407 {
1408 nopt++;
1409 if (nopt+1 >= argc) RA_error ("need 2 arguments after -bucket ");
1410 sscanf (argv[nopt], "%d", &ival);
1411 if ((ival < 0) || (ival > MAX_BRICKS))
1412 RA_error ("illegal argument after -bucket ");
1413 option_data->numbricks = ival;
1414 nopt++;
1415
1416 option_data->bucket_filename =
1417 malloc (sizeof(char) * MAX_NAME_LENGTH);
1418 if (option_data->bucket_filename == NULL)
1419 RA_error ("Unable to allocate memory for bucket_filename");
1420 strcpy (option_data->bucket_filename, argv[nopt]);
1421
1422 /*----- set number of sub-bricks in the bucket -----*/
1423 p = regmodel->p;
1424 if (ival == 0)
1425 nbricks = 2*p + 2;
1426 else
1427 nbricks = ival;
1428 option_data->numbricks = nbricks;
1429
1430 /*----- allocate memory and initialize bucket dataset options -----*/
1431 option_data->brick_type = malloc (sizeof(int) * nbricks);
1432 option_data->brick_coef = malloc (sizeof(int) * nbricks);
1433 option_data->brick_label = malloc (sizeof(char *) * nbricks);
1434 for (ibrick = 0; ibrick < nbricks; ibrick++)
1435 {
1436 option_data->brick_type[ibrick] = -1;
1437 option_data->brick_coef[ibrick] = -1;
1438 option_data->brick_label[ibrick] =
1439 malloc (sizeof(char) * MAX_NAME_LENGTH);
1440 }
1441
1442 if (ival == 0) /*----- throw everything into the bucket -----*/
1443 {
1444 for (ip = 0; ip < p; ip++)
1445 {
1446 ix = regmodel->flist[ip];
1447
1448 ibrick = 2*ip;
1449 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1450 option_data->brick_coef[ibrick] = ix;
1451 sprintf (label, "Coef #%d est.", ix);
1452 strcpy (option_data->brick_label[ibrick], label);
1453
1454 ibrick = 2*ip + 1;
1455 option_data->brick_type[ibrick] = FUNC_TT_TYPE;
1456 option_data->brick_coef[ibrick] = ix;
1457 sprintf (label, "Coef #%d t-stat", ix);
1458 strcpy (option_data->brick_label[ibrick], label);
1459 }
1460
1461 ibrick = 2*p;
1462 option_data->brick_type[ibrick] = FUNC_FT_TYPE;
1463 strcpy (option_data->brick_label[ibrick], "F-stat Regression");
1464
1465 ibrick = 2*p + 1;
1466 option_data->brick_type[ibrick] = FUNC_THR_TYPE;
1467 strcpy (option_data->brick_label[ibrick], "R^2 Regression");
1468
1469 }
1470
1471 nopt++;
1472 continue;
1473 }
1474
1475
1476 /*----- -brick m type k label -----*/
1477 if (strncmp(argv[nopt], "-brick", 6) == 0)
1478 {
1479 nopt++;
1480 if (nopt+2 >= argc) RA_error ("need more arguments after -brick ");
1481 sscanf (argv[nopt], "%d", &ibrick);
1482 if ((ibrick < 0) || (ibrick >= option_data->numbricks))
1483 RA_error ("illegal argument after -brick ");
1484 nopt++;
1485
1486 if (strncmp(argv[nopt], "coef", 4) == 0)
1487 {
1488 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
1489
1490 nopt++;
1491 sscanf (argv[nopt], "%d", &ival);
1492 if ((ival < 0) || (ival > cols))
1493 RA_error ("illegal argument after coef ");
1494 option_data->brick_coef[ibrick] = ival;
1495
1496 nopt++;
1497 if (nopt >= argc) RA_error ("need more arguments after -brick ");
1498 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1499
1500 }
1501 else if (strncmp(argv[nopt], "tstat", 4) == 0)
1502 {
1503 option_data->brick_type[ibrick] = FUNC_TT_TYPE;
1504
1505 nopt++;
1506 sscanf (argv[nopt], "%d", &ival);
1507 if ((ival < 0) || (ival > cols))
1508 RA_error ("illegal argument after tstat ");
1509 option_data->brick_coef[ibrick] = ival;
1510
1511 nopt++;
1512 if (nopt >= argc) RA_error ("need more arguments after -brick ");
1513 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1514
1515 }
1516 else if (strncmp(argv[nopt], "fstat", 4) == 0)
1517 {
1518 option_data->brick_type[ibrick] = FUNC_FT_TYPE;
1519
1520 nopt++;
1521 if (nopt >= argc) RA_error ("need more arguments after -brick ");
1522 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1523
1524 }
1525 else if (strncmp(argv[nopt], "rstat", 4) == 0)
1526 {
1527 option_data->brick_type[ibrick] = FUNC_THR_TYPE;
1528
1529 nopt++;
1530 if (nopt >= argc) RA_error ("need more arguments after -brick ");
1531 strcpy (option_data->brick_label[ibrick], argv[nopt]);
1532
1533 }
1534 else
1535 {
1536 sprintf(message,"Unrecognized option after -brick: %s\n",
1537 argv[nopt]);
1538 RA_error (message);
1539 }
1540
1541 nopt++;
1542 continue;
1543 }
1544
1545
1546 /*----- unknown command -----*/
1547 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
1548 RA_error (message);
1549
1550 }
1551
1552 /*----- check for fewer than expected datasets -----*/
1553 if (irows < rows-1)
1554 RA_error ("Fewer than expected datasets were entered");
1555 }
1556
1557
1558 /*---------------------------------------------------------------------------*/
1559 /*
1560 Count the number of output data volumes and output files required.
1561 */
1562
count_volumes_and_files(model * regmodel,RA_options * option_data)1563 void count_volumes_and_files
1564 (
1565 model * regmodel, /* regression model */
1566 RA_options * option_data /* user input options */
1567 )
1568
1569 {
1570 int ip; /* parameter index */
1571 int ix; /* x-variable index */
1572 int ib; /* sub-brick index */
1573 int p; /* number of parameters in the model */
1574 int nbricks; /* number of bricks in bucket dataset */
1575
1576
1577 /*----- initialize local variables -----*/
1578 nbricks = option_data->numbricks;
1579 p = regmodel->p;
1580
1581
1582 /*----- count number of volumes and files for F-statistics -----*/
1583 for (ip = 0; ip < p; ip++)
1584 {
1585 ix = regmodel->flist[ip];
1586
1587 if (option_data->fcoef_filename[ix] != NULL)
1588 {
1589 option_data->numf = 1;
1590 option_data->numfiles += 1;
1591 }
1592 }
1593
1594 for (ib = 0; ib < nbricks; ib++)
1595 if (option_data->brick_type[ib] == FUNC_FT_TYPE)
1596 option_data->numf = 1;
1597
1598
1599 /*----- count number of volumes and files for R^2 -----*/
1600 for (ip = 0; ip < p; ip++)
1601 {
1602 ix = regmodel->flist[ip];
1603
1604 if (option_data->rcoef_filename[ix] != NULL)
1605 {
1606 option_data->numr = 1;
1607 option_data->numfiles += 1;
1608 }
1609 }
1610
1611 for (ib = 0; ib < nbricks; ib++)
1612 if (option_data->brick_type[ib] == FUNC_THR_TYPE)
1613 option_data->numr = 1;
1614
1615
1616 /*----- count number of volumes and files for t-statistics -----*/
1617 for (ip = 0; ip < p; ip++)
1618 {
1619 ix = regmodel->flist[ip];
1620
1621 if (option_data->tcoef_filename[ix] != NULL)
1622 {
1623 option_data->numt += 1;
1624 option_data->numfiles += 1;
1625 }
1626
1627 else
1628 for (ib = 0; ib < nbricks; ib++)
1629 if ((option_data->brick_type[ib] == FUNC_TT_TYPE)
1630 && (option_data->brick_coef[ib] == ix))
1631 option_data->numt += 1;
1632 }
1633
1634
1635 /*----- count number of volumes for regression coefficients -----*/
1636 for (ip = 0; ip < p; ip++)
1637 {
1638 ix = regmodel->flist[ip];
1639
1640 if ( (option_data->fcoef_filename[ix] != NULL)
1641 || (option_data->rcoef_filename[ix] != NULL)
1642 || (option_data->tcoef_filename[ix] != NULL))
1643 option_data->numc += 1;
1644
1645 else
1646 for (ib = 0; ib < nbricks; ib++)
1647 if ((option_data->brick_type[ib] == FUNC_FIM_TYPE)
1648 && (option_data->brick_coef[ib] == ix))
1649 option_data->numc += 1;
1650 }
1651
1652 }
1653
1654
1655 /*---------------------------------------------------------------------------*/
1656 /*
1657 Break problem into smaller pieces.
1658 */
1659
break_into_pieces(int num_datasets,RA_options * option_data)1660 void break_into_pieces
1661 (
1662 int num_datasets, /* number of input datasets */
1663 RA_options * option_data /* user input options */
1664 )
1665
1666 {
1667 int num_vols; /* number of output volumes */
1668
1669
1670 /*----- count number of distinct floating point volumes required -----*/
1671 num_vols = option_data->numf + option_data->numr
1672 + option_data->numt + option_data->numc;
1673
1674 /*----- calculate number of voxels per piece -----*/
1675 option_data->piece_size = option_data->workmem * MEGA
1676 / ((num_datasets + num_vols) * sizeof(float));
1677 if (option_data->piece_size > option_data->nxyz)
1678 option_data->piece_size = option_data->nxyz;
1679
1680 /*----- calculate number of pieces per dataset -----*/
1681 option_data->num_pieces = (option_data->nxyz + option_data->piece_size - 1)
1682 / option_data->piece_size;
1683
1684 printf ("num_pieces = %d piece_size = %d \n",
1685 option_data->num_pieces, option_data->piece_size);
1686
1687 }
1688
1689
1690 /*---------------------------------------------------------------------------*/
1691 /*
1692 Identify repeat observations. (Repeat observations are necessary for
1693 performing the F-test for lack of fit.) Convert -flof input from alpha to
1694 corresponding F value.
1695 */
1696
identify_repeats(matrix * xdata,model * regmodel,RA_options * option_data)1697 void identify_repeats
1698 (
1699 matrix * xdata, /* independent variable matrix */
1700 model * regmodel, /* linear regression model */
1701 RA_options * option_data /* user input options */
1702 )
1703
1704 {
1705 int i, k; /* matrix row indices */
1706 int j; /* matrix column index */
1707 int match; /* boolean for repeat observation found */
1708
1709 int which; /* indicator for F-distribution calculation */
1710 double p, q; /* cumulative probabilities under F-dist. */
1711 double f; /* F value corresponding to probability p */
1712 double dfn, dfd; /* numerator and denominator dof */
1713 int status; /* calculation status */
1714 double bound; /* search bound */
1715
1716
1717 /*----- allocate memory -----*/
1718 option_data->levels = (int *) malloc (sizeof(int) * (xdata->rows));
1719 MTEST (option_data->levels);
1720 option_data->counts = (int *) malloc (sizeof(int) * (xdata->rows));
1721 MTEST (option_data->counts);
1722
1723
1724 /*----- initialization -----*/
1725 for (i = 0; i < xdata->rows; i++)
1726 option_data->counts[i] = 0;
1727 option_data->levels[0] = 0;
1728 option_data->counts[0] = 1;
1729 option_data->c = 1;
1730
1731
1732 /*----- loop over observations -----*/
1733 for (i = 1; i < xdata->rows; i++)
1734 {
1735 option_data->levels[i] = option_data->c;
1736
1737 /*----- determine if this is a repeat observation -----*/
1738 for (k = 0; k < i; k++)
1739 {
1740 match = 1;
1741 for (j = 1; j < xdata->cols; j++)
1742 if (xdata->elts[i][j] != xdata->elts[k][j]) match = 0;
1743
1744 if (match)
1745 {
1746 option_data->levels[i] = option_data->levels[k];
1747 break;
1748 }
1749 }
1750
1751 /*----- increment count of repeat observations -----*/
1752 k = option_data->levels[i];
1753 option_data->counts[k] ++;
1754
1755 /*----- increment count of distinct observation levels -----*/
1756 if (k == option_data->c) option_data->c++;
1757 }
1758
1759
1760 /*----- determine F value corresponding to alpha for lack of fit -----*/
1761 which = 2;
1762 q = option_data->flofmax;
1763 p = 1.0 - q;
1764 dfn = option_data->c - regmodel->p;
1765 dfd = xdata->rows - option_data->c;
1766 cdff (&which, &p, &q, &f, &dfn, &dfd, &status, &bound);
1767 if (status != 0) RA_error ("Error in calculating F - lack of fit ");
1768 option_data->flofmax = f;
1769
1770
1771 }
1772
1773
1774 /*---------------------------------------------------------------------------*/
1775 /*
1776 Routine to check for valid inputs.
1777 */
1778
check_for_valid_inputs(matrix * xdata,model * regmodel,RA_options * option_data)1779 void check_for_valid_inputs
1780 (
1781 matrix * xdata, /* independent variable matrix */
1782 model * regmodel, /* regression model */
1783 RA_options * option_data /* user input options */
1784 )
1785
1786 {
1787 int ib; /* sub-brick index */
1788
1789
1790 /*----- check data set dimensions -----*/
1791 if (xdata->cols < 2)
1792 RA_error ("Too few X variables ");
1793 if (xdata->rows < 3)
1794 RA_error ("Too few data sets for Y-observations ");
1795
1796 /*----- check model dimensions -----*/
1797 if (regmodel->q < 1)
1798 RA_error ("Reduced regression model is too small");
1799 if (regmodel->p <= regmodel->q)
1800 RA_error ("Full regression model is too small");
1801 if (xdata->rows <= regmodel->p)
1802 RA_error ("Too few data sets for fitting full regression model ");
1803 if (regmodel->rlist[0] != 0)
1804 RA_error ("Reduced model must include variable index 0");
1805
1806
1807 /*----- check for repeat observations -----*/
1808 if (option_data->flofmax >= 0.0)
1809 {
1810 if (xdata->rows <= option_data->c)
1811 RA_error ("Cannot conduct F-test for lack of fit (no repeat obs.) ");
1812 if (option_data->c <= regmodel->p)
1813 RA_error ("Cannot conduct F-test for lack of fit (c <= p) ");
1814 }
1815
1816 /*----- check bucket dataset parameters -----*/
1817 if (option_data->numbricks > 0)
1818 for (ib = 0; ib < option_data->numbricks; ib++)
1819 {
1820 if (option_data->brick_type[ib] < 0)
1821 RA_error ("Must specify contents of each brick in the bucket");
1822 }
1823 }
1824
1825
1826 /*---------------------------------------------------------------------------*/
1827 /*
1828 Check whether one output file already exists.
1829 */
1830
check_one_output_file(RA_options * option_data,char * filename)1831 void check_one_output_file
1832 (
1833 RA_options * option_data, /* user input options */
1834 char * filename /* output file name */
1835 )
1836
1837 {
1838 THD_3dim_dataset * dset=NULL; /* input afni data set pointer */
1839 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
1840 int ierror; /* number of errors in editing data */
1841
1842
1843 /*----- read first dataset -----*/
1844 dset = THD_open_dataset (option_data->first_dataset ) ;
1845 CHECK_OPEN_ERROR(dset,option_data->first_dataset) ;
1846
1847 /*-- make an empty copy of this dataset, for eventual output --*/
1848 new_dset = EDIT_empty_copy( dset ) ;
1849
1850
1851 ierror = EDIT_dset_items( new_dset ,
1852 ADN_prefix , filename ,
1853 ADN_label1 , filename ,
1854 ADN_directory_name , option_data->session ,
1855 ADN_self_name , filename ,
1856 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE :
1857 GEN_FUNC_TYPE ,
1858 ADN_none ) ;
1859
1860 if( ierror > 0 ){
1861 fprintf(stderr,
1862 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
1863 exit(1) ;
1864 }
1865
1866 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
1867 fprintf(stderr,
1868 "*** Output dataset file %s already exists--cannot continue!\a\n",
1869 new_dset->dblk->diskptr->header_name ) ;
1870 exit(1) ;
1871 }
1872
1873 /*----- deallocate memory -----*/
1874 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
1875 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
1876
1877 }
1878
1879
1880 /*---------------------------------------------------------------------------*/
1881 /*
1882 Routine to check whether any output files already exist.
1883 */
1884
check_output_files(RA_options * option_data)1885 void check_output_files
1886 (
1887 RA_options * option_data /* user input options */
1888 )
1889
1890 {
1891 int ix; /* x-variable index */
1892
1893
1894 for (ix = 0; ix < MAX_XVARS; ix++)
1895 {
1896 if (option_data->fcoef_filename[ix] != NULL)
1897 check_one_output_file (option_data, option_data->fcoef_filename[ix]);
1898 if (option_data->rcoef_filename[ix] != NULL)
1899 check_one_output_file (option_data, option_data->rcoef_filename[ix]);
1900 if (option_data->tcoef_filename[ix] != NULL)
1901 check_one_output_file (option_data, option_data->tcoef_filename[ix]);
1902 }
1903
1904 if (option_data->bucket_filename != NULL)
1905 check_one_output_file (option_data, option_data->bucket_filename);
1906
1907 }
1908
1909
1910 /*---------------------------------------------------------------------------*/
1911 /*
1912 Check whether any temporary data files already exist.
1913 */
1914
check_temporary_files(matrix * xdata,model * regmodel,RA_options * option_data)1915 void check_temporary_files
1916 (
1917 matrix * xdata, /* independent variable matrix */
1918 model * regmodel, /* linear regression model */
1919 RA_options * option_data /* user input options */
1920 )
1921
1922 {
1923 int p; /* number of parameters in full model */
1924 int ip, ix; /* parameter index */
1925 int num_pieces; /* dataset is divided into this many pieces */
1926 char filename[MAX_NAME_LENGTH]; /* name for temporary data file */
1927
1928
1929 /*----- initialize local variables -----*/
1930 p = regmodel->p;
1931 num_pieces = option_data->num_pieces;
1932
1933
1934 /*----- check for F-statistics data file -----*/
1935 if (option_data->numf > 0)
1936 {
1937 strcpy (filename, "freg");
1938 check_volume (filename, num_pieces);
1939 }
1940
1941
1942 /*----- check for R^2 data file -----*/
1943 if (option_data->numr > 0)
1944 {
1945 strcpy (filename, "rsqr");
1946 check_volume (filename, num_pieces);
1947 }
1948
1949
1950 /*----- check for t-statistics data files -----*/
1951 if (option_data->numt > 0)
1952 {
1953 for (ip = 0; ip < p; ip++)
1954 {
1955 ix = regmodel->flist[ip];
1956
1957 if (option_data->tcoef_filename[ix] != NULL)
1958 {
1959 sprintf (filename, "tcoef.%d", ix);
1960 check_volume (filename, num_pieces);
1961 }
1962 }
1963 }
1964
1965
1966 /*----- check for regression coefficients data files -----*/
1967 if (option_data->numc > 0)
1968 {
1969 for (ip = 0; ip < p; ip++)
1970 {
1971 ix = regmodel->flist[ip];
1972
1973 if ( (option_data->fcoef_filename[ix] != NULL)
1974 || (option_data->rcoef_filename[ix] != NULL)
1975 || (option_data->tcoef_filename[ix] != NULL) )
1976 {
1977 sprintf (filename, "coef.%d", ix);
1978 check_volume (filename, num_pieces);
1979 }
1980 }
1981 }
1982
1983
1984 }
1985
1986
1987 /*---------------------------------------------------------------------------*/
1988 /*
1989 Routine to determine if sufficient disk space exists for storing
1990 the temporary data files.
1991 */
1992
check_disk_space(RA_options * option_data)1993 void check_disk_space
1994 (
1995 RA_options * option_data /* user input options */
1996 )
1997
1998 {
1999 char ch; /* user response */
2000 int nxyz; /* number of voxels per image */
2001 int nmax; /* maximum number of disk files */
2002 char filename[MAX_NAME_LENGTH]; /* output file name */
2003
2004
2005 /*----- initialize local variables -----*/
2006 nxyz = option_data->nxyz;
2007
2008 /*----- first, determine space required for temporary volume data -----*/
2009 nmax = 4 * nxyz * (option_data->numf + option_data->numr + option_data->numt
2010 + option_data->numc);
2011
2012 /*----- determine additional space required for output files -----*/
2013 nmax += 4 * nxyz * option_data->numfiles;
2014
2015 printf("\nThis problem requires approximately %d MB of free disk space.\n",
2016 (nmax/1000000) + 1);
2017
2018
2019 /*----- Determine how much disk space is available. -----*/
2020 printf ("Summary of available disk space: \n\n");
2021 system (DF);
2022 printf ("\nDo you wish to proceed? ");
2023 ch = getchar();
2024 if ((ch != 'Y') && (ch != 'y')) exit(0);
2025 printf ("\n");
2026 }
2027
2028
2029 /*---------------------------------------------------------------------------*/
2030 /*
2031 Program initialization.
2032 */
2033
initialize_program(int argc,char ** argv,matrix * xdata,model * regmodel,RA_options * option_data)2034 void initialize_program
2035 (
2036 int argc, /* number of input arguments */
2037 char ** argv, /* array of input arguments */
2038 matrix * xdata, /* independent variable matrix */
2039 model * regmodel, /* regression model */
2040 RA_options * option_data /* user input options */
2041 )
2042
2043 {
2044
2045 /*----- save command line for history notes -----*/
2046 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
2047
2048
2049 /*----- create independent variable data matrix -----*/
2050 matrix_initialize (xdata);
2051
2052
2053 /*----- get all operator inputs -----*/
2054 get_inputs (argc, argv, xdata, regmodel, option_data);
2055
2056
2057 /*----- use first data set to get data set dimensions -----*/
2058 option_data->first_dataset = option_data->yname[0];
2059 get_dimensions (option_data);
2060 printf ("Data set dimensions: nx = %d ny = %d nz = %d nxyz = %d \n",
2061 option_data->nx, option_data->ny, option_data->nz, option_data->nxyz);
2062
2063
2064 /*----- count the number of output volumes and files required -----*/
2065 count_volumes_and_files (regmodel, option_data);
2066
2067
2068 /*----- break problem into smaller pieces -----*/
2069 break_into_pieces (xdata->rows, option_data);
2070
2071
2072 /*----- identify repeat observations -----*/
2073 if (option_data->flofmax >= 0.0)
2074 identify_repeats (xdata, regmodel, option_data);
2075
2076
2077 /*----- check for valid inputs -----*/
2078 check_for_valid_inputs (xdata, regmodel, option_data);
2079
2080
2081 /*----- check whether any of the output files already exist -----*/
2082 if( THD_deathcon() ) check_output_files (option_data);
2083
2084
2085 /*----- check whether temporary files already exist -----*/
2086 check_temporary_files (xdata, regmodel, option_data);
2087
2088
2089 /*----- check whether there is sufficient disk space -----*/
2090 if (option_data->diskspace) check_disk_space (option_data);
2091
2092 }
2093
2094
2095 /*---------------------------------------------------------------------------*/
2096 /*
2097 Perform initialization required for the regression analysis.
2098 */
2099
init_regression_analysis(int p,int q,int * flist,int * rlist,matrix xdata,matrix * x_full,matrix * xtxinv_full,matrix * xtxinvxt_full,matrix * x_base,matrix * xtxinvxt_base,matrix * x_rdcd,matrix * xtxinvxt_rdcd)2100 void init_regression_analysis
2101 (
2102 int p, /* number of parameters in the full model */
2103 int q, /* number of parameters in the rdcd model */
2104 int * flist, /* list of parameters in the full model */
2105 int * rlist, /* list of parameters in the rdcd model */
2106 matrix xdata, /* independent variable matrix */
2107 matrix * x_full, /* extracted X matrix for full model */
2108 matrix * xtxinv_full, /* matrix: 1/(X'X) for full model */
2109 matrix * xtxinvxt_full, /* matrix: (1/(X'X))X' for full model */
2110 matrix * x_base, /* extracted X matrix for baseline model */
2111 matrix * xtxinvxt_base, /* matrix: (1/(X'X))X' for baseline model */
2112 matrix * x_rdcd, /* extracted X matrices for reduced models */
2113 matrix * xtxinvxt_rdcd /* matrix: (1/(X'X))X' for reduced models */
2114 )
2115
2116 {
2117 int zlist[MAX_XVARS]; /* list of parameters in constant model */
2118 int ip; /* parameter index */
2119 matrix xtxinv_temp; /* intermediate results */
2120
2121
2122 /*----- Initialize matrix -----*/
2123 matrix_initialize (&xtxinv_temp);
2124
2125
2126 /*----- Initialize list of parameters in the constant model -----*/
2127 for (ip = 0; ip < MAX_XVARS; ip++)
2128 zlist[ip] = 0;
2129
2130
2131 /*----- Calculate constant matrices which will be needed later -----*/
2132 calc_matrices (xdata, 1, zlist, x_base, &xtxinv_temp, xtxinvxt_base);
2133 calc_matrices (xdata, q, rlist, x_rdcd, &xtxinv_temp, xtxinvxt_rdcd);
2134 calc_matrices (xdata, p, flist, x_full, xtxinv_full, xtxinvxt_full);
2135
2136
2137 /*----- Destroy matrix -----*/
2138 matrix_destroy (&xtxinv_temp);
2139
2140 }
2141
2142
2143 /*---------------------------------------------------------------------------*/
2144 /*
2145 Perform the entire regression analysis for one voxel.
2146 */
2147
regression_analysis(int N,int p,int q,matrix x_full,matrix xtxinv_full,matrix xtxinvxt_full,matrix x_base,matrix xtxinvxt_base,matrix x_rdcd,matrix xtxinvxt_rdcd,vector y,float rms_min,int * levels,int * counts,int c,float flofmax,float * flof,vector * coef_full,vector * scoef_full,vector * tcoef_full,float * freg,float * rsqr)2148 void regression_analysis
2149 (
2150 int N, /* number of data points */
2151 int p, /* number of parameters in the full model */
2152 int q, /* number of parameters in the rdcd model */
2153 matrix x_full, /* extracted X matrix for full model */
2154 matrix xtxinv_full, /* matrix: 1/(X'X) for full model */
2155 matrix xtxinvxt_full, /* matrix: (1/(X'X))X' for full model */
2156 matrix x_base, /* extracted X matrix for baseline model */
2157 matrix xtxinvxt_base, /* matrix: (1/(X'X))X' for baseline model */
2158 matrix x_rdcd, /* extracted X matrix for reduced model */
2159 matrix xtxinvxt_rdcd, /* matrix: (1/(X'X))X' for reduced model */
2160 vector y, /* vector of measured data */
2161 float rms_min, /* minimum rms error to reject zero model */
2162 int * levels, /* indices for repeat observations */
2163 int * counts, /* number of observations at each level */
2164 int c, /* number of unique rows in ind. var. matrix */
2165 float flofmax, /* max. allowed F-stat due to lack of fit */
2166 float * flof, /* F-statistic for lack of fit */
2167 vector * coef_full, /* regression parameters */
2168 vector * scoef_full, /* std. devs. for regression parameters */
2169 vector * tcoef_full, /* t-statistics for regression parameters */
2170 float * freg, /* F-statistic for the full regression model */
2171 float * rsqr /* coeff. of multiple determination R^2 */
2172 )
2173
2174 {
2175 float sse_base; /* error sum of squares, baseline model */
2176 float sse_rdcd; /* error sum of squares, reduced model */
2177 float sse_full; /* error sum of squares, full model */
2178 float sspe; /* pure error sum of squares */
2179 vector coef_temp; /* intermediate results */
2180
2181
2182 /*----- Initialization -----*/
2183 vector_initialize (&coef_temp);
2184
2185
2186 /*----- Calculate regression coefficients for baseline model -----*/
2187 calc_coef (xtxinvxt_base, y, &coef_temp);
2188
2189
2190 /*----- Calculate the error sum of squares for the baseline model -----*/
2191 sse_base = calc_sse (x_base, coef_temp, y);
2192
2193
2194 /*----- Stop here if variation about baseline is sufficiently low -----*/
2195 if (sqrt(sse_base/N) < rms_min)
2196 {
2197 vector_create (p, coef_full);
2198 vector_create (p, scoef_full);
2199 vector_create (p, tcoef_full);
2200 *freg = 0.0;
2201 *rsqr = 0.0;
2202 vector_destroy (&coef_temp);
2203 return;
2204 }
2205
2206
2207 /*----- Calculate regression coefficients for the full model -----*/
2208 calc_coef (xtxinvxt_full, y, coef_full);
2209
2210
2211 /*----- Calculate the error sum of squares for the full model -----*/
2212 sse_full = calc_sse (x_full, *coef_full, y);
2213
2214
2215 /*----- Calculate t-statistics for the regression coefficients -----*/
2216 calc_tcoef (N, p, sse_full, xtxinv_full,
2217 *coef_full, scoef_full, tcoef_full);
2218
2219
2220 /*----- Test for lack of fit -----*/
2221 if (flofmax > 0.0)
2222 {
2223 /*----- Calculate the pure error sum of squares -----*/
2224 sspe = calc_sspe (y, levels, counts, c);
2225
2226 /*----- Calculate F-statistic for lack of fit -----*/
2227 *flof = calc_flof (N, p, c, sse_full, sspe);
2228
2229 if (*flof > flofmax)
2230 {
2231 vector_create (p, coef_full);
2232 vector_create (p, scoef_full);
2233 vector_create (p, tcoef_full);
2234 *freg = 0.0;
2235 *rsqr = 0.0;
2236 vector_destroy (&coef_temp);
2237 return;
2238 }
2239 }
2240 else
2241 *flof = -1.0;
2242
2243
2244 /*----- Calculate regression coefficients for reduced model -----*/
2245 calc_coef (xtxinvxt_rdcd, y, &coef_temp);
2246
2247
2248 /*----- Calculate the error sum of squares for the reduced model -----*/
2249 sse_rdcd = calc_sse (x_rdcd, coef_temp, y);
2250
2251
2252 /*----- Calculate F-statistic for significance of the regression -----*/
2253 *freg = calc_freg (N, p, q, sse_full, sse_rdcd);
2254
2255
2256 /*----- Calculate coefficient of multiple determination R^2 -----*/
2257 *rsqr = calc_rsqr (sse_full, sse_base);
2258
2259
2260 /*----- Dispose of vector -----*/
2261 vector_destroy (&coef_temp);
2262
2263 }
2264
2265
2266 /*---------------------------------------------------------------------------*/
2267 /*
2268 Save results for current voxel into piece data for output later.
2269 */
2270
save_voxel(int iv,vector y,float fdisp,model * regmodel,float flof,vector coef,vector tcoef,float freg,float rsqr,float * freg_piece,float * rsqr_piece,float ** coef_piece,float ** tcoef_piece)2271 void save_voxel
2272 (
2273 int iv, /* current voxel number within piece */
2274 vector y, /* vector of measured data */
2275 float fdisp, /* minimum F-statistic for display */
2276 model * regmodel, /* linear regression model */
2277 float flof, /* F-statistic for lack of fit */
2278 vector coef, /* regression parameters */
2279 vector tcoef, /* t-statistics for regression parameters */
2280 float freg, /* F-statistic for the full regression model */
2281 float rsqr, /* coeff. of multiple determination R^2 */
2282
2283 float * freg_piece, /* piece F-statistic for the full regression model */
2284 float * rsqr_piece, /* piece coeff. of multiple determination R^2 */
2285 float ** coef_piece, /* piece regression parameters */
2286 float ** tcoef_piece /* piece t-statistics for regression parameters */
2287 )
2288
2289 {
2290 int ip; /* parameter index */
2291 int ix; /* x-variable index */
2292
2293
2294 /*----- save regression results into piece data -----*/
2295 if (freg_piece != NULL) freg_piece[iv] = freg;
2296 if (rsqr_piece != NULL) rsqr_piece[iv] = rsqr;
2297
2298
2299 /*----- save regression parameter estimates -----*/
2300 for (ip = 0; ip < regmodel->p; ip++)
2301 {
2302 ix = regmodel->flist[ip];
2303
2304 if (coef_piece[ix] != NULL) coef_piece[ix][iv] = coef.elts[ip];
2305
2306 if (tcoef_piece[ix] != NULL) tcoef_piece[ix][iv] = tcoef.elts[ip];
2307
2308 }
2309
2310
2311 /*----- if so requested, display results for this voxel -----*/
2312 if ((fdisp >= 0.0) && (freg >= fdisp))
2313 {
2314 printf ("\n\nVoxel #%d: \n", iv);
2315 printf ("\nY data: \n");
2316 for (ip = 0; ip < y.dim; ip++)
2317 printf ("Y[%d] = %f \n", ip, y.elts[ip]);
2318
2319 if (flof >= 0.0) printf ("\nF lack of fit = %f \n", flof);
2320 printf ("\nF regression = %f \n", freg);
2321 printf ("R-squared = %f \n", rsqr);
2322
2323 printf ("\nRegression Coefficients: \n");
2324 for (ip = 0; ip < coef.dim; ip++)
2325 {
2326 ix = regmodel->flist[ip];
2327 printf ("b[%d] = %f ", ix, coef.elts[ip]);
2328 printf ("t[%d] = %f \n", ix, tcoef.elts[ip]);
2329 }
2330
2331 }
2332
2333 }
2334
2335 /*---------------------------------------------------------------------------*/
2336 /*
2337 Calculate the multiple linear regression analysis for all voxels
2338 (by breaking the datasets into smaller pieces, if necessary).
2339 */
2340
calculate_results(matrix xdata,model * regmodel,RA_options * option_data)2341 void calculate_results
2342 (
2343 matrix xdata, /* independent variable matrix */
2344 model * regmodel, /* linear regression model */
2345 RA_options * option_data /* user input options */
2346 )
2347
2348 {
2349 int * flist = NULL; /* list of parameters in the full model */
2350 int * rlist = NULL; /* list of parameters in the rdcd model */
2351 int n; /* number of data points */
2352 int p; /* number of parameters in the full model */
2353 int q; /* number of parameters in the rdcd model */
2354 float flof; /* F-statistic for lack of fit */
2355 float freg; /* F-statistic for the full regression model */
2356 float rsqr; /* coeff. of multiple determination R^2 */
2357 vector coef; /* regression parameters */
2358 vector scoef; /* std. devs. for regression parameters */
2359 vector tcoef; /* t-statistics for regression parameters */
2360
2361 matrix x_full; /* extracted X matrix for full model */
2362 matrix xtxinv_full; /* matrix: 1/(X'X) for full model */
2363 matrix xtxinvxt_full; /* matrix: (1/(X'X))X' for full model */
2364 matrix x_base; /* extracted X matrix for baseline model */
2365 matrix xtxinvxt_base; /* matrix: (1/(X'X))X' for baseline model */
2366 matrix x_rdcd; /* extracted X matrix for reduced model */
2367 matrix xtxinvxt_rdcd; /* matrix: (1/(X'X))X' for reduced model */
2368 vector y; /* vector of measured data */
2369
2370 int i; /* dataset index */
2371 int nxyz; /* number of voxels per dataset */
2372 int num_datasets; /* total number of datasets */
2373 int piece_size; /* number of voxels in dataset piece */
2374 int num_pieces; /* dataset is divided into this many pieces */
2375 int piece; /* piece index */
2376 int piece_len; /* number of voxels in current piece */
2377 int fim_offset; /* array offset to current piece */
2378 int ivox; /* index of voxel in current piece */
2379 int nvox; /* index of voxel within entire volume */
2380
2381 float ** yfimar = NULL; /* array of pieces of Y-datasets */
2382 float * freg_piece = NULL; /* piece for F-statistics */
2383 float * rsqr_piece = NULL; /* piece for R^2 */
2384 float ** coef_piece = NULL; /* pieces for regression coefficients */
2385 float ** tcoef_piece = NULL; /* pieces for t-statistics */
2386
2387
2388 /*----- Initialize matrices and vectors -----*/
2389 matrix_initialize (&x_full);
2390 matrix_initialize (&xtxinv_full);
2391 matrix_initialize (&xtxinvxt_full);
2392 matrix_initialize (&x_base);
2393 matrix_initialize (&xtxinvxt_base);
2394 matrix_initialize (&x_rdcd);
2395 matrix_initialize (&xtxinvxt_rdcd);
2396 vector_initialize (&coef);
2397 vector_initialize (&scoef);
2398 vector_initialize (&tcoef);
2399 vector_initialize (&y);
2400
2401
2402 /*----- initialize local variables -----*/
2403 n = xdata.rows;
2404 p = regmodel->p;
2405 flist = regmodel->flist;
2406 q = regmodel->q;
2407 rlist = regmodel->rlist;
2408 nxyz = option_data->nxyz;
2409 piece_size = option_data->piece_size;
2410 num_pieces = option_data->num_pieces;
2411
2412
2413 /*----- allocate memory space for pieces -----*/
2414 allocate_pieces (xdata, regmodel, option_data, &yfimar,
2415 &freg_piece, &rsqr_piece, &coef_piece, &tcoef_piece);
2416
2417
2418 /*----- initialization for the regression analysis -----*/
2419 init_regression_analysis (p, q, flist, rlist, xdata,
2420 &x_full, &xtxinv_full, &xtxinvxt_full,
2421 &x_base, &xtxinvxt_base, &x_rdcd, &xtxinvxt_rdcd);
2422
2423
2424 vector_create (n, &y);
2425
2426
2427 if (option_data->fdisp >= 0)
2428 {
2429 printf ("\n");
2430 printf ("X matrix: \n");
2431 matrix_print (xdata);
2432 }
2433
2434
2435 /*----- loop over the pieces of the input datasets -----*/
2436 nvox = -1;
2437 for (piece = 0; piece < num_pieces; piece++)
2438 {
2439 printf ("piece = %d \n", piece);
2440 fim_offset = piece * piece_size;
2441 if (piece < num_pieces-1)
2442 piece_len = piece_size;
2443 else
2444 piece_len = nxyz - fim_offset;
2445
2446 /*----- read in the Y-data -----*/
2447 for (i = 0; i < n; i++)
2448 read_afni_data (option_data, option_data->yname[i],
2449 piece_len, fim_offset, yfimar[i]);
2450
2451
2452 /*----- loop over voxels in this piece -----*/
2453 for (ivox = 0; ivox < piece_len; ivox++)
2454 {
2455 nvox += 1;
2456
2457
2458 /*----- extract Y-data for this voxel -----*/
2459 for (i = 0; i < n; i++)
2460 y.elts[i] = yfimar[i][ivox];
2461
2462
2463 /*----- calculate results for this voxel -----*/
2464 regression_analysis (n, p, q,
2465 x_full, xtxinv_full, xtxinvxt_full, x_base,
2466 xtxinvxt_base, x_rdcd, xtxinvxt_rdcd,
2467 y, option_data->rms_min, option_data->levels,
2468 option_data->counts, option_data->c,
2469 option_data->flofmax, &flof,
2470 &coef, &scoef, &tcoef, &freg, &rsqr);
2471
2472
2473 /*----- save results for this voxel -----*/
2474 save_voxel (ivox, y, option_data->fdisp,
2475 regmodel, flof, coef, tcoef, freg, rsqr,
2476 freg_piece, rsqr_piece, coef_piece, tcoef_piece);
2477
2478
2479 } /* loop over voxels within this piece */
2480
2481
2482 /*----- save piece data into external files -----*/
2483 save_pieces (piece, piece_len,
2484 freg_piece, rsqr_piece, coef_piece, tcoef_piece);
2485
2486
2487 } /* loop over pieces */
2488
2489
2490 /*----- deallocate memory for pieces -----*/
2491 deallocate_pieces (n, &yfimar, &freg_piece, &rsqr_piece,
2492 &coef_piece, &tcoef_piece);
2493
2494
2495 /*----- Dispose of matrices -----*/
2496 vector_destroy (&y);
2497 vector_destroy (&tcoef);
2498 vector_destroy (&scoef);
2499 vector_destroy (&coef);
2500 matrix_destroy (&xtxinvxt_rdcd);
2501 matrix_destroy (&x_rdcd);
2502 matrix_destroy (&xtxinvxt_base);
2503 matrix_destroy (&x_base);
2504 matrix_destroy (&xtxinvxt_full);
2505 matrix_destroy (&xtxinv_full);
2506 matrix_destroy (&x_full);
2507
2508 }
2509
2510 /*---------------------------------------------------------------------------*/
2511 /*
2512 Routine to write one AFNI data set.
2513
2514 This data set may be either `fith' type (intensity + threshold)
2515 or `fitt' type (intensity + t-statistic)
2516 or `fift' type (intensity + F-statistic).
2517
2518 The intensity data is in ffim, and the corresponding statistic is in ftr.
2519
2520 */
2521
write_afni_data(RA_options * option_data,char * filename,float * ffim,float * ftr,int func_type,int numdof,int dendof)2522 void write_afni_data
2523 (
2524 RA_options * option_data, /* user input options */
2525 char * filename, /* output file name */
2526 float * ffim, /* volume of intensity data */
2527 float * ftr, /* volume of test statistics */
2528 int func_type, /* afni data set type */
2529 int numdof, /* numerator degrees of freedom */
2530 int dendof /* denominator degrees of freedom */
2531 )
2532
2533 {
2534 int nxyz; /* number of voxels */
2535 int ii; /* voxel index */
2536 THD_3dim_dataset * dset=NULL; /* input afni data set pointer */
2537 THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
2538 int ierror; /* number of errors in editing data */
2539 int ibuf[32]; /* integer buffer */
2540 float fbuf[MAX_STAT_AUX]; /* float buffer */
2541 float fimfac; /* scale factor for short data */
2542 int output_datum; /* data type for output data */
2543 short * tsp = NULL; /* 2nd sub-brick data pointer */
2544 void * vdif = NULL; /* 1st sub-brick data pointer */
2545 float top, bot, func_scale_short=0.0; /* parameters for scaling data */
2546 int top_ss, bot_ss=0; /* 2nd sub-brick value limits */
2547 char label[80]; /* label for output file history */
2548
2549
2550 /*----- initialize local variables -----*/
2551 nxyz = option_data->nxyz;
2552
2553 /*----- read first dataset -----*/
2554 dset = THD_open_dataset (option_data->first_dataset) ;
2555 CHECK_OPEN_ERROR(dset,option_data->first_dataset) ;
2556
2557
2558 /*-- make an empty copy of this dataset, for eventual output --*/
2559 new_dset = EDIT_empty_copy( dset ) ;
2560
2561
2562 /*----- Record history of dataset -----*/
2563
2564 sprintf (label, "Output prefix: %s", filename);
2565 if( commandline != NULL )
2566 tross_multi_Append_History( new_dset , commandline,label,NULL ) ;
2567 else
2568 tross_Append_History ( new_dset, label);
2569
2570
2571 output_datum = MRI_short ;
2572 ibuf[0] = output_datum ; ibuf[1] = MRI_short ;
2573
2574
2575 ierror = EDIT_dset_items( new_dset ,
2576 ADN_prefix , filename ,
2577 ADN_label1 , filename ,
2578 ADN_directory_name , option_data->session ,
2579 ADN_self_name , filename ,
2580 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE :
2581 GEN_FUNC_TYPE ,
2582 ADN_func_type , func_type ,
2583 ADN_nvals , FUNC_nvals[func_type] ,
2584 ADN_datum_array , ibuf ,
2585 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
2586 ADN_none ) ;
2587
2588 if( ierror > 0 ){
2589 fprintf(stderr,
2590 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
2591 exit(1) ;
2592 }
2593
2594 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
2595 fprintf(stderr,
2596 "*** Output dataset file %s already exists--cannot continue!\a\n",
2597 new_dset->dblk->diskptr->header_name ) ;
2598 exit(1) ;
2599 }
2600
2601 /*----- deleting exemplar dataset -----*/
2602 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
2603
2604
2605 /*----- allocate memory for output data -----*/
2606 vdif = (void *) malloc( mri_datum_size(output_datum) * nxyz ) ;
2607 tsp = (short *) malloc( sizeof(short) * nxyz ) ;
2608
2609 /*----- attach bricks to new data set -----*/
2610 mri_fix_data_pointer (vdif, DSET_BRICK(new_dset,0));
2611 mri_fix_data_pointer (tsp, DSET_BRICK(new_dset,1));
2612
2613
2614 /*----- convert data type to output specification -----*/
2615 fimfac = EDIT_coerce_autoscale_new (nxyz,
2616 MRI_float, ffim,
2617 output_datum, vdif);
2618 if (fimfac != 0.0) fimfac = 1.0 / fimfac;
2619
2620
2621 top_ss = 32700;
2622
2623 if (func_type == FUNC_THR_TYPE) /* threshold */
2624 {
2625 func_scale_short = FUNC_THR_SCALE_SHORT;
2626 bot_ss = 0;
2627 }
2628 else if (func_type == FUNC_TT_TYPE) /* t-statistic */
2629 {
2630 func_scale_short = FUNC_TT_SCALE_SHORT;
2631 bot_ss = -top_ss;
2632 }
2633 else if (func_type == FUNC_FT_TYPE) /* F-statistic */
2634 {
2635 func_scale_short = FUNC_FT_SCALE_SHORT;
2636 bot_ss = 0;
2637 }
2638 else
2639 RA_error ("Illegal ouput dataset function type");
2640
2641 top = top_ss / func_scale_short;
2642 bot = bot_ss / func_scale_short;
2643
2644
2645 for (ii = 0; ii < nxyz; ii++)
2646 {
2647 if (ftr[ii] > top)
2648 tsp[ii] = top_ss;
2649 else if (ftr[ii] < bot)
2650 tsp[ii] = bot_ss;
2651 else
2652 tsp[ii] = (short) (func_scale_short * ftr[ii] + 0.5);
2653 }
2654
2655
2656 /*----- write afni data set -----*/
2657 if (func_type == FUNC_THR_TYPE) /* threshold */
2658 printf("++ Writing `fith' dataset ");
2659 else if (func_type == FUNC_TT_TYPE) /* t-statistic */
2660 printf("++ Writing `fitt' dataset ");
2661 else if (func_type == FUNC_FT_TYPE) /* F-statistic */
2662 printf("++ Writing `fift' dataset ");
2663
2664 printf("into %s\n", DSET_BRIKNAME(new_dset) ) ;
2665
2666 fbuf[0] = numdof;
2667 fbuf[1] = dendof;
2668 for( ii=2 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
2669 (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
2670
2671 fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
2672 fbuf[1] = 1.0 / func_scale_short ;
2673 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
2674
2675 if( !AFNI_noenv("AFNI_AUTOMATIC_FDR") )
2676 { int ii = THD_create_all_fdrcurves( new_dset ) ;
2677 if( ii > 0 ) ININFO_message("created %d FDR curves in header",ii) ;
2678 }
2679
2680 THD_load_statistics( new_dset ) ;
2681 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
2682
2683
2684 /*----- deallocate memory -----*/
2685 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
2686
2687 }
2688
2689
2690 /*---------------------------------------------------------------------------*/
2691 /*
2692 Routine to write one bucket data set.
2693 */
2694
write_bucket_data(matrix xdata,model * regmodel,RA_options * option_data)2695 void write_bucket_data
2696 (
2697 matrix xdata, /* independent variable matrix */
2698 model * regmodel, /* linear regression model */
2699 RA_options * option_data /* user input options */
2700 )
2701
2702 {
2703 const float EPSILON = 1.0e-10;
2704 int p; /* number of parameters in full model */
2705 int q; /* number of parameters in reduced model */
2706 int n; /* number of data points */
2707 int nxyz; /* number of voxels */
2708 THD_3dim_dataset * old_dset = NULL; /* prototype dataset */
2709 THD_3dim_dataset * new_dset = NULL; /* output bucket dataset */
2710 char * output_prefix; /* prefix name for bucket dataset */
2711 char * output_session; /* directory for bucket dataset */
2712 int nbricks, ib; /* number of sub-bricks in bucket dataset */
2713 short ** bar = NULL; /* bar[ib] points to data for sub-brick #ib */
2714 float ** far = NULL;
2715 float factor; /* factor is new scale factor for sub-brick #ib */
2716 int brick_type; /* indicates statistical type of sub-brick */
2717 int brick_coef; /* regression coefficient index for sub-brick */
2718 char * brick_label; /* character string label for sub-brick */
2719 int ierror; /* number of errors in editing data */
2720 char filename[MAX_NAME_LENGTH]; /* name for temporary data file */
2721 int piece_size; /* number of voxels in dataset piece */
2722 int num_pieces; /* dataset is divided into this many pieces */
2723 float * volume = NULL; /* volume of floating point data */
2724 char label[80]; /* label for output file history */
2725
2726 /*----- initialize local variables -----*/
2727 p = regmodel->p;
2728 q = regmodel->q;
2729 n = xdata.rows;
2730 nxyz = option_data->nxyz;
2731 piece_size = option_data->piece_size;
2732 num_pieces = option_data->num_pieces;
2733 nbricks = option_data->numbricks;
2734 output_prefix = option_data->bucket_filename;
2735 output_session = option_data->session;
2736
2737
2738 /*----- allocate memory -----*/
2739 volume = (float *) malloc (sizeof(float) * nxyz);
2740 MTEST (volume);
2741 if (option_data->datum == MRI_float) {
2742 far = (float **) malloc (sizeof(float *) * nbricks);
2743 MTEST (far);
2744 } else {
2745 bar = (short **) malloc (sizeof(short *) * nbricks);
2746 MTEST (bar);
2747 }
2748
2749 /*----- read first dataset -----*/
2750 old_dset = THD_open_dataset (option_data->first_dataset) ;
2751
2752
2753 /*-- make an empty copy of this dataset, for eventual output --*/
2754 new_dset = EDIT_empty_copy (old_dset);
2755
2756
2757 /*----- Record history of dataset -----*/
2758 if( commandline != NULL )
2759 tross_Append_History( new_dset , commandline ) ;
2760 sprintf (label, "Output prefix: %s", output_prefix);
2761 tross_Append_History ( new_dset, label);
2762
2763
2764 /*----- Modify some structural properties. Note that the nbricks
2765 just make empty sub-bricks, without any data attached. -----*/
2766 ierror = EDIT_dset_items (new_dset,
2767 ADN_prefix, output_prefix,
2768 ADN_directory_name, output_session,
2769 ADN_type, HEAD_FUNC_TYPE,
2770 ADN_func_type, FUNC_BUCK_TYPE,
2771 ADN_ntt, 0, /* no time */
2772 ADN_nvals, nbricks,
2773 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
2774 ADN_none ) ;
2775
2776 if( ierror > 0 )
2777 {
2778 fprintf(stderr,
2779 "*** %d errors in attempting to create output dataset!\n",
2780 ierror);
2781 exit(1);
2782 }
2783
2784 if (THD_is_file(DSET_HEADNAME(new_dset)))
2785 {
2786 fprintf(stderr,
2787 "*** Output dataset file %s already exists--cannot continue!\n",
2788 DSET_HEADNAME(new_dset));
2789 exit(1);
2790 }
2791
2792
2793 /*----- deleting exemplar dataset -----*/
2794 THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ;
2795
2796
2797 /*----- loop over new sub-brick index, attach data array with
2798 EDIT_substitute_brick then put some strings into the labels and
2799 keywords, and modify the sub-brick scaling factor -----*/
2800 for (ib = 0; ib < nbricks; ib++)
2801 {
2802 /*----- get information about this sub-brick -----*/
2803 brick_type = option_data->brick_type[ib];
2804 brick_coef = option_data->brick_coef[ib];
2805 brick_label = option_data->brick_label[ib];
2806
2807 if (brick_type == FUNC_FIM_TYPE)
2808 {
2809 sprintf (filename, "coef.%d", brick_coef);
2810 }
2811 else if (brick_type == FUNC_THR_TYPE)
2812 {
2813 sprintf (filename, "rsqr");
2814 }
2815 else if (brick_type == FUNC_TT_TYPE)
2816 {
2817 sprintf (filename, "tcoef.%d", brick_coef);
2818 EDIT_BRICK_TO_FITT (new_dset, ib, n-p);
2819 }
2820 else if (brick_type == FUNC_FT_TYPE)
2821 {
2822 sprintf (filename, "freg");
2823 EDIT_BRICK_TO_FIFT (new_dset, ib, p-q, n-p);
2824 }
2825
2826 /*----- allocate memory for output sub-brick -----*/
2827 read_volume (filename, volume, nxyz, piece_size, num_pieces);
2828 delete_volume (filename, nxyz, piece_size, num_pieces);
2829
2830 if (option_data->datum == MRI_float) {
2831 far[ib] = (float *) malloc (sizeof(float) * nxyz);
2832 MTEST (far[ib]);
2833 memcpy((void *)far[ib], (void *)volume, sizeof(float) * nxyz);
2834 /*----- attach far[ib] to be sub-brick #ib -----*/
2835 EDIT_substitute_brick (new_dset, ib, MRI_float, far[ib]);
2836 } else {
2837 bar[ib] = (short *) malloc (sizeof(short) * nxyz);
2838 MTEST (bar[ib]);
2839 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
2840 MRI_short, bar[ib]);
2841
2842 if (factor < EPSILON) factor = 0.0;
2843 else factor = 1.0 / factor;
2844 EDIT_BRICK_FACTOR (new_dset, ib, factor);
2845 /*----- attach bar[ib] to be sub-brick #ib -----*/
2846 EDIT_substitute_brick (new_dset, ib, MRI_short, bar[ib]);
2847 }
2848
2849 /*----- edit the sub-brick -----*/
2850 EDIT_BRICK_LABEL (new_dset, ib, brick_label);
2851
2852
2853
2854 }
2855
2856
2857 /*----- write bucket data set -----*/
2858
2859 INFO_message("Writing bucket dataset: %s", DSET_BRIKNAME(new_dset)) ;
2860
2861 if( !AFNI_noenv("AFNI_AUTOMATIC_FDR") )
2862 { int ii = THD_create_all_fdrcurves( new_dset ) ;
2863 if( ii > 0 ) ININFO_message("created %d FDR curves in header",ii) ;
2864 }
2865
2866 THD_load_statistics (new_dset);
2867 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
2868
2869
2870 /*----- deallocate memory -----*/
2871 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
2872 free (volume); volume = NULL;
2873
2874 }
2875
2876
2877 /*---------------------------------------------------------------------------*/
2878 /*
2879 Write out the user requested output files.
2880 */
2881
output_results(matrix xdata,model * regmodel,RA_options * option_data)2882 void output_results
2883 (
2884 matrix xdata, /* independent variable matrix */
2885 model * regmodel, /* linear regression model */
2886 RA_options * option_data /* user input options */
2887 )
2888
2889 {
2890 int p; /* number of parameters in full model */
2891 int q; /* number of parameters in reduced model */
2892 int n; /* number of data points */
2893 int nxyz; /* number of voxels */
2894 int ip, ix; /* parameter index */
2895 int numdof, dendof; /* numerator and denominator degrees of freedom */
2896 int piece_size; /* number of voxels in dataset piece */
2897 int num_pieces; /* dataset is divided into this many pieces */
2898 float * volume1 = NULL; /* volume data for 1st sub-brick */
2899 float * volume2 = NULL; /* volume data for 2nd sub-brick */
2900 char filename[MAX_NAME_LENGTH]; /* name for temporary data file */
2901
2902
2903 /*----- initialize local variables -----*/
2904 p = regmodel->p;
2905 q = regmodel->q;
2906 n = xdata.rows;
2907 nxyz = option_data->nxyz;
2908 piece_size = option_data->piece_size;
2909 num_pieces = option_data->num_pieces;
2910
2911
2912 /*----- allocate memory space for volume data -----*/
2913 volume1 = (float *) malloc (sizeof(float) * nxyz);
2914 MTEST (volume1);
2915 volume2 = (float *) malloc (sizeof(float) * nxyz);
2916 MTEST (volume2);
2917
2918
2919 /*----- write t-statistics data files -----*/
2920 if (option_data->numt > 0)
2921 {
2922 numdof = n - p;
2923 dendof = 0;
2924
2925 for (ip = 0; ip < p; ip++)
2926 {
2927 ix = regmodel->flist[ip];
2928
2929 if (option_data->tcoef_filename[ix] != NULL)
2930 {
2931 sprintf (filename, "coef.%d", ix);
2932 read_volume (filename, volume1, nxyz, piece_size, num_pieces);
2933
2934 sprintf (filename, "tcoef.%d", ix);
2935 read_volume (filename, volume2, nxyz, piece_size, num_pieces);
2936
2937 write_afni_data (option_data,
2938 option_data->tcoef_filename[ix],
2939 volume1, volume2,
2940 FUNC_TT_TYPE, numdof, dendof);
2941 }
2942 }
2943 }
2944
2945
2946 /*----- write R^2 data files -----*/
2947 if (option_data->numr > 0)
2948 {
2949 strcpy (filename, "rsqr");
2950 read_volume (filename, volume2, nxyz, piece_size, num_pieces);
2951
2952 for (ip = 0; ip < p; ip++)
2953 {
2954 ix = regmodel->flist[ip];
2955
2956 if (option_data->rcoef_filename[ix] != NULL)
2957 {
2958 sprintf (filename, "coef.%d", ix);
2959 read_volume (filename, volume1, nxyz, piece_size, num_pieces);
2960
2961 write_afni_data (option_data,
2962 option_data->rcoef_filename[ix],
2963 volume1, volume2,
2964 FUNC_THR_TYPE, 0, 0);
2965 }
2966 }
2967 }
2968
2969
2970 /*----- write F-statistics data files -----*/
2971 if (option_data->numf > 0)
2972 {
2973 strcpy (filename, "freg");
2974 read_volume (filename, volume2, nxyz, piece_size, num_pieces);
2975
2976 numdof = p - q;
2977 dendof = n - p;
2978
2979 for (ip = 0; ip < p; ip++)
2980 {
2981 ix = regmodel->flist[ip];
2982
2983 if (option_data->fcoef_filename[ix] != NULL)
2984 {
2985 sprintf (filename, "coef.%d", ix);
2986 read_volume (filename, volume1, nxyz, piece_size, num_pieces);
2987
2988 write_afni_data (option_data,
2989 option_data->fcoef_filename[ix],
2990 volume1, volume2,
2991 FUNC_FT_TYPE, numdof, dendof);
2992 }
2993 }
2994 }
2995
2996
2997 /*----- deallocate memory space for volume data -----*/
2998 free (volume1); volume1 = NULL;
2999 free (volume2); volume2 = NULL;
3000
3001
3002 /*----- write the bucket dataset -----*/
3003 if (option_data->numbricks > 0)
3004 write_bucket_data (xdata, regmodel, option_data);
3005
3006
3007 }
3008
3009
3010 /*---------------------------------------------------------------------------*/
3011 /*
3012 Deallocate memory and delete any remaining temporary data files.
3013 */
3014
terminate_program(matrix * xdata,model * regmodel,RA_options * option_data)3015 void terminate_program
3016 (
3017 matrix * xdata, /* independent variable matrix */
3018 model * regmodel, /* linear regression model */
3019 RA_options * option_data /* user input options */
3020 )
3021
3022 {
3023 int p; /* number of parameters in full model */
3024 int ip, ix; /* parameter index */
3025 int ib; /* sub-brick index */
3026 int nxyz; /* number of voxels */
3027 int piece_size; /* number of voxels in dataset piece */
3028 int num_pieces; /* dataset is divided into this many pieces */
3029 char filename[MAX_NAME_LENGTH]; /* name for temporary data file */
3030
3031
3032 /*----- initialize local variables -----*/
3033 p = regmodel->p;
3034 nxyz = option_data->nxyz;
3035 piece_size = option_data->piece_size;
3036 num_pieces = option_data->num_pieces;
3037
3038
3039 /*----- delete F-statistics data files -----*/
3040 if (option_data->numf > 0)
3041 {
3042 strcpy (filename, "freg");
3043 delete_volume (filename, nxyz, piece_size, num_pieces);
3044 }
3045
3046
3047 /*----- delete R^2 data files -----*/
3048 if (option_data->numr > 0)
3049 {
3050 strcpy (filename, "rsqr");
3051 delete_volume (filename, nxyz, piece_size, num_pieces);
3052 }
3053
3054
3055 /*----- delete t-statistics data files -----*/
3056 if (option_data->numt > 0)
3057 {
3058 for (ip = 0; ip < p; ip++)
3059 {
3060 ix = regmodel->flist[ip];
3061 sprintf (filename, "tcoef.%d", ix);
3062 delete_volume (filename, nxyz, piece_size, num_pieces);
3063 }
3064 }
3065
3066
3067 /*----- delete regression coefficients data files -----*/
3068 if (option_data->numc > 0)
3069 {
3070 for (ip = 0; ip < p; ip++)
3071 {
3072 ix = regmodel->flist[ip];
3073 sprintf (filename, "coef.%d", ix);
3074 delete_volume (filename, nxyz, piece_size, num_pieces);
3075 }
3076 }
3077
3078
3079 /*----- dispose of data matrix -----*/
3080 matrix_destroy (xdata);
3081
3082
3083 /*----- deallocate memory -----*/
3084 if (option_data->yname != NULL)
3085 {
3086 for (ip = 0; ip < MAX_OBSERVATIONS; ip++)
3087 {
3088 if (option_data->yname[ip] != NULL)
3089 {
3090 free (option_data->yname[ip]);
3091 option_data->yname[ip] = NULL;
3092 }
3093 }
3094 free (option_data->yname);
3095 option_data->yname = NULL;
3096 }
3097
3098 if (option_data->fcoef_filename != NULL)
3099 {
3100 for (ip = 0; ip < MAX_XVARS; ip++)
3101 {
3102 if (option_data->fcoef_filename[ip] != NULL)
3103 {
3104 free (option_data->fcoef_filename[ip]);
3105 option_data->fcoef_filename[ip] = NULL;
3106 }
3107 }
3108 free (option_data->fcoef_filename);
3109 option_data->fcoef_filename = NULL;
3110 }
3111
3112 if (option_data->rcoef_filename != NULL)
3113 {
3114 for (ip = 0; ip < MAX_XVARS; ip++)
3115 {
3116 if (option_data->rcoef_filename[ip] != NULL)
3117 {
3118 free (option_data->rcoef_filename[ip]);
3119 option_data->rcoef_filename[ip] = NULL;
3120 }
3121 }
3122 free (option_data->rcoef_filename);
3123 option_data->rcoef_filename = NULL;
3124 }
3125
3126 if (option_data->tcoef_filename != NULL)
3127 {
3128 for (ip = 0; ip < MAX_XVARS; ip++)
3129 {
3130 if (option_data->tcoef_filename[ip] != NULL)
3131 {
3132 free (option_data->tcoef_filename[ip]);
3133 option_data->tcoef_filename[ip] = NULL;
3134 }
3135 }
3136 free (option_data->tcoef_filename);
3137 option_data->tcoef_filename = NULL;
3138 }
3139
3140 if (option_data->levels != NULL)
3141 {
3142 free (option_data->levels);
3143 option_data->levels = NULL;
3144 }
3145
3146 if (option_data->counts != NULL)
3147 {
3148 free (option_data->counts);
3149 option_data->counts = NULL;
3150 }
3151
3152
3153 if (option_data->bucket_filename != NULL)
3154 {
3155 free (option_data->bucket_filename);
3156 option_data->bucket_filename = NULL;
3157 }
3158
3159 if (option_data->brick_type != NULL)
3160 {
3161 free (option_data->brick_type);
3162 option_data->brick_type = NULL;
3163 }
3164
3165 if (option_data->brick_coef != NULL)
3166 {
3167 free (option_data->brick_coef);
3168 option_data->brick_coef = NULL;
3169 }
3170
3171 if (option_data->brick_label != NULL)
3172 {
3173 for (ib = 0; ib < option_data->numbricks; ib++)
3174 {
3175 if (option_data->brick_label[ib] != NULL)
3176 {
3177 free (option_data->brick_label[ib]);
3178 option_data->brick_label[ib] = NULL;
3179 }
3180 }
3181 free (option_data->brick_label);
3182 option_data->brick_label = NULL;
3183 }
3184
3185 }
3186
3187
3188 /*---------------------------------------------------------------------------*/
3189 /*
3190 Multiple linear regression analysis (3dRegAna)
3191 */
3192
main(int argc,char ** argv)3193 int main
3194 (
3195 int argc, /* number of input arguments */
3196 char ** argv /* array of input arguments */
3197 )
3198
3199 {
3200 RA_options option_data; /* user input options */
3201 matrix xdata; /* independent variable matrix */
3202 model regmodel; /* linear regression model */
3203 int piece_size; /* number of voxels in dataset piece */
3204
3205
3206 /*----- Identify software -----*/
3207 #if 0
3208 printf ("\n\n");
3209 printf ("Program: %s \n", PROGRAM_NAME);
3210 printf ("Author: %s \n", PROGRAM_AUTHOR);
3211 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
3212 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
3213 printf ("\n");
3214 #endif
3215
3216 UNIQ_idcode_fill(SUFFIX+strlen(SUFFIX)) ; /* 27 Apr 2012 */
3217
3218 /*-- 22 Feb 1999: addto the arglist, if user wants to --*/
3219
3220 PRINT_VERSION("3dRegAna") ; AUTHOR(PROGRAM_AUTHOR);
3221 mainENTRY("3dRegAna main") ; machdep() ;
3222
3223 { int new_argc ; char ** new_argv ;
3224 addto_args( argc , argv , &new_argc , &new_argv ) ;
3225 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
3226 }
3227
3228 /*----- program initialization -----*/
3229 initialize_program (argc, argv, &xdata, ®model, &option_data);
3230
3231
3232 /*----- perform regression analysis -----*/
3233 calculate_results (xdata, ®model, &option_data);
3234
3235
3236 /*----- write requested output files -----*/
3237 output_results (xdata, ®model, &option_data);
3238
3239
3240 /*----- end of program -----*/
3241 terminate_program (&xdata, ®model, &option_data);
3242
3243 exit(0);
3244 }
3245
3246
3247 /*---------------------------------------------------------------------------*/
3248