00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059 #define PROGRAM_NAME "3dRegAna"
00060 #define PROGRAM_AUTHOR "B. Douglas Ward"
00061 #define PROGRAM_INITIAL "10 Oct 1997"
00062 #define PROGRAM_LATEST "02 Dec 2002"
00063
00064
00065
00066 #define SUFFIX ".3dregana"
00067
00068 #include <stdio.h>
00069 #include <math.h>
00070 #include "mrilib.h"
00071 #include "matrix.h"
00072
00073
00074 #include "RegAna.c"
00075
00076
00077
00078
00079 #ifdef HP
00080 # define DF "bdf ."
00081 #endif
00082
00083
00084 #ifdef OSF1
00085 # define DF "df -k ."
00086 #endif
00087
00088
00089 #ifdef SGI
00090 # define DF "df -k ."
00091 #endif
00092
00093
00094 #if defined(SOLARIS) || defined(SUN)
00095 # define DF "df -k"
00096 #endif
00097
00098
00099 #ifdef RS6000
00100 #endif
00101
00102
00103 #ifdef LINUX
00104 # define DF "df -k ."
00105 #endif
00106
00107
00108 #ifndef DF
00109 # define DF "df"
00110 #endif
00111
00112
00113
00114
00115 #define MAX_XVARS 101
00116 #define MAX_OBSERVATIONS 1000
00117 #define MAX_NAME_LENGTH THD_MAX_NAME
00118 #define MEGA 1048576
00119
00120 static char * commandline = NULL ;
00121
00122
00123 typedef struct model
00124 {
00125 int p;
00126 int flist[MAX_XVARS];
00127 int q;
00128 int rlist[MAX_XVARS];
00129 } model;
00130
00131
00132 typedef struct RA_options
00133 {
00134 int datum;
00135 char session[MAX_NAME_LENGTH];
00136
00137 char ** yname;
00138 char * first_dataset;
00139
00140 int nx, ny, nz;
00141 int nxyz;
00142
00143 int diskspace;
00144
00145 int workmem;
00146 int piece_size;
00147 int num_pieces;
00148
00149 float rms_min;
00150 float fdisp;
00151
00152 int * levels;
00153 int * counts;
00154 int c;
00155 float flofmax;
00156
00157 int numf;
00158 int numr;
00159 int numt;
00160 int numc;
00161
00162 char ** fcoef_filename;
00163 char ** rcoef_filename;
00164 char ** tcoef_filename;
00165
00166 int numfiles;
00167
00168 char * bucket_filename;
00169 int numbricks;
00170 int * brick_type;
00171 int * brick_coef;
00172 char ** brick_label;
00173
00174 } RA_options;
00175
00176
00177
00178
00179
00180
00181
00182 void RA_error (char * message)
00183 {
00184 fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00185 exit(1);
00186 }
00187
00188
00189
00190
00191
00192
00193
00194 void display_help_menu()
00195 {
00196 printf
00197 (
00198 "This program performs multiple linear regression analysis. \n\n"
00199 "Usage: \n"
00200 "3dRegAna \n"
00201 "-rows n number of input datasets \n"
00202 "-cols m number of X variables \n"
00203 "-xydata X11 X12 ... X1m filename X variables and Y observations \n"
00204 " . . \n"
00205 " . . \n"
00206 " . . \n"
00207 "-xydata Xn1 Xn2 ... Xnm filename X variables and Y observations \n"
00208 " \n"
00209 "-model i1 ... iq : j1 ... jr definition of linear regression model; \n"
00210 " reduced model: \n"
00211 " Y = f(Xj1,...,Xjr) \n"
00212 " full model: \n"
00213 " Y = f(Xj1,...,Xjr,Xi1,...,Xiq) \n"
00214 " \n"
00215 "[-diskspace] print out disk space required for program execution\n"
00216 "[-workmem mega] number of megabytes of RAM to use for statistical \n"
00217 " workspace (default = 12) \n"
00218 "[-rmsmin r] r = minimum rms error to reject constant model \n"
00219 "[-fdisp fval] display (to screen) results for those voxels \n"
00220 " whose F-statistic is > fval \n"
00221 " \n"
00222 "[-flof alpha] alpha = minimum p value for F due to lack of fit \n"
00223 " \n"
00224 " \n"
00225 "The following commands generate individual AFNI 2 sub-brick datasets: \n"
00226 " \n"
00227 "[-fcoef k prefixname] estimate of kth regression coefficient \n"
00228 " along with F-test for the regression \n"
00229 " is written to AFNI `fift' dataset \n"
00230 "[-rcoef k prefixname] estimate of kth regression coefficient \n"
00231 " along with coef. of mult. deter. R^2 \n"
00232 " is written to AFNI `fith' dataset \n"
00233 "[-tcoef k prefixname] estimate of kth regression coefficient \n"
00234 " along with t-test for the coefficient \n"
00235 " is written to AFNI `fitt' dataset \n"
00236 " \n"
00237 " \n"
00238 "The following commands generate one AFNI 'bucket' type dataset: \n"
00239 " \n"
00240 "[-bucket n prefixname] create one AFNI 'bucket' dataset having \n"
00241 " n sub-bricks; n=0 creates default output;\n"
00242 " output 'bucket' is written to prefixname \n"
00243 "The mth sub-brick will contain: \n"
00244 "[-brick m coef k label] kth parameter regression coefficient \n"
00245 "[-brick m fstat label] F-stat for significance of regression \n"
00246 "[-brick m rstat label] coefficient of multiple determination R^2 \n"
00247 "[-brick m tstat k label] t-stat for kth regression coefficient \n"
00248 "\n" );
00249
00250 printf
00251 (
00252 "\n"
00253 "N.B.: For this program, the user must specify 1 and only 1 sub-brick \n"
00254 " with each -xydata command. That is, if an input dataset contains\n"
00255 " more than 1 sub-brick, a sub-brick selector must be used, e.g.: \n"
00256 " -xydata 2.17 4.59 7.18 'fred+orig[3]' \n"
00257 );
00258
00259 printf("\n" MASTER_SHORTHELP_STRING ) ;
00260
00261 exit(0);
00262 }
00263
00264
00265
00266
00267
00268
00269 #define DOPEN(ds,name) \
00270 do{ int pv ; (ds) = THD_open_dataset((name)) ; \
00271 if( !ISVALID_3DIM_DATASET((ds)) ){ \
00272 fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \
00273 if( (ds)->daxes->nxx!=nx || (ds)->daxes->nyy!=ny || \
00274 (ds)->daxes->nzz!=nz ){ \
00275 fprintf(stderr,"*** Axes mismatch: %s\n",(name)) ; exit(1) ; } \
00276 if( DSET_NVALS((ds)) != 1 ){ \
00277 fprintf(stderr,"*** Must specify 1 sub-brick: %s\n",(name));exit(1);}\
00278 THD_load_datablock( (ds)->dblk ) ; \
00279 pv = DSET_PRINCIPAL_VALUE((ds)) ; \
00280 if( DSET_ARRAY((ds),pv) == NULL ){ \
00281 fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); } \
00282 if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){ \
00283 fprintf(stderr,"*** Can't use complex data: %s\n",(name)); exit(1); \
00284 } \
00285 break ; } while (0)
00286
00287
00288
00289
00290
00291
00292 #define SUB_POINTER(ds,vv,ind,ptr) \
00293 do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){ \
00294 default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1); \
00295 case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ; \
00296 (ptr) = (void *)( fim + (ind) ) ; \
00297 } break ; \
00298 case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ; \
00299 (ptr) = (void *)( fim + (ind) ) ; \
00300 } break ; \
00301 case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ; \
00302 (ptr) = (void *)( fim + (ind) ) ; \
00303 } break ; } break ; } while(0)
00304
00305
00306
00307
00308
00309
00310
00311 void get_dimensions
00312 (
00313 RA_options * option_data
00314 )
00315
00316 {
00317
00318 THD_3dim_dataset * dset=NULL;
00319
00320
00321
00322 dset = THD_open_dataset( option_data->first_dataset ) ;
00323 if( ! ISVALID_3DIM_DATASET(dset) ){
00324 fprintf(stderr,"*** Unable to open dataset file %s\n",
00325 option_data->first_dataset);
00326 exit(1) ;
00327 }
00328
00329
00330 option_data->nx = dset->daxes->nxx ;
00331 option_data->ny = dset->daxes->nyy ;
00332 option_data->nz = dset->daxes->nzz ;
00333 option_data->nxyz = option_data->nx * option_data->ny * option_data->nz ;
00334
00335 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00336
00337 }
00338
00339
00340
00341
00342
00343
00344
00345
00346 void read_afni_data
00347 (
00348 RA_options * option_data,
00349 char * filename,
00350 int piece_len,
00351 int fim_offset,
00352 float * ffim
00353 )
00354
00355 {
00356 int iv;
00357 THD_3dim_dataset * dset=NULL;
00358 void * vfim = NULL;
00359 int nx, ny, nz, nxyz;
00360
00361 nx = option_data->nx;
00362 ny = option_data->ny;
00363 nz = option_data->nz;
00364 nxyz = option_data->nxyz;
00365
00366
00367
00368 DOPEN (dset, filename) ;
00369 iv = DSET_PRINCIPAL_VALUE(dset) ;
00370
00371
00372 SUB_POINTER (dset, iv, fim_offset, vfim) ;
00373 EDIT_coerce_scale_type (piece_len, DSET_BRICK_FACTOR(dset,iv),
00374 DSET_BRICK_TYPE(dset,iv), vfim,
00375 MRI_float ,ffim ) ;
00376
00377 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00378 }
00379
00380
00381
00382
00383
00384
00385
00386 void check_piece
00387 (
00388 char * filename
00389 )
00390
00391 {
00392 FILE * far = NULL;
00393 char sfilename[MAX_NAME_LENGTH];
00394 char message[MAX_NAME_LENGTH];
00395
00396
00397
00398 strcpy (sfilename, filename);
00399 strcat (sfilename, SUFFIX);
00400 far = fopen (sfilename, "r");
00401 if (far != NULL)
00402 {
00403 sprintf (message, "temporary file %s already exists. ", sfilename);
00404 RA_error (message);
00405 }
00406
00407 }
00408
00409
00410
00411
00412
00413
00414
00415 void read_piece
00416 (
00417 char * filename,
00418 float * fin,
00419 int size
00420 )
00421
00422 {
00423 char sfilename[MAX_NAME_LENGTH];
00424 char message[MAX_NAME_LENGTH];
00425 FILE * far = NULL;
00426 int count;
00427
00428
00429
00430 strcpy (sfilename, filename);
00431 strcat (sfilename, SUFFIX);
00432
00433
00434 far = fopen (sfilename, "r");
00435 if (far == NULL)
00436 {
00437 sprintf (message, "Unable to open temporary file %s", sfilename);
00438 RA_error (message);
00439 }
00440
00441
00442 count = fread (fin, sizeof(float), size, far);
00443 fclose (far);
00444
00445
00446 if (count != size)
00447 {
00448 sprintf (message, "Error in reading temporary file %s", sfilename);
00449 RA_error (message);
00450 }
00451 }
00452
00453
00454
00455
00456
00457
00458
00459
00460 void write_piece
00461 (
00462 char * filename,
00463 float * fout,
00464 int size
00465 )
00466
00467 {
00468 char sfilename[MAX_NAME_LENGTH];
00469 char message[MAX_NAME_LENGTH];
00470 FILE * far = NULL;
00471 int count;
00472
00473
00474
00475 strcpy (sfilename, filename);
00476 strcat (sfilename, SUFFIX);
00477
00478
00479 far = fopen (sfilename, "r");
00480 if (far != NULL)
00481 {
00482 sprintf (message, "Temporary file %s already exists. ", sfilename);
00483 RA_error (message);
00484 }
00485
00486
00487 far = fopen (sfilename, "w");
00488 if (far == NULL)
00489 {
00490 sprintf (message, "Unable to write temporary file %s ", sfilename);
00491 RA_error (message);
00492 }
00493
00494
00495 count = fwrite (fout, sizeof(float), size, far);
00496 fclose (far);
00497
00498
00499 if (count != size)
00500 {
00501 sprintf (message, "Error in writing temporary file %s ", sfilename);
00502 RA_error (message);
00503 }
00504
00505 }
00506
00507
00508
00509
00510
00511
00512
00513 void delete_piece
00514 (
00515 char * filename
00516 )
00517
00518 {
00519 char sfilename[MAX_NAME_LENGTH];
00520
00521
00522 strcpy (sfilename, filename);
00523 strcat (sfilename, SUFFIX);
00524
00525
00526 remove (sfilename);
00527
00528 }
00529
00530
00531
00532
00533
00534
00535
00536 void allocate_pieces
00537 (
00538 matrix xdata,
00539 model * regmodel,
00540 RA_options * option_data,
00541
00542 float *** yfimar,
00543 float ** freg_piece,
00544 float ** rsqr_piece,
00545 float *** coef_piece,
00546 float *** tcoef_piece
00547 )
00548
00549 {
00550 int n;
00551 int p;
00552 int * flist = NULL;
00553 int i;
00554 int ip;
00555 int ix;
00556 int piece_size;
00557 int ib;
00558 int nbricks;
00559
00560
00561
00562 n = xdata.rows;
00563 p = regmodel->p;
00564 flist = regmodel->flist;
00565 piece_size = option_data->piece_size;
00566 nbricks = option_data->numbricks;
00567
00568
00569
00570 *yfimar = (float **) malloc (sizeof(float *) * n);
00571 MTEST (*yfimar);
00572 for (i = 0; i < n; i++)
00573 {
00574 (*yfimar)[i] = (float *) malloc(sizeof(float) * piece_size);
00575 MTEST ((*yfimar)[i]);
00576 }
00577
00578
00579
00580 for (ip = 0; ip < p; ip++)
00581 {
00582 ix = flist[ip];
00583
00584 if (option_data->fcoef_filename[ix] != NULL)
00585 {
00586 if (*freg_piece == NULL)
00587 {
00588 *freg_piece = (float *) malloc (sizeof(float) * piece_size);
00589 MTEST (*freg_piece);
00590 }
00591 }
00592 }
00593
00594 for (ib = 0; ib < nbricks; ib++)
00595 {
00596 if (option_data->brick_type[ib] == FUNC_FT_TYPE)
00597 {
00598 if (*freg_piece == NULL)
00599 {
00600 *freg_piece = (float *) malloc (sizeof(float) * piece_size);
00601 MTEST (*freg_piece);
00602 }
00603 }
00604 }
00605
00606
00607
00608 for (ip = 0; ip < p; ip++)
00609 {
00610 ix = flist[ip];
00611
00612 if (option_data->rcoef_filename[ix] != NULL)
00613 {
00614 if (*rsqr_piece == NULL)
00615 {
00616 *rsqr_piece = (float *) malloc (sizeof(float) * piece_size);
00617 MTEST (*rsqr_piece);
00618 }
00619 }
00620 }
00621
00622 for (ib = 0; ib < nbricks; ib++)
00623 {
00624 if (option_data->brick_type[ib] == FUNC_THR_TYPE)
00625 {
00626 if (*rsqr_piece == NULL)
00627 {
00628 *rsqr_piece = (float *) malloc (sizeof(float) * piece_size);
00629 MTEST (*rsqr_piece);
00630 }
00631 }
00632 }
00633
00634
00635
00636 *coef_piece = (float **) malloc (sizeof(float *) * MAX_XVARS);
00637 MTEST (*coef_piece);
00638
00639 *tcoef_piece = (float **) malloc (sizeof(float *) * MAX_XVARS);
00640 MTEST (*tcoef_piece);
00641
00642 for (ix = 0; ix < MAX_XVARS; ix++)
00643 {
00644 (*coef_piece)[ix] = NULL;
00645 (*tcoef_piece)[ix] = NULL;
00646 }
00647
00648 for (ip = 0; ip < p; ip++)
00649 {
00650 ix = flist[ip];
00651
00652 if ( (option_data->fcoef_filename[ix] != NULL)
00653 || (option_data->rcoef_filename[ix] != NULL)
00654 || (option_data->tcoef_filename[ix] != NULL))
00655 {
00656 (*coef_piece)[ix] =
00657 (float *) malloc (sizeof(float) * piece_size);
00658 MTEST ((*coef_piece)[ix]);
00659 }
00660
00661 if (option_data->tcoef_filename[ix] != NULL)
00662 {
00663 (*tcoef_piece)[ix] =
00664 (float *) malloc (sizeof(float) * piece_size);
00665 MTEST ((*tcoef_piece)[ix]);
00666 }
00667 }
00668
00669 for (ib = 0; ib < nbricks; ib++)
00670 {
00671 if (option_data->brick_type[ib] == FUNC_FIM_TYPE)
00672 {
00673 ix = option_data->brick_coef[ib];
00674 if ((*coef_piece)[ix] == NULL)
00675 {
00676 (*coef_piece)[ix] = (float *) malloc (sizeof(float)*piece_size);
00677 MTEST ((*coef_piece)[ix]);
00678 }
00679 }
00680 }
00681
00682 for (ib = 0; ib < nbricks; ib++)
00683 {
00684 if (option_data->brick_type[ib] == FUNC_TT_TYPE)
00685 {
00686 ix = option_data->brick_coef[ib];
00687 if ((*tcoef_piece)[ix] == NULL)
00688 {
00689 (*tcoef_piece)[ix] = (float *) malloc (sizeof(float)*piece_size);
00690 MTEST ((*tcoef_piece)[ix]);
00691 }
00692 }
00693 }
00694
00695 }
00696
00697
00698
00699
00700
00701
00702
00703 void save_pieces
00704 (
00705 int piece,
00706 int piece_len,
00707
00708 float * freg_piece,
00709 float * rsqr_piece,
00710 float ** coef_piece,
00711 float ** tcoef_piece
00712 )
00713
00714 {
00715 int ip;
00716 char filename[MAX_NAME_LENGTH];
00717
00718
00719
00720 if (freg_piece != NULL)
00721 {
00722 sprintf (filename, "freg.p%d", piece);
00723 write_piece (filename, freg_piece, piece_len);
00724 }
00725
00726
00727
00728 if (rsqr_piece != NULL)
00729 {
00730 sprintf (filename, "rsqr.p%d", piece);
00731 write_piece (filename, rsqr_piece, piece_len);
00732 }
00733
00734
00735
00736 if (coef_piece != NULL)
00737 {
00738 for (ip = 0; ip < MAX_XVARS; ip++)
00739 if (coef_piece[ip] != NULL)
00740 {
00741 sprintf (filename, "coef.%d.p%d", ip, piece);
00742 write_piece (filename, coef_piece[ip], piece_len);
00743 }
00744 }
00745
00746
00747
00748 if (tcoef_piece != NULL)
00749 {
00750 for (ip = 0; ip < MAX_XVARS; ip++)
00751 if (tcoef_piece[ip] != NULL)
00752 {
00753 sprintf (filename, "tcoef.%d.p%d", ip, piece);
00754 write_piece (filename, tcoef_piece[ip], piece_len);
00755 }
00756 }
00757
00758 }
00759
00760
00761
00762
00763
00764
00765
00766 void deallocate_pieces
00767 (
00768 int n,
00769 float *** yfimar,
00770 float ** freg_piece,
00771 float ** rsqr_piece,
00772 float *** coef_piece,
00773 float *** tcoef_piece
00774 )
00775
00776 {
00777 int i;
00778 int ip;
00779
00780
00781
00782 if (*yfimar != NULL)
00783 {
00784 for (i = 0; i < n; i++)
00785 {
00786 free ((*yfimar)[i]); (*yfimar)[i] = NULL;
00787 }
00788 free (*yfimar);
00789 *yfimar = NULL;
00790 }
00791
00792
00793 if (*freg_piece != NULL)
00794 {
00795 free (*freg_piece);
00796 *freg_piece = NULL;
00797 }
00798
00799
00800 if (*rsqr_piece != NULL)
00801 {
00802 free (*rsqr_piece);
00803 *rsqr_piece = NULL;
00804 }
00805
00806
00807 if (*coef_piece != NULL)
00808 {
00809 for (ip = 0; ip < MAX_XVARS; ip++)
00810 if ((*coef_piece)[ip] != NULL)
00811 {
00812 free ((*coef_piece)[ip]);
00813 (*coef_piece)[ip] = NULL;
00814 }
00815 free (*coef_piece);
00816 *coef_piece = NULL;
00817 }
00818
00819
00820 if (*tcoef_piece != NULL)
00821 {
00822 for (ip = 0; ip < MAX_XVARS; ip++)
00823 if ((*tcoef_piece)[ip] != NULL)
00824 {
00825 free ((*tcoef_piece)[ip]);
00826 (*tcoef_piece)[ip] = NULL;
00827 }
00828 free (*tcoef_piece);
00829 *tcoef_piece = NULL;
00830 }
00831
00832 }
00833
00834
00835
00836
00837
00838
00839
00840 void check_volume
00841 (
00842 char * filename,
00843 int num_pieces
00844 )
00845
00846 {
00847 int piece;
00848 char sfilename[MAX_NAME_LENGTH];
00849
00850
00851
00852 for (piece = 0; piece < num_pieces; piece++)
00853 {
00854
00855 sprintf (sfilename, "%s.p%d", filename, piece);
00856 check_piece (sfilename);
00857
00858 }
00859
00860 }
00861
00862
00863
00864
00865
00866
00867
00868 void read_volume
00869 (
00870 char * filename,
00871 float * volume,
00872 int nxyz,
00873 int piece_size,
00874 int num_pieces
00875 )
00876
00877 {
00878 int piece;
00879 int piece_len;
00880 int fim_offset;
00881 char sfilename[MAX_NAME_LENGTH];
00882
00883
00884
00885 for (piece = 0; piece < num_pieces; piece++)
00886 {
00887
00888
00889 fim_offset = piece * piece_size;
00890
00891
00892 if (piece < num_pieces-1)
00893 piece_len = piece_size;
00894 else
00895 piece_len = nxyz - fim_offset;
00896
00897
00898 sprintf (sfilename, "%s.p%d", filename, piece);
00899 read_piece (sfilename, volume + fim_offset, piece_len);
00900
00901 }
00902
00903 }
00904
00905
00906
00907
00908
00909
00910
00911 void delete_volume
00912 (
00913 char * filename,
00914 int nxyz,
00915 int piece_size,
00916 int num_pieces
00917 )
00918
00919 {
00920 int piece;
00921 char sfilename[MAX_NAME_LENGTH];
00922
00923
00924
00925 for (piece = 0; piece < num_pieces; piece++)
00926 {
00927
00928
00929 sprintf (sfilename, "%s.p%d", filename, piece);
00930 delete_piece (sfilename);
00931
00932 }
00933
00934 }
00935
00936
00937
00938
00939
00940
00941
00942 void initialize_options
00943 (
00944 model * regmodel,
00945 RA_options * option_data
00946 )
00947
00948 {
00949 int ip;
00950
00951
00952 regmodel->p = 0;
00953 regmodel->q = 0;
00954
00955 option_data->datum = ILLEGAL_TYPE;
00956 strcpy (option_data->session, "./");
00957
00958 option_data->yname = NULL;
00959 option_data->first_dataset = NULL;
00960
00961 option_data->nx = 0;
00962 option_data->ny = 0;
00963 option_data->nz = 0;
00964 option_data->nxyz = 0;
00965
00966 option_data->diskspace = 0;
00967 option_data->workmem = 12;
00968 option_data->piece_size = 0;
00969 option_data->num_pieces = 0;
00970
00971 option_data->rms_min = 0.0;
00972 option_data->fdisp = -1.0;
00973
00974 option_data->levels = NULL;
00975 option_data->counts = NULL;
00976 option_data->c = 0;
00977 option_data->flofmax = -1;
00978
00979 option_data->numf = 0;
00980 option_data->numr = 0;
00981 option_data->numc = 0;
00982 option_data->numt = 0;
00983
00984 option_data->fcoef_filename = NULL;
00985 option_data->rcoef_filename = NULL;
00986 option_data->tcoef_filename = NULL;
00987
00988 option_data->numfiles = 0;
00989
00990
00991 option_data->yname
00992 = (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
00993 for (ip = 0; ip < MAX_OBSERVATIONS; ip++)
00994 option_data->yname[ip] = NULL;
00995
00996
00997
00998 option_data->fcoef_filename =
00999 (char **) malloc (sizeof(char *) * MAX_XVARS);
01000 if (option_data->fcoef_filename == NULL)
01001 RA_error ("Unable to allocate memory for fcoef_filename");
01002
01003 option_data->rcoef_filename =
01004 (char **) malloc (sizeof(char *) * MAX_XVARS);
01005 if (option_data->rcoef_filename == NULL)
01006 RA_error ("Unable to allocate memory for rcoef_filename");
01007
01008 option_data->tcoef_filename =
01009 (char **) malloc (sizeof(char *) * MAX_XVARS);
01010 if (option_data->tcoef_filename == NULL)
01011 RA_error ("Unable to allocate memory for tcoef_filename");
01012
01013 for (ip = 0; ip < MAX_XVARS; ip++)
01014 {
01015 option_data->fcoef_filename[ip] = NULL;
01016 option_data->rcoef_filename[ip] = NULL;
01017 option_data->tcoef_filename[ip] = NULL;
01018 }
01019
01020
01021
01022 option_data->bucket_filename = NULL;
01023 option_data->numbricks = -1;
01024 option_data->brick_type = NULL;
01025 option_data->brick_coef = NULL;
01026 option_data->brick_label = NULL;
01027
01028 }
01029
01030
01031
01032
01033
01034
01035
01036 void sort_model_indices
01037 (
01038 model * regmodel
01039 )
01040
01041 {
01042 int i, j, temp;
01043
01044
01045
01046 for (i = 0; i < regmodel->p - 1; i++)
01047 for (j = i+1; j < regmodel->p; j++)
01048 if (regmodel->flist[i] > regmodel->flist[j])
01049 {
01050 temp = regmodel->flist[i];
01051 regmodel->flist[i] = regmodel->flist[j];
01052 regmodel->flist[j] = temp;
01053 }
01054 else if (regmodel->flist[i] == regmodel->flist[j])
01055 RA_error ("Duplicate variable indices in model definition");
01056
01057
01058
01059 for (i = 0; i < regmodel->q - 1; i++)
01060 for (j = i+1; j < regmodel->q; j++)
01061 if (regmodel->rlist[i] > regmodel->rlist[j])
01062 {
01063 temp = regmodel->rlist[i];
01064 regmodel->rlist[i] = regmodel->rlist[j];
01065 regmodel->rlist[j] = temp;
01066 }
01067
01068 }
01069
01070
01071
01072
01073
01074
01075
01076 void get_inputs
01077 (
01078 int argc,
01079 char ** argv,
01080 matrix * xdata,
01081 model * regmodel,
01082 RA_options * option_data
01083 )
01084
01085 {
01086 const MAX_BRICKS = 100;
01087 int nopt = 1;
01088 int ival, index;
01089 float fval;
01090 int rows, cols;
01091 int irows, jcols;
01092 THD_3dim_dataset * dset=NULL;
01093 char message[MAX_NAME_LENGTH];
01094 char label[MAX_NAME_LENGTH];
01095 int ibrick;
01096 int nbricks;
01097 int p;
01098 int ip, ix;
01099
01100
01101
01102 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
01103
01104
01105
01106 AFNI_logger (PROGRAM_NAME,argc,argv);
01107
01108
01109
01110 initialize_options (regmodel, option_data);
01111
01112
01113
01114 while (nopt < argc && argv[nopt][0] == '-')
01115 {
01116
01117
01118 if( strncmp(argv[nopt],"-datum",6) == 0 ){
01119 if( ++nopt >= argc ) RA_error("need an argument after -datum!") ;
01120
01121 if( strcmp(argv[nopt],"short") == 0 ){
01122 option_data->datum = MRI_short ;
01123 } else if( strcmp(argv[nopt],"float") == 0 ){
01124 option_data->datum = MRI_float ;
01125 } else {
01126 char buf[256] ;
01127 sprintf(buf,"-datum of type '%s' is not supported in 3dRegAna.",
01128 argv[nopt] ) ;
01129 RA_error(buf) ;
01130 }
01131 nopt++ ; continue ;
01132 }
01133
01134
01135
01136 if( strncmp(argv[nopt],"-session",8) == 0 ){
01137 nopt++ ;
01138 if( nopt >= argc ) RA_error("need argument after -session!") ;
01139 strcpy(option_data->session , argv[nopt++]) ;
01140 continue ;
01141 }
01142
01143
01144
01145 if( strncmp(argv[nopt],"-diskspace",10) == 0 )
01146 {
01147 option_data->diskspace = 1;
01148 nopt++ ; continue ;
01149 }
01150
01151
01152
01153
01154 if( strncmp(argv[nopt],"-workmem",6) == 0 ){
01155 nopt++ ;
01156 if( nopt >= argc ) RA_error ("need argument after -workmem!") ;
01157 sscanf (argv[nopt], "%d", &ival);
01158 if( ival <= 0 ) RA_error ("illegal argument after -workmem!") ;
01159 option_data->workmem = ival ;
01160 nopt++ ; continue ;
01161 }
01162
01163
01164
01165 if (strncmp(argv[nopt], "-rmsmin", 7) == 0)
01166 {
01167 nopt++;
01168 if (nopt >= argc) RA_error ("need argument after -rmsmin ");
01169 sscanf (argv[nopt], "%f", &fval);
01170 if (fval < 0.0)
01171 RA_error ("illegal argument after -rmsmin ");
01172 option_data->rms_min = fval;
01173 nopt++;
01174 continue;
01175 }
01176
01177
01178
01179 if (strncmp(argv[nopt], "-flof", 6) == 0)
01180 {
01181 nopt++;
01182 if (nopt >= argc) RA_error ("need argument after -flof ");
01183 sscanf (argv[nopt], "%f", &fval);
01184 if ((fval < 0.0) || (fval > 1.0))
01185 RA_error ("illegal argument after -flof ");
01186 option_data->flofmax = fval;
01187 nopt++;
01188 continue;
01189 }
01190
01191
01192
01193 if (strncmp(argv[nopt], "-fdisp", 6) == 0)
01194 {
01195 nopt++;
01196 if (nopt >= argc) RA_error ("need argument after -fdisp ");
01197 sscanf (argv[nopt], "%f", &fval);
01198 option_data->fdisp = fval;
01199 nopt++;
01200 continue;
01201 }
01202
01203
01204
01205 if (strncmp(argv[nopt], "-rows", 5) == 0)
01206 {
01207 nopt++;
01208 if (nopt >= argc) RA_error ("need argument after -rows ");
01209 sscanf (argv[nopt], "%d", &ival);
01210 if ((ival <= 0) || (ival > MAX_OBSERVATIONS))
01211 RA_error ("illegal argument after -rows ");
01212 rows = ival;
01213 nopt++;
01214 continue;
01215 }
01216
01217
01218
01219 if (strncmp(argv[nopt], "-cols", 5) == 0)
01220 {
01221 nopt++;
01222 if (nopt >= argc) RA_error ("need argument after -cols ");
01223 sscanf (argv[nopt], "%d", &ival);
01224 if ((ival <= 0) || (ival > MAX_XVARS))
01225 RA_error ("illegal argument after -cols ");
01226 cols = ival;
01227 nopt++;
01228
01229 matrix_create (rows, cols+1, xdata);
01230 irows = -1;
01231 continue;
01232 }
01233
01234
01235
01236 if (strncmp(argv[nopt], "-xydata", 7) == 0)
01237 {
01238 nopt++;
01239 if (nopt + cols >= argc)
01240 RA_error ("need cols+1 arguments after -xydata ");
01241
01242 irows++;
01243 if (irows > rows-1) RA_error ("too many data files");
01244
01245 xdata->elts[irows][0] = 1.0;
01246 for (jcols = 1; jcols <= cols; jcols++)
01247 {
01248 sscanf (argv[nopt], "%f", &fval);
01249 xdata->elts[irows][jcols] = fval;
01250 nopt++;
01251 }
01252
01253
01254 dset = THD_open_dataset( argv[nopt] ) ;
01255 if( ! ISVALID_3DIM_DATASET(dset) )
01256 {
01257 sprintf(message,"Unable to open dataset file %s\n", argv[nopt]);
01258 RA_error (message);
01259 }
01260 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
01261
01262 option_data->yname[irows]
01263 = malloc (sizeof(char) * MAX_NAME_LENGTH);
01264 strcpy (option_data->yname[irows], argv[nopt]);
01265 nopt++;
01266 continue;
01267 }
01268
01269
01270
01271 if (strncmp(argv[nopt], "-model", 6) == 0)
01272 {
01273 nopt++;
01274 if (nopt >= argc) RA_error ("need arguments after -model ");
01275
01276 while ((nopt < argc)
01277 && (strncmp(argv[nopt], "-", 1) != 0)
01278 && (strncmp(argv[nopt], ":", 1) != 0))
01279 {
01280 sscanf (argv[nopt], "%d", &ival);
01281 if ((ival <= 0) || (ival > cols))
01282 RA_error ("illegal argument after -model ");
01283 regmodel->flist[regmodel->p] = ival;
01284
01285 regmodel->p++;
01286 nopt++;
01287 }
01288
01289 if (strncmp(argv[nopt], ":", 1) == 0)
01290 {
01291 nopt++;
01292
01293 while ((nopt < argc)
01294 && (strncmp(argv[nopt], "-", 1) != 0))
01295 {
01296 sscanf (argv[nopt], "%d", &ival);
01297 if ((ival < 0) || (ival > cols))
01298 RA_error ("illegal argument after -model ");
01299 regmodel->flist[regmodel->p] = ival;
01300 regmodel->rlist[regmodel->q] = ival;
01301 regmodel->p++;
01302 regmodel->q++;
01303 nopt++;
01304 }
01305 }
01306
01307
01308 sort_model_indices (regmodel);
01309
01310 continue;
01311 }
01312
01313
01314
01315 if (strncmp(argv[nopt], "-fcoef", 6) == 0)
01316 {
01317 nopt++;
01318 if (nopt+1 >= argc) RA_error ("need 2 arguments after -fcoef ");
01319 sscanf (argv[nopt], "%d", &ival);
01320 if ((ival < 0) || (ival > cols))
01321 RA_error ("illegal argument after -fcoef ");
01322 index = ival;
01323 nopt++;
01324
01325 option_data->fcoef_filename[index] =
01326 malloc (sizeof(char) * MAX_NAME_LENGTH);
01327 if (option_data->fcoef_filename[index] == NULL)
01328 RA_error ("Unable to allocate memory for fcoef_filename");
01329 strcpy (option_data->fcoef_filename[index], argv[nopt]);
01330
01331 nopt++;
01332 continue;
01333 }
01334
01335
01336
01337 if (strncmp(argv[nopt], "-rcoef", 6) == 0)
01338 {
01339 nopt++;
01340 if (nopt+1 >= argc) RA_error ("need 2 arguments after -rcoef ");
01341 sscanf (argv[nopt], "%d", &ival);
01342 if ((ival < 0) || (ival > cols))
01343 RA_error ("illegal argument after -rcoef ");
01344 index = ival;
01345 nopt++;
01346
01347 option_data->rcoef_filename[index] =
01348 malloc (sizeof(char) * MAX_NAME_LENGTH);
01349 if (option_data->rcoef_filename[index] == NULL)
01350 RA_error ("Unable to allocate memory for rcoef_filename");
01351 strcpy (option_data->rcoef_filename[index], argv[nopt]);
01352
01353 nopt++;
01354 continue;
01355 }
01356
01357
01358
01359 if (strncmp(argv[nopt], "-tcoef", 6) == 0)
01360 {
01361 nopt++;
01362 if (nopt+1 >= argc) RA_error ("need 2 arguments after -tcoef ");
01363 sscanf (argv[nopt], "%d", &ival);
01364 if ((ival < 0) || (ival > cols))
01365 RA_error ("illegal argument after -tcoef ");
01366 index = ival;
01367 nopt++;
01368
01369 option_data->tcoef_filename[index] =
01370 malloc (sizeof(char) * MAX_NAME_LENGTH);
01371 if (option_data->tcoef_filename[index] == NULL)
01372 RA_error ("Unable to allocate memory for tcoef_filename");
01373 strcpy (option_data->tcoef_filename[index], argv[nopt]);
01374
01375 nopt++;
01376 continue;
01377 }
01378
01379
01380
01381 if (strncmp(argv[nopt], "-bucket", 7) == 0)
01382 {
01383 nopt++;
01384 if (nopt+1 >= argc) RA_error ("need 2 arguments after -bucket ");
01385 sscanf (argv[nopt], "%d", &ival);
01386 if ((ival < 0) || (ival > MAX_BRICKS))
01387 RA_error ("illegal argument after -bucket ");
01388 option_data->numbricks = ival;
01389 nopt++;
01390
01391 option_data->bucket_filename =
01392 malloc (sizeof(char) * MAX_NAME_LENGTH);
01393 if (option_data->bucket_filename == NULL)
01394 RA_error ("Unable to allocate memory for bucket_filename");
01395 strcpy (option_data->bucket_filename, argv[nopt]);
01396
01397
01398 p = regmodel->p;
01399 if (ival == 0)
01400 nbricks = 2*p + 2;
01401 else
01402 nbricks = ival;
01403 option_data->numbricks = nbricks;
01404
01405
01406 option_data->brick_type = malloc (sizeof(int) * nbricks);
01407 option_data->brick_coef = malloc (sizeof(int) * nbricks);
01408 option_data->brick_label = malloc (sizeof(char *) * nbricks);
01409 for (ibrick = 0; ibrick < nbricks; ibrick++)
01410 {
01411 option_data->brick_type[ibrick] = -1;
01412 option_data->brick_coef[ibrick] = -1;
01413 option_data->brick_label[ibrick] =
01414 malloc (sizeof(char) * MAX_NAME_LENGTH);
01415 }
01416
01417 if (ival == 0)
01418 {
01419 for (ip = 0; ip < p; ip++)
01420 {
01421 ix = regmodel->flist[ip];
01422
01423 ibrick = 2*ip;
01424 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01425 option_data->brick_coef[ibrick] = ix;
01426 sprintf (label, "Coef #%d est.", ix);
01427 strcpy (option_data->brick_label[ibrick], label);
01428
01429 ibrick = 2*ip + 1;
01430 option_data->brick_type[ibrick] = FUNC_TT_TYPE;
01431 option_data->brick_coef[ibrick] = ix;
01432 sprintf (label, "Coef #%d t-stat", ix);
01433 strcpy (option_data->brick_label[ibrick], label);
01434 }
01435
01436 ibrick = 2*p;
01437 option_data->brick_type[ibrick] = FUNC_FT_TYPE;
01438 strcpy (option_data->brick_label[ibrick], "F-stat Regression");
01439
01440 ibrick = 2*p + 1;
01441 option_data->brick_type[ibrick] = FUNC_THR_TYPE;
01442 strcpy (option_data->brick_label[ibrick], "R^2 Regression");
01443
01444 }
01445
01446 nopt++;
01447 continue;
01448 }
01449
01450
01451
01452 if (strncmp(argv[nopt], "-brick", 6) == 0)
01453 {
01454 nopt++;
01455 if (nopt+2 >= argc) RA_error ("need more arguments after -brick ");
01456 sscanf (argv[nopt], "%d", &ibrick);
01457 if ((ibrick < 0) || (ibrick >= option_data->numbricks))
01458 RA_error ("illegal argument after -brick ");
01459 nopt++;
01460
01461 if (strncmp(argv[nopt], "coef", 4) == 0)
01462 {
01463 option_data->brick_type[ibrick] = FUNC_FIM_TYPE;
01464
01465 nopt++;
01466 sscanf (argv[nopt], "%d", &ival);
01467 if ((ival < 0) || (ival > cols))
01468 RA_error ("illegal argument after coef ");
01469 option_data->brick_coef[ibrick] = ival;
01470
01471 nopt++;
01472 if (nopt >= argc) RA_error ("need more arguments after -brick ");
01473 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01474
01475 }
01476 else if (strncmp(argv[nopt], "tstat", 4) == 0)
01477 {
01478 option_data->brick_type[ibrick] = FUNC_TT_TYPE;
01479
01480 nopt++;
01481 sscanf (argv[nopt], "%d", &ival);
01482 if ((ival < 0) || (ival > cols))
01483 RA_error ("illegal argument after tstat ");
01484 option_data->brick_coef[ibrick] = ival;
01485
01486 nopt++;
01487 if (nopt >= argc) RA_error ("need more arguments after -brick ");
01488 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01489
01490 }
01491 else if (strncmp(argv[nopt], "fstat", 4) == 0)
01492 {
01493 option_data->brick_type[ibrick] = FUNC_FT_TYPE;
01494
01495 nopt++;
01496 if (nopt >= argc) RA_error ("need more arguments after -brick ");
01497 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01498
01499 }
01500 else if (strncmp(argv[nopt], "rstat", 4) == 0)
01501 {
01502 option_data->brick_type[ibrick] = FUNC_THR_TYPE;
01503
01504 nopt++;
01505 if (nopt >= argc) RA_error ("need more arguments after -brick ");
01506 strcpy (option_data->brick_label[ibrick], argv[nopt]);
01507
01508 }
01509 else
01510 {
01511 sprintf(message,"Unrecognized option after -brick: %s\n",
01512 argv[nopt]);
01513 RA_error (message);
01514 }
01515
01516 nopt++;
01517 continue;
01518 }
01519
01520
01521
01522 sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
01523 RA_error (message);
01524
01525 }
01526
01527
01528 if (irows < rows-1)
01529 RA_error ("Fewer than expected datasets were entered");
01530 }
01531
01532
01533
01534
01535
01536
01537
01538 void count_volumes_and_files
01539 (
01540 model * regmodel,
01541 RA_options * option_data
01542 )
01543
01544 {
01545 int ip;
01546 int ix;
01547 int ib;
01548 int p;
01549 int nbricks;
01550
01551
01552
01553 nbricks = option_data->numbricks;
01554 p = regmodel->p;
01555
01556
01557
01558 for (ip = 0; ip < p; ip++)
01559 {
01560 ix = regmodel->flist[ip];
01561
01562 if (option_data->fcoef_filename[ix] != NULL)
01563 {
01564 option_data->numf = 1;
01565 option_data->numfiles += 1;
01566 }
01567 }
01568
01569 for (ib = 0; ib < nbricks; ib++)
01570 if (option_data->brick_type[ib] == FUNC_FT_TYPE)
01571 option_data->numf = 1;
01572
01573
01574
01575 for (ip = 0; ip < p; ip++)
01576 {
01577 ix = regmodel->flist[ip];
01578
01579 if (option_data->rcoef_filename[ix] != NULL)
01580 {
01581 option_data->numr = 1;
01582 option_data->numfiles += 1;
01583 }
01584 }
01585
01586 for (ib = 0; ib < nbricks; ib++)
01587 if (option_data->brick_type[ib] == FUNC_THR_TYPE)
01588 option_data->numr = 1;
01589
01590
01591
01592 for (ip = 0; ip < p; ip++)
01593 {
01594 ix = regmodel->flist[ip];
01595
01596 if (option_data->tcoef_filename[ix] != NULL)
01597 {
01598 option_data->numt += 1;
01599 option_data->numfiles += 1;
01600 }
01601
01602 else
01603 for (ib = 0; ib < nbricks; ib++)
01604 if ((option_data->brick_type[ib] == FUNC_TT_TYPE)
01605 && (option_data->brick_coef[ib] == ix))
01606 option_data->numt += 1;
01607 }
01608
01609
01610
01611 for (ip = 0; ip < p; ip++)
01612 {
01613 ix = regmodel->flist[ip];
01614
01615 if ( (option_data->fcoef_filename[ix] != NULL)
01616 || (option_data->rcoef_filename[ix] != NULL)
01617 || (option_data->tcoef_filename[ix] != NULL))
01618 option_data->numc += 1;
01619
01620 else
01621 for (ib = 0; ib < nbricks; ib++)
01622 if ((option_data->brick_type[ib] == FUNC_FIM_TYPE)
01623 && (option_data->brick_coef[ib] == ix))
01624 option_data->numc += 1;
01625 }
01626
01627 }
01628
01629
01630
01631
01632
01633
01634
01635 void break_into_pieces
01636 (
01637 int num_datasets,
01638 RA_options * option_data
01639 )
01640
01641 {
01642 int num_vols;
01643
01644
01645
01646 num_vols = option_data->numf + option_data->numr
01647 + option_data->numt + option_data->numc;
01648
01649
01650 option_data->piece_size = option_data->workmem * MEGA
01651 / ((num_datasets + num_vols) * sizeof(float));
01652 if (option_data->piece_size > option_data->nxyz)
01653 option_data->piece_size = option_data->nxyz;
01654
01655
01656 option_data->num_pieces = (option_data->nxyz + option_data->piece_size - 1)
01657 / option_data->piece_size;
01658
01659 printf ("num_pieces = %d piece_size = %d \n",
01660 option_data->num_pieces, option_data->piece_size);
01661
01662 }
01663
01664
01665
01666
01667
01668
01669
01670
01671
01672 void identify_repeats
01673 (
01674 matrix * xdata,
01675 model * regmodel,
01676 RA_options * option_data
01677 )
01678
01679 {
01680 int i, k;
01681 int j;
01682 int match;
01683
01684 int which;
01685 double p, q;
01686 double f;
01687 double dfn, dfd;
01688 int status;
01689 double bound;
01690
01691
01692
01693 option_data->levels = (int *) malloc (sizeof(int) * (xdata->rows));
01694 MTEST (option_data->levels);
01695 option_data->counts = (int *) malloc (sizeof(int) * (xdata->rows));
01696 MTEST (option_data->counts);
01697
01698
01699
01700 for (i = 0; i < xdata->rows; i++)
01701 option_data->counts[i] = 0;
01702 option_data->levels[0] = 0;
01703 option_data->counts[0] = 1;
01704 option_data->c = 1;
01705
01706
01707
01708 for (i = 1; i < xdata->rows; i++)
01709 {
01710 option_data->levels[i] = option_data->c;
01711
01712
01713 for (k = 0; k < i; k++)
01714 {
01715 match = 1;
01716 for (j = 1; j < xdata->cols; j++)
01717 if (xdata->elts[i][j] != xdata->elts[k][j]) match = 0;
01718
01719 if (match)
01720 {
01721 option_data->levels[i] = option_data->levels[k];
01722 break;
01723 }
01724 }
01725
01726
01727 k = option_data->levels[i];
01728 option_data->counts[k] ++;
01729
01730
01731 if (k == option_data->c) option_data->c++;
01732 }
01733
01734
01735
01736 which = 2;
01737 q = option_data->flofmax;
01738 p = 1.0 - q;
01739 dfn = option_data->c - regmodel->p;
01740 dfd = xdata->rows - option_data->c;
01741 cdff (&which, &p, &q, &f, &dfn, &dfd, &status, &bound);
01742 if (status != 0) RA_error ("Error in calculating F - lack of fit ");
01743 option_data->flofmax = f;
01744
01745
01746 }
01747
01748
01749
01750
01751
01752
01753
01754 void check_for_valid_inputs
01755 (
01756 matrix * xdata,
01757 model * regmodel,
01758 RA_options * option_data
01759 )
01760
01761 {
01762 int ib;
01763
01764
01765
01766 if (xdata->cols < 2)
01767 RA_error ("Too few X variables ");
01768 if (xdata->rows < 3)
01769 RA_error ("Too few data sets for Y-observations ");
01770
01771
01772 if (regmodel->q < 1)
01773 RA_error ("Reduced regression model is too small");
01774 if (regmodel->p <= regmodel->q)
01775 RA_error ("Full regression model is too small");
01776 if (xdata->rows <= regmodel->p)
01777 RA_error ("Too few data sets for fitting full regression model ");
01778 if (regmodel->rlist[0] != 0)
01779 RA_error ("Reduced model must include variable index 0");
01780
01781
01782
01783 if (option_data->flofmax >= 0.0)
01784 {
01785 if (xdata->rows <= option_data->c)
01786 RA_error ("Cannot conduct F-test for lack of fit (no repeat obs.) ");
01787 if (option_data->c <= regmodel->p)
01788 RA_error ("Cannot conduct F-test for lack of fit (c <= p) ");
01789 }
01790
01791
01792 if (option_data->numbricks > 0)
01793 for (ib = 0; ib < option_data->numbricks; ib++)
01794 {
01795 if (option_data->brick_type[ib] < 0)
01796 RA_error ("Must specify contents of each brick in the bucket");
01797 }
01798 }
01799
01800
01801
01802
01803
01804
01805
01806 void check_one_output_file
01807 (
01808 RA_options * option_data,
01809 char * filename
01810 )
01811
01812 {
01813 THD_3dim_dataset * dset=NULL;
01814 THD_3dim_dataset * new_dset=NULL;
01815 int ierror;
01816
01817
01818
01819 dset = THD_open_dataset (option_data->first_dataset ) ;
01820 if( ! ISVALID_3DIM_DATASET(dset) ){
01821 fprintf(stderr,"*** Unable to open dataset file %s\n",
01822 option_data->first_dataset);
01823 exit(1) ;
01824 }
01825
01826
01827 new_dset = EDIT_empty_copy( dset ) ;
01828
01829
01830 ierror = EDIT_dset_items( new_dset ,
01831 ADN_prefix , filename ,
01832 ADN_label1 , filename ,
01833 ADN_directory_name , option_data->session ,
01834 ADN_self_name , filename ,
01835 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE :
01836 GEN_FUNC_TYPE ,
01837 ADN_none ) ;
01838
01839 if( ierror > 0 ){
01840 fprintf(stderr,
01841 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
01842 exit(1) ;
01843 }
01844
01845 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
01846 fprintf(stderr,
01847 "*** Output dataset file %s already exists--cannot continue!\a\n",
01848 new_dset->dblk->diskptr->header_name ) ;
01849 exit(1) ;
01850 }
01851
01852
01853 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
01854 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01855
01856 }
01857
01858
01859
01860
01861
01862
01863
01864 void check_output_files
01865 (
01866 RA_options * option_data
01867 )
01868
01869 {
01870 int ix;
01871
01872
01873 for (ix = 0; ix < MAX_XVARS; ix++)
01874 {
01875 if (option_data->fcoef_filename[ix] != NULL)
01876 check_one_output_file (option_data, option_data->fcoef_filename[ix]);
01877 if (option_data->rcoef_filename[ix] != NULL)
01878 check_one_output_file (option_data, option_data->rcoef_filename[ix]);
01879 if (option_data->tcoef_filename[ix] != NULL)
01880 check_one_output_file (option_data, option_data->tcoef_filename[ix]);
01881 }
01882
01883 if (option_data->bucket_filename != NULL)
01884 check_one_output_file (option_data, option_data->bucket_filename);
01885
01886 }
01887
01888
01889
01890
01891
01892
01893
01894 void check_temporary_files
01895 (
01896 matrix * xdata,
01897 model * regmodel,
01898 RA_options * option_data
01899 )
01900
01901 {
01902 int p;
01903 int ip, ix;
01904 int num_pieces;
01905 char filename[MAX_NAME_LENGTH];
01906
01907
01908
01909 p = regmodel->p;
01910 num_pieces = option_data->num_pieces;
01911
01912
01913
01914 if (option_data->numf > 0)
01915 {
01916 strcpy (filename, "freg");
01917 check_volume (filename, num_pieces);
01918 }
01919
01920
01921
01922 if (option_data->numr > 0)
01923 {
01924 strcpy (filename, "rsqr");
01925 check_volume (filename, num_pieces);
01926 }
01927
01928
01929
01930 if (option_data->numt > 0)
01931 {
01932 for (ip = 0; ip < p; ip++)
01933 {
01934 ix = regmodel->flist[ip];
01935
01936 if (option_data->tcoef_filename[ix] != NULL)
01937 {
01938 sprintf (filename, "tcoef.%d", ix);
01939 check_volume (filename, num_pieces);
01940 }
01941 }
01942 }
01943
01944
01945
01946 if (option_data->numc > 0)
01947 {
01948 for (ip = 0; ip < p; ip++)
01949 {
01950 ix = regmodel->flist[ip];
01951
01952 if ( (option_data->fcoef_filename[ix] != NULL)
01953 || (option_data->rcoef_filename[ix] != NULL)
01954 || (option_data->tcoef_filename[ix] != NULL) )
01955 {
01956 sprintf (filename, "coef.%d", ix);
01957 check_volume (filename, num_pieces);
01958 }
01959 }
01960 }
01961
01962
01963 }
01964
01965
01966
01967
01968
01969
01970
01971
01972 void check_disk_space
01973 (
01974 RA_options * option_data
01975 )
01976
01977 {
01978 char ch;
01979 int nxyz;
01980 int nmax;
01981 char filename[MAX_NAME_LENGTH];
01982
01983
01984
01985 nxyz = option_data->nxyz;
01986
01987
01988 nmax = 4 * nxyz * (option_data->numf + option_data->numr + option_data->numt
01989 + option_data->numc);
01990
01991
01992 nmax += 4 * nxyz * option_data->numfiles;
01993
01994 printf("\nThis problem requires approximately %d MB of free disk space.\n",
01995 (nmax/1000000) + 1);
01996
01997
01998
01999 printf ("Summary of available disk space: \n\n");
02000 system (DF);
02001 printf ("\nDo you wish to proceed? ");
02002 ch = getchar();
02003 if ((ch != 'Y') && (ch != 'y')) exit(0);
02004 printf ("\n");
02005 }
02006
02007
02008
02009
02010
02011
02012
02013 void initialize_program
02014 (
02015 int argc,
02016 char ** argv,
02017 matrix * xdata,
02018 model * regmodel,
02019 RA_options * option_data
02020 )
02021
02022 {
02023
02024
02025 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
02026
02027
02028
02029 matrix_initialize (xdata);
02030
02031
02032
02033 get_inputs (argc, argv, xdata, regmodel, option_data);
02034
02035
02036
02037 option_data->first_dataset = option_data->yname[0];
02038 get_dimensions (option_data);
02039 printf ("Data set dimensions: nx = %d ny = %d nz = %d nxyz = %d \n",
02040 option_data->nx, option_data->ny, option_data->nz, option_data->nxyz);
02041
02042
02043
02044 count_volumes_and_files (regmodel, option_data);
02045
02046
02047
02048 break_into_pieces (xdata->rows, option_data);
02049
02050
02051
02052 if (option_data->flofmax >= 0.0)
02053 identify_repeats (xdata, regmodel, option_data);
02054
02055
02056
02057 check_for_valid_inputs (xdata, regmodel, option_data);
02058
02059
02060
02061 check_output_files (option_data);
02062
02063
02064
02065 check_temporary_files (xdata, regmodel, option_data);
02066
02067
02068
02069 if (option_data->diskspace) check_disk_space (option_data);
02070
02071 }
02072
02073
02074
02075
02076
02077
02078
02079 void init_regression_analysis
02080 (
02081 int p,
02082 int q,
02083 int * flist,
02084 int * rlist,
02085 matrix xdata,
02086 matrix * x_full,
02087 matrix * xtxinv_full,
02088 matrix * xtxinvxt_full,
02089 matrix * x_base,
02090 matrix * xtxinvxt_base,
02091 matrix * x_rdcd,
02092 matrix * xtxinvxt_rdcd
02093 )
02094
02095 {
02096 int zlist[MAX_XVARS];
02097 int ip;
02098 matrix xtxinv_temp;
02099
02100
02101
02102 matrix_initialize (&xtxinv_temp);
02103
02104
02105
02106 for (ip = 0; ip < MAX_XVARS; ip++)
02107 zlist[ip] = 0;
02108
02109
02110
02111 calc_matrices (xdata, 1, zlist, x_base, &xtxinv_temp, xtxinvxt_base);
02112 calc_matrices (xdata, q, rlist, x_rdcd, &xtxinv_temp, xtxinvxt_rdcd);
02113 calc_matrices (xdata, p, flist, x_full, xtxinv_full, xtxinvxt_full);
02114
02115
02116
02117 matrix_destroy (&xtxinv_temp);
02118
02119 }
02120
02121
02122
02123
02124
02125
02126
02127 void regression_analysis
02128 (
02129 int N,
02130 int p,
02131 int q,
02132 matrix x_full,
02133 matrix xtxinv_full,
02134 matrix xtxinvxt_full,
02135 matrix x_base,
02136 matrix xtxinvxt_base,
02137 matrix x_rdcd,
02138 matrix xtxinvxt_rdcd,
02139 vector y,
02140 float rms_min,
02141 int * levels,
02142 int * counts,
02143 int c,
02144 float flofmax,
02145 float * flof,
02146 vector * coef_full,
02147 vector * scoef_full,
02148 vector * tcoef_full,
02149 float * freg,
02150 float * rsqr
02151 )
02152
02153 {
02154 float sse_base;
02155 float sse_rdcd;
02156 float sse_full;
02157 float sspe;
02158 vector coef_temp;
02159
02160
02161
02162 vector_initialize (&coef_temp);
02163
02164
02165
02166 calc_coef (xtxinvxt_base, y, &coef_temp);
02167
02168
02169
02170 sse_base = calc_sse (x_base, coef_temp, y);
02171
02172
02173
02174 if (sqrt(sse_base/N) < rms_min)
02175 {
02176 vector_create (p, coef_full);
02177 vector_create (p, scoef_full);
02178 vector_create (p, tcoef_full);
02179 *freg = 0.0;
02180 *rsqr = 0.0;
02181 vector_destroy (&coef_temp);
02182 return;
02183 }
02184
02185
02186
02187 calc_coef (xtxinvxt_full, y, coef_full);
02188
02189
02190
02191 sse_full = calc_sse (x_full, *coef_full, y);
02192
02193
02194
02195 calc_tcoef (N, p, sse_full, xtxinv_full,
02196 *coef_full, scoef_full, tcoef_full);
02197
02198
02199
02200 if (flofmax > 0.0)
02201 {
02202
02203 sspe = calc_sspe (y, levels, counts, c);
02204
02205
02206 *flof = calc_flof (N, p, c, sse_full, sspe);
02207
02208 if (*flof > flofmax)
02209 {
02210 vector_create (p, coef_full);
02211 vector_create (p, scoef_full);
02212 vector_create (p, tcoef_full);
02213 *freg = 0.0;
02214 *rsqr = 0.0;
02215 vector_destroy (&coef_temp);
02216 return;
02217 }
02218 }
02219 else
02220 *flof = -1.0;
02221
02222
02223
02224 calc_coef (xtxinvxt_rdcd, y, &coef_temp);
02225
02226
02227
02228 sse_rdcd = calc_sse (x_rdcd, coef_temp, y);
02229
02230
02231
02232 *freg = calc_freg (N, p, q, sse_full, sse_rdcd);
02233
02234
02235
02236 *rsqr = calc_rsqr (sse_full, sse_base);
02237
02238
02239
02240 vector_destroy (&coef_temp);
02241
02242 }
02243
02244
02245
02246
02247
02248
02249
02250 void save_voxel
02251 (
02252 int iv,
02253 vector y,
02254 float fdisp,
02255 model * regmodel,
02256 float flof,
02257 vector coef,
02258 vector tcoef,
02259 float freg,
02260 float rsqr,
02261
02262 float * freg_piece,
02263 float * rsqr_piece,
02264 float ** coef_piece,
02265 float ** tcoef_piece
02266 )
02267
02268 {
02269 int ip;
02270 int ix;
02271
02272
02273
02274 if (freg_piece != NULL) freg_piece[iv] = freg;
02275 if (rsqr_piece != NULL) rsqr_piece[iv] = rsqr;
02276
02277
02278
02279 for (ip = 0; ip < regmodel->p; ip++)
02280 {
02281 ix = regmodel->flist[ip];
02282
02283 if (coef_piece[ix] != NULL) coef_piece[ix][iv] = coef.elts[ip];
02284
02285 if (tcoef_piece[ix] != NULL) tcoef_piece[ix][iv] = tcoef.elts[ip];
02286
02287 }
02288
02289
02290
02291 if ((fdisp >= 0.0) && (freg >= fdisp))
02292 {
02293 printf ("\n\nVoxel #%d: \n", iv);
02294 printf ("\nY data: \n");
02295 for (ip = 0; ip < y.dim; ip++)
02296 printf ("Y[%d] = %f \n", ip, y.elts[ip]);
02297
02298 if (flof >= 0.0) printf ("\nF lack of fit = %f \n", flof);
02299 printf ("\nF regression = %f \n", freg);
02300 printf ("R-squared = %f \n", rsqr);
02301
02302 printf ("\nRegression Coefficients: \n");
02303 for (ip = 0; ip < coef.dim; ip++)
02304 {
02305 ix = regmodel->flist[ip];
02306 printf ("b[%d] = %f ", ix, coef.elts[ip]);
02307 printf ("t[%d] = %f \n", ix, tcoef.elts[ip]);
02308 }
02309
02310 }
02311
02312 }
02313
02314
02315
02316
02317
02318
02319
02320 void calculate_results
02321 (
02322 matrix xdata,
02323 model * regmodel,
02324 RA_options * option_data
02325 )
02326
02327 {
02328 int * flist = NULL;
02329 int * rlist = NULL;
02330 int n;
02331 int p;
02332 int q;
02333 float flof;
02334 float freg;
02335 float rsqr;
02336 vector coef;
02337 vector scoef;
02338 vector tcoef;
02339
02340 matrix x_full;
02341 matrix xtxinv_full;
02342 matrix xtxinvxt_full;
02343 matrix x_base;
02344 matrix xtxinvxt_base;
02345 matrix x_rdcd;
02346 matrix xtxinvxt_rdcd;
02347 vector y;
02348
02349 int i;
02350 int nxyz;
02351 int num_datasets;
02352 int piece_size;
02353 int num_pieces;
02354 int piece;
02355 int piece_len;
02356 int fim_offset;
02357 int ivox;
02358 int nvox;
02359
02360 float ** yfimar = NULL;
02361 float * freg_piece = NULL;
02362 float * rsqr_piece = NULL;
02363 float ** coef_piece = NULL;
02364 float ** tcoef_piece = NULL;
02365
02366
02367
02368 matrix_initialize (&x_full);
02369 matrix_initialize (&xtxinv_full);
02370 matrix_initialize (&xtxinvxt_full);
02371 matrix_initialize (&x_base);
02372 matrix_initialize (&xtxinvxt_base);
02373 matrix_initialize (&x_rdcd);
02374 matrix_initialize (&xtxinvxt_rdcd);
02375 vector_initialize (&coef);
02376 vector_initialize (&scoef);
02377 vector_initialize (&tcoef);
02378 vector_initialize (&y);
02379
02380
02381
02382 n = xdata.rows;
02383 p = regmodel->p;
02384 flist = regmodel->flist;
02385 q = regmodel->q;
02386 rlist = regmodel->rlist;
02387 nxyz = option_data->nxyz;
02388 piece_size = option_data->piece_size;
02389 num_pieces = option_data->num_pieces;
02390
02391
02392
02393 allocate_pieces (xdata, regmodel, option_data, &yfimar,
02394 &freg_piece, &rsqr_piece, &coef_piece, &tcoef_piece);
02395
02396
02397
02398 init_regression_analysis (p, q, flist, rlist, xdata,
02399 &x_full, &xtxinv_full, &xtxinvxt_full,
02400 &x_base, &xtxinvxt_base, &x_rdcd, &xtxinvxt_rdcd);
02401
02402
02403 vector_create (n, &y);
02404
02405
02406 if (option_data->fdisp >= 0)
02407 {
02408 printf ("\n");
02409 printf ("X matrix: \n");
02410 matrix_print (xdata);
02411 }
02412
02413
02414
02415 nvox = -1;
02416 for (piece = 0; piece < num_pieces; piece++)
02417 {
02418 printf ("piece = %d \n", piece);
02419 fim_offset = piece * piece_size;
02420 if (piece < num_pieces-1)
02421 piece_len = piece_size;
02422 else
02423 piece_len = nxyz - fim_offset;
02424
02425
02426 for (i = 0; i < n; i++)
02427 read_afni_data (option_data, option_data->yname[i],
02428 piece_len, fim_offset, yfimar[i]);
02429
02430
02431
02432 for (ivox = 0; ivox < piece_len; ivox++)
02433 {
02434 nvox += 1;
02435
02436
02437
02438 for (i = 0; i < n; i++)
02439 y.elts[i] = yfimar[i][ivox];
02440
02441
02442
02443 regression_analysis (n, p, q,
02444 x_full, xtxinv_full, xtxinvxt_full, x_base,
02445 xtxinvxt_base, x_rdcd, xtxinvxt_rdcd,
02446 y, option_data->rms_min, option_data->levels,
02447 option_data->counts, option_data->c,
02448 option_data->flofmax, &flof,
02449 &coef, &scoef, &tcoef, &freg, &rsqr);
02450
02451
02452
02453 save_voxel (ivox, y, option_data->fdisp,
02454 regmodel, flof, coef, tcoef, freg, rsqr,
02455 freg_piece, rsqr_piece, coef_piece, tcoef_piece);
02456
02457
02458 }
02459
02460
02461
02462 save_pieces (piece, piece_len,
02463 freg_piece, rsqr_piece, coef_piece, tcoef_piece);
02464
02465
02466 }
02467
02468
02469
02470 deallocate_pieces (n, &yfimar, &freg_piece, &rsqr_piece,
02471 &coef_piece, &tcoef_piece);
02472
02473
02474
02475 vector_destroy (&y);
02476 vector_destroy (&tcoef);
02477 vector_destroy (&scoef);
02478 vector_destroy (&coef);
02479 matrix_destroy (&xtxinvxt_rdcd);
02480 matrix_destroy (&x_rdcd);
02481 matrix_destroy (&xtxinvxt_base);
02482 matrix_destroy (&x_base);
02483 matrix_destroy (&xtxinvxt_full);
02484 matrix_destroy (&xtxinv_full);
02485 matrix_destroy (&x_full);
02486
02487 }
02488
02489
02490
02491
02492
02493
02494
02495
02496 float EDIT_coerce_autoscale_new
02497 (
02498 int nxyz,
02499 int itype,
02500 void *ivol,
02501 int otype,
02502 void *ovol
02503 )
02504
02505 {
02506 float fac=0.0 , top ;
02507
02508 if( MRI_IS_INT_TYPE(otype) ){
02509 top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
02510 if (top == 0.0) fac = 0.0;
02511 else fac = MRI_TYPE_maxval[otype]/top;
02512 }
02513
02514 EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
02515 return ( fac );
02516 }
02517
02518
02519
02520
02521
02522
02523
02524
02525
02526
02527
02528
02529
02530
02531 void write_afni_data
02532 (
02533 RA_options * option_data,
02534 char * filename,
02535 float * ffim,
02536 float * ftr,
02537 int func_type,
02538 int numdof,
02539 int dendof
02540 )
02541
02542 {
02543 int nxyz;
02544 int ii;
02545 THD_3dim_dataset * dset=NULL;
02546 THD_3dim_dataset * new_dset=NULL;
02547 int ierror;
02548 int ibuf[32];
02549 float fbuf[MAX_STAT_AUX];
02550 float fimfac;
02551 int output_datum;
02552 short * tsp = NULL;
02553 void * vdif = NULL;
02554 float top, bot, func_scale_short;
02555 int top_ss, bot_ss;
02556 char label[80];
02557
02558
02559
02560 nxyz = option_data->nxyz;
02561
02562
02563 dset = THD_open_dataset (option_data->first_dataset) ;
02564 if( ! ISVALID_3DIM_DATASET(dset) ){
02565 fprintf(stderr,"*** Unable to open dataset file %s\n",
02566 option_data->first_dataset); exit(1) ;
02567 }
02568
02569
02570
02571 new_dset = EDIT_empty_copy( dset ) ;
02572
02573
02574
02575
02576 sprintf (label, "Output prefix: %s", filename);
02577 if( commandline != NULL )
02578 tross_multi_Append_History( new_dset , commandline,label,NULL ) ;
02579 else
02580 tross_Append_History ( new_dset, label);
02581
02582
02583 output_datum = MRI_short ;
02584 ibuf[0] = output_datum ; ibuf[1] = MRI_short ;
02585
02586
02587 ierror = EDIT_dset_items( new_dset ,
02588 ADN_prefix , filename ,
02589 ADN_label1 , filename ,
02590 ADN_directory_name , option_data->session ,
02591 ADN_self_name , filename ,
02592 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE :
02593 GEN_FUNC_TYPE ,
02594 ADN_func_type , func_type ,
02595 ADN_nvals , FUNC_nvals[func_type] ,
02596 ADN_datum_array , ibuf ,
02597 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
02598 ADN_none ) ;
02599
02600 if( ierror > 0 ){
02601 fprintf(stderr,
02602 "*** %d errors in attempting to create output dataset!\n", ierror ) ;
02603 exit(1) ;
02604 }
02605
02606 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
02607 fprintf(stderr,
02608 "*** Output dataset file %s already exists--cannot continue!\a\n",
02609 new_dset->dblk->diskptr->header_name ) ;
02610 exit(1) ;
02611 }
02612
02613
02614 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
02615
02616
02617
02618 vdif = (void *) malloc( mri_datum_size(output_datum) * nxyz ) ;
02619 tsp = (short *) malloc( sizeof(short) * nxyz ) ;
02620
02621
02622 mri_fix_data_pointer (vdif, DSET_BRICK(new_dset,0));
02623 mri_fix_data_pointer (tsp, DSET_BRICK(new_dset,1));
02624
02625
02626
02627 fimfac = EDIT_coerce_autoscale_new (nxyz,
02628 MRI_float, ffim,
02629 output_datum, vdif);
02630 if (fimfac != 0.0) fimfac = 1.0 / fimfac;
02631
02632
02633 top_ss = 32700;
02634
02635 if (func_type == FUNC_THR_TYPE)
02636 {
02637 func_scale_short = FUNC_THR_SCALE_SHORT;
02638 bot_ss = 0;
02639 }
02640 else if (func_type == FUNC_TT_TYPE)
02641 {
02642 func_scale_short = FUNC_TT_SCALE_SHORT;
02643 bot_ss = -top_ss;
02644 }
02645 else if (func_type == FUNC_FT_TYPE)
02646 {
02647 func_scale_short = FUNC_FT_SCALE_SHORT;
02648 bot_ss = 0;
02649 }
02650 else
02651 RA_error ("Illegal ouput dataset function type");
02652
02653 top = top_ss / func_scale_short;
02654 bot = bot_ss / func_scale_short;
02655
02656
02657 for (ii = 0; ii < nxyz; ii++)
02658 {
02659 if (ftr[ii] > top)
02660 tsp[ii] = top_ss;
02661 else if (ftr[ii] < bot)
02662 tsp[ii] = bot_ss;
02663 else
02664 tsp[ii] = (short) (func_scale_short * ftr[ii] + 0.5);
02665 }
02666
02667
02668
02669 if (func_type == FUNC_THR_TYPE)
02670 printf("----- Writing `fith' dataset ");
02671 else if (func_type == FUNC_TT_TYPE)
02672 printf("----- Writing `fitt' dataset ");
02673 else if (func_type == FUNC_FT_TYPE)
02674 printf("----- Writing `fift' dataset ");
02675
02676 printf("into %s\n", DSET_BRIKNAME(new_dset) ) ;
02677
02678 fbuf[0] = numdof;
02679 fbuf[1] = dendof;
02680 for( ii=2 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
02681 (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
02682
02683 fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
02684 fbuf[1] = 1.0 / func_scale_short ;
02685 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
02686
02687 THD_load_statistics( new_dset ) ;
02688 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
02689
02690
02691
02692 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
02693
02694 }
02695
02696
02697
02698
02699
02700
02701
02702 void write_bucket_data
02703 (
02704 matrix xdata,
02705 model * regmodel,
02706 RA_options * option_data
02707 )
02708
02709 {
02710 const float EPSILON = 1.0e-10;
02711 int p;
02712 int q;
02713 int n;
02714 int nxyz;
02715 THD_3dim_dataset * old_dset = NULL;
02716 THD_3dim_dataset * new_dset = NULL;
02717 char * output_prefix;
02718 char * output_session;
02719 int nbricks, ib;
02720 short ** bar = NULL;
02721 float factor;
02722 int brick_type;
02723 int brick_coef;
02724 char * brick_label;
02725 int ierror;
02726 char filename[MAX_NAME_LENGTH];
02727 int piece_size;
02728 int num_pieces;
02729 float * volume = NULL;
02730 char label[80];
02731
02732
02733
02734 p = regmodel->p;
02735 q = regmodel->q;
02736 n = xdata.rows;
02737 nxyz = option_data->nxyz;
02738 piece_size = option_data->piece_size;
02739 num_pieces = option_data->num_pieces;
02740 nbricks = option_data->numbricks;
02741 output_prefix = option_data->bucket_filename;
02742 output_session = option_data->session;
02743
02744
02745
02746 volume = (float *) malloc (sizeof(float) * nxyz);
02747 MTEST (volume);
02748 bar = (short **) malloc (sizeof(short *) * nbricks);
02749 MTEST (bar);
02750
02751
02752
02753 old_dset = THD_open_dataset (option_data->first_dataset) ;
02754
02755
02756
02757 new_dset = EDIT_empty_copy (old_dset);
02758
02759
02760
02761 if( commandline != NULL )
02762 tross_Append_History( new_dset , commandline ) ;
02763 sprintf (label, "Output prefix: %s", output_prefix);
02764 tross_Append_History ( new_dset, label);
02765
02766
02767
02768
02769 ierror = EDIT_dset_items (new_dset,
02770 ADN_prefix, output_prefix,
02771 ADN_directory_name, output_session,
02772 ADN_type, HEAD_FUNC_TYPE,
02773 ADN_func_type, FUNC_BUCK_TYPE,
02774 ADN_ntt, 0,
02775 ADN_nvals, nbricks,
02776 ADN_malloc_type, DATABLOCK_MEM_MALLOC ,
02777 ADN_none ) ;
02778
02779 if( ierror > 0 )
02780 {
02781 fprintf(stderr,
02782 "*** %d errors in attempting to create output dataset!\n",
02783 ierror);
02784 exit(1);
02785 }
02786
02787 if (THD_is_file(DSET_HEADNAME(new_dset)))
02788 {
02789 fprintf(stderr,
02790 "*** Output dataset file %s already exists--cannot continue!\n",
02791 DSET_HEADNAME(new_dset));
02792 exit(1);
02793 }
02794
02795
02796
02797 THD_delete_3dim_dataset( old_dset , False ); old_dset = NULL ;
02798
02799
02800
02801
02802
02803 for (ib = 0; ib < nbricks; ib++)
02804 {
02805
02806 brick_type = option_data->brick_type[ib];
02807 brick_coef = option_data->brick_coef[ib];
02808 brick_label = option_data->brick_label[ib];
02809
02810 if (brick_type == FUNC_FIM_TYPE)
02811 {
02812 sprintf (filename, "coef.%d", brick_coef);
02813 }
02814 else if (brick_type == FUNC_THR_TYPE)
02815 {
02816 sprintf (filename, "rsqr");
02817 }
02818 else if (brick_type == FUNC_TT_TYPE)
02819 {
02820 sprintf (filename, "tcoef.%d", brick_coef);
02821 EDIT_BRICK_TO_FITT (new_dset, ib, n-p);
02822 }
02823 else if (brick_type == FUNC_FT_TYPE)
02824 {
02825 sprintf (filename, "freg");
02826 EDIT_BRICK_TO_FIFT (new_dset, ib, p-q, n-p);
02827 }
02828
02829
02830 bar[ib] = (short *) malloc (sizeof(short) * nxyz);
02831 MTEST (bar[ib]);
02832
02833 read_volume (filename, volume, nxyz, piece_size, num_pieces);
02834 delete_volume (filename, nxyz, piece_size, num_pieces);
02835
02836 factor = EDIT_coerce_autoscale_new (nxyz, MRI_float, volume,
02837 MRI_short, bar[ib]);
02838
02839 if (factor < EPSILON) factor = 0.0;
02840 else factor = 1.0 / factor;
02841
02842
02843 EDIT_BRICK_LABEL (new_dset, ib, brick_label);
02844 EDIT_BRICK_FACTOR (new_dset, ib, factor);
02845
02846
02847
02848 EDIT_substitute_brick (new_dset, ib, MRI_short, bar[ib]);
02849
02850 }
02851
02852
02853
02854 THD_load_statistics (new_dset);
02855 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
02856 fprintf(stderr,"----- Wrote bucket dataset ");
02857 fprintf(stderr,"into %s\n", DSET_BRIKNAME(new_dset));
02858
02859
02860
02861 THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
02862 free (volume); volume = NULL;
02863
02864 }
02865
02866
02867
02868
02869
02870
02871
02872 void output_results
02873 (
02874 matrix xdata,
02875 model * regmodel,
02876 RA_options * option_data
02877 )
02878
02879 {
02880 int p;
02881 int q;
02882 int n;
02883 int nxyz;
02884 int ip, ix;
02885 int numdof, dendof;
02886 int piece_size;
02887 int num_pieces;
02888 float * volume1 = NULL;
02889 float * volume2 = NULL;
02890 char filename[MAX_NAME_LENGTH];
02891
02892
02893
02894 p = regmodel->p;
02895 q = regmodel->q;
02896 n = xdata.rows;
02897 nxyz = option_data->nxyz;
02898 piece_size = option_data->piece_size;
02899 num_pieces = option_data->num_pieces;
02900
02901
02902
02903 volume1 = (float *) malloc (sizeof(float) * nxyz);
02904 MTEST (volume1);
02905 volume2 = (float *) malloc (sizeof(float) * nxyz);
02906 MTEST (volume2);
02907
02908
02909
02910 if (option_data->numt > 0)
02911 {
02912 numdof = n - p;
02913 dendof = 0;
02914
02915 for (ip = 0; ip < p; ip++)
02916 {
02917 ix = regmodel->flist[ip];
02918
02919 if (option_data->tcoef_filename[ix] != NULL)
02920 {
02921 sprintf (filename, "coef.%d", ix);
02922 read_volume (filename, volume1, nxyz, piece_size, num_pieces);
02923
02924 sprintf (filename, "tcoef.%d", ix);
02925 read_volume (filename, volume2, nxyz, piece_size, num_pieces);
02926
02927 write_afni_data (option_data,
02928 option_data->tcoef_filename[ix],
02929 volume1, volume2,
02930 FUNC_TT_TYPE, numdof, dendof);
02931 }
02932 }
02933 }
02934
02935
02936
02937 if (option_data->numr > 0)
02938 {
02939 strcpy (filename, "rsqr");
02940 read_volume (filename, volume2, nxyz, piece_size, num_pieces);
02941
02942 for (ip = 0; ip < p; ip++)
02943 {
02944 ix = regmodel->flist[ip];
02945
02946 if (option_data->rcoef_filename[ix] != NULL)
02947 {
02948 sprintf (filename, "coef.%d", ix);
02949 read_volume (filename, volume1, nxyz, piece_size, num_pieces);
02950
02951 write_afni_data (option_data,
02952 option_data->rcoef_filename[ix],
02953 volume1, volume2,
02954 FUNC_THR_TYPE, 0, 0);
02955 }
02956 }
02957 }
02958
02959
02960
02961 if (option_data->numf > 0)
02962 {
02963 strcpy (filename, "freg");
02964 read_volume (filename, volume2, nxyz, piece_size, num_pieces);
02965
02966 numdof = p - q;
02967 dendof = n - p;
02968
02969 for (ip = 0; ip < p; ip++)
02970 {
02971 ix = regmodel->flist[ip];
02972
02973 if (option_data->fcoef_filename[ix] != NULL)
02974 {
02975 sprintf (filename, "coef.%d", ix);
02976 read_volume (filename, volume1, nxyz, piece_size, num_pieces);
02977
02978 write_afni_data (option_data,
02979 option_data->fcoef_filename[ix],
02980 volume1, volume2,
02981 FUNC_FT_TYPE, numdof, dendof);
02982 }
02983 }
02984 }
02985
02986
02987
02988 free (volume1); volume1 = NULL;
02989 free (volume2); volume2 = NULL;
02990
02991
02992
02993 if (option_data->numbricks > 0)
02994 write_bucket_data (xdata, regmodel, option_data);
02995
02996
02997 }
02998
02999
03000
03001
03002
03003
03004
03005 void terminate_program
03006 (
03007 matrix * xdata,
03008 model * regmodel,
03009 RA_options * option_data
03010 )
03011
03012 {
03013 int p;
03014 int ip, ix;
03015 int ib;
03016 int nxyz;
03017 int piece_size;
03018 int num_pieces;
03019 char filename[MAX_NAME_LENGTH];
03020
03021
03022
03023 p = regmodel->p;
03024 nxyz = option_data->nxyz;
03025 piece_size = option_data->piece_size;
03026 num_pieces = option_data->num_pieces;
03027
03028
03029
03030 if (option_data->numf > 0)
03031 {
03032 strcpy (filename, "freg");
03033 delete_volume (filename, nxyz, piece_size, num_pieces);
03034 }
03035
03036
03037
03038 if (option_data->numr > 0)
03039 {
03040 strcpy (filename, "rsqr");
03041 delete_volume (filename, nxyz, piece_size, num_pieces);
03042 }
03043
03044
03045
03046 if (option_data->numt > 0)
03047 {
03048 for (ip = 0; ip < p; ip++)
03049 {
03050 ix = regmodel->flist[ip];
03051 sprintf (filename, "tcoef.%d", ix);
03052 delete_volume (filename, nxyz, piece_size, num_pieces);
03053 }
03054 }
03055
03056
03057
03058 if (option_data->numc > 0)
03059 {
03060 for (ip = 0; ip < p; ip++)
03061 {
03062 ix = regmodel->flist[ip];
03063 sprintf (filename, "coef.%d", ix);
03064 delete_volume (filename, nxyz, piece_size, num_pieces);
03065 }
03066 }
03067
03068
03069
03070 matrix_destroy (xdata);
03071
03072
03073
03074 if (option_data->yname != NULL)
03075 {
03076 for (ip = 0; ip < MAX_OBSERVATIONS; ip++)
03077 {
03078 if (option_data->yname[ip] != NULL)
03079 {
03080 free (option_data->yname[ip]);
03081 option_data->yname[ip] = NULL;
03082 }
03083 }
03084 free (option_data->yname);
03085 option_data->yname = NULL;
03086 }
03087
03088 if (option_data->fcoef_filename != NULL)
03089 {
03090 for (ip = 0; ip < MAX_XVARS; ip++)
03091 {
03092 if (option_data->fcoef_filename[ip] != NULL)
03093 {
03094 free (option_data->fcoef_filename[ip]);
03095 option_data->fcoef_filename[ip] = NULL;
03096 }
03097 }
03098 free (option_data->fcoef_filename);
03099 option_data->fcoef_filename = NULL;
03100 }
03101
03102 if (option_data->rcoef_filename != NULL)
03103 {
03104 for (ip = 0; ip < MAX_XVARS; ip++)
03105 {
03106 if (option_data->rcoef_filename[ip] != NULL)
03107 {
03108 free (option_data->rcoef_filename[ip]);
03109 option_data->rcoef_filename[ip] = NULL;
03110 }
03111 }
03112 free (option_data->rcoef_filename);
03113 option_data->rcoef_filename = NULL;
03114 }
03115
03116 if (option_data->tcoef_filename != NULL)
03117 {
03118 for (ip = 0; ip < MAX_XVARS; ip++)
03119 {
03120 if (option_data->tcoef_filename[ip] != NULL)
03121 {
03122 free (option_data->tcoef_filename[ip]);
03123 option_data->tcoef_filename[ip] = NULL;
03124 }
03125 }
03126 free (option_data->tcoef_filename);
03127 option_data->tcoef_filename = NULL;
03128 }
03129
03130 if (option_data->levels != NULL)
03131 {
03132 free (option_data->levels);
03133 option_data->levels = NULL;
03134 }
03135
03136 if (option_data->counts != NULL)
03137 {
03138 free (option_data->counts);
03139 option_data->counts = NULL;
03140 }
03141
03142
03143 if (option_data->bucket_filename != NULL)
03144 {
03145 free (option_data->bucket_filename);
03146 option_data->bucket_filename = NULL;
03147 }
03148
03149 if (option_data->brick_type != NULL)
03150 {
03151 free (option_data->brick_type);
03152 option_data->brick_type = NULL;
03153 }
03154
03155 if (option_data->brick_coef != NULL)
03156 {
03157 free (option_data->brick_coef);
03158 option_data->brick_coef = NULL;
03159 }
03160
03161 if (option_data->brick_label != NULL)
03162 {
03163 for (ib = 0; ib < option_data->numbricks; ib++)
03164 {
03165 if (option_data->brick_label[ib] != NULL)
03166 {
03167 free (option_data->brick_label[ib]);
03168 option_data->brick_label[ib] = NULL;
03169 }
03170 }
03171 free (option_data->brick_label);
03172 option_data->brick_label = NULL;
03173 }
03174
03175 }
03176
03177
03178
03179
03180
03181
03182
03183 int main
03184 (
03185 int argc,
03186 char ** argv
03187 )
03188
03189 {
03190 RA_options option_data;
03191 matrix xdata;
03192 model regmodel;
03193 int piece_size;
03194
03195
03196
03197 printf ("\n\n");
03198 printf ("Program: %s \n", PROGRAM_NAME);
03199 printf ("Author: %s \n", PROGRAM_AUTHOR);
03200 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
03201 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
03202 printf ("\n");
03203
03204
03205
03206 machdep() ;
03207 { int new_argc ; char ** new_argv ;
03208 addto_args( argc , argv , &new_argc , &new_argv ) ;
03209 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
03210 }
03211
03212
03213 initialize_program (argc, argv, &xdata, ®model, &option_data);
03214
03215
03216
03217 calculate_results (xdata, ®model, &option_data);
03218
03219
03220
03221 output_results (xdata, ®model, &option_data);
03222
03223
03224
03225 terminate_program (&xdata, ®model, &option_data);
03226
03227 exit(0);
03228 }
03229
03230
03231
