00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1998-2001, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   Plugin to calculate the deconvolution of a measured 3d+time dataset
00010   with a specified input stimulus time series.  Output consists of the
00011   least squares estimate of the impulse response function, and the fitted 
00012   time series (input stimulus convolved with estimated impulse response
00013   function).
00014 
00015   File:    plug_deconvolve.c
00016   Author:  B. Douglas Ward
00017   Date:    09 September 1998
00018 
00019   Mod:     Generate error message for time series of insufficient length.
00020   Date:    29 October 1998
00021 
00022   Mod:     Restructured matrix calculations to improve execution speed.
00023   Date:    16 December 1998
00024 
00025   Mod:     Minor correction to stim_label option.
00026   Date:    17 December 1998
00027 
00028   Mod:     Removed restriction on length of input time series.
00029   Date:    31 December 1998
00030 
00031   Mod:     Accept mean square error from full model.
00032   Date:    04 January 1999
00033 
00034   Mod:     Earlier termination if unable to invert X matrix.
00035            (Avoids redundant error messages.)
00036   Date:    06 January 1999
00037 
00038   Mod:     Change NLast default value.
00039   Date:    27 May 1999
00040 
00041   Mod:     Added option for matrix calculation of general linear tests.
00042   Date:    02 July 1999
00043 
00044   Mod:     Increased max. allowed number of input stimulus functions.
00045   Date:    24 August 1999
00046 
00047   Mod:     Additional statistical output (partial R^2 statistics).
00048   Date:    07 September 1999
00049 
00050   Mod:     Allow reading of multiple input stimulus functions from a single
00051            file by selection of individual columns.
00052   Date:    09 November 1999
00053 
00054   Mod:     Modifications for compatibility with 3dDeconvolve options for
00055            writing the fitted full model time series (-fitts) and the 
00056            residual error time series (-errts) to 3d+time datasets.
00057   Date:    22 November 1999
00058 
00059   Mod:     Added test for maximum number of full model parameters.
00060   Date:    24 November 1999
00061 
00062   Mod:     Small change to order of print out to screen.
00063   Date:    29 November 1999
00064 
00065   Mod:     Increased maximum number of GLT matrices and linear constraints.
00066   Date:    03 January 2000
00067 
00068   Mod:     Modified matrix_file_read to use mri_read_ascii routine.
00069   Date:    12 January 2000
00070 
00071   Mod:     Added censor input option to allow operator to eliminate individual
00072            time points from the analysis.
00073   Date:    21 June 2000
00074 
00075   Mod:     Added screen display of p-values.
00076   Date:    22 June 2000
00077 
00078   Mod:     Added -glt_label option for labeling the GLT matrix.
00079   Date:    23 June 2000
00080 
00081   Mod:     Added Concat Runs option for analysis of a concatenated 3d+time 
00082            dataset.
00083   Date:    27 June 2000
00084 
00085   Mod:     Increased size of screen output buffer.
00086   Date:    27 July 2000
00087 
00088   Mod:     Increased maximum number of linear constraints.
00089   Date:    07 December 2000
00090 
00091   Mod:     Changes required for compatibility with -nodata option 
00092            (norm.std.dev.'s for GLT linear constraints).
00093   Date:    21 December 2000
00094 
00095   Mod:     Added NPTR option, to allow input stim functions to be sampled
00096            at a multiple of the 1/TR rate.
00097   Date:    02 January 2001 
00098 
00099   Mod:     Increased maximum degree of polynomial ort.
00100   Date:    16 May 2001
00101 
00102   Mod:     Removed automatic override of NFirst option.
00103   Date:    08 June 2001
00104 
00105   Mod:     Enhanced screen output:  Display of p-values for individual stim
00106            function regression coefficients.  Display of t-stats and p-values
00107            for individual linear constraints within a GLT.
00108   Date:    29 January 2002
00109 
00110   Mod:     Allow user to specify no baseline parameters in the model 
00111            by setting Baseline to "None".           
00112   Date:    26 February 2002
00113 
00114   Mod:     Allow user to specify which input stimulus functions are part of
00115            the baseline model.
00116   Date:    02 May 2002
00117 
00118   Mod:     Improved graphical representation of the estimated impulse response
00119            function when the user selects option DC_IRF under Tran 1D of the 
00120            graph options menu.
00121   Date:    06 May 2002
00122 
00123   Mod:     Fixed bug in show_options rountine that would cause program crash 
00124            when the baseline was set to "None".
00125   Date:    03 Oct 2002
00126 
00127   Mod:     Increased size of screen output buffer (again).
00128   Date:    28 October 2002
00129 
00130   Mod:     Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
00131   Date:    02 December 2002
00132 
00133   Mod:     Additional input error testing for -censor and -concat options.
00134   Date:    17 March 2003
00135 */
00136 
00137 /*---------------------------------------------------------------------------*/
00138 
00139 #define PROGRAM_NAME    "plug_deconvolve"            /* name of this program */
00140 #define PROGRAM_AUTHOR  "B. Douglas Ward"                  /* program author */
00141 #define PROGRAM_INITIAL "09 September 1998"  /* initial program release date */
00142 #define PROGRAM_LATEST  "18 March 2003"      /* latest program revision date */
00143 
00144 /*---------------------------------------------------------------------------*/
00145 
00146 #include "afni.h"
00147 #include "matrix.h"
00148 
00149 #define MAX_NAME_LENGTH THD_MAX_NAME    /* max. string length for file names */
00150 #define MAX_XVARS 250                           /* max. number of parameters */
00151 #define MAX_STIMTS 20                 /* max. number of stimulus time series */
00152 #define MAX_GLT 20                    /* max. number of general linear tests */
00153 #define MAX_CONSTR 20                 /* max. number of linear constraints   */
00154 
00155 #define RA_error DC_error
00156 
00157 /*---------------------------------------------------------------------------*/
00158 /*
00159    Print error message and return.
00160 */
00161 
00162 static void DC_error (char * message)
00163 {
00164   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00165 
00166 }
00167 
00168 
00169 /*---------------------------------------------------------------------------*/
00170 
00171 #ifndef ALLOW_PLUGINS
00172 #  error "Plugins not properly set up -- see machdep.h"
00173 #endif
00174 
00175 
00176 /*------------- string to 'help' the user -------------*/
00177 
00178 static char helpstring[] =
00179    " Purpose: Control DC_Fit, DC_Err, and DC_IRF  Deconvolution Functions.  \n"
00180    "                                                                        \n"
00181    " Control      Base      = Baseline polynomial: None, Constant, Linear,  \n"
00182    "                            Quadratic, Cubic, Quartic, or Quintic       \n"
00183    "              NFirst    = Number of first time series point to use      \n"
00184    "              NLast     = Number of last time series point to use       \n"
00185    "              IRF       = Label of stimulus fnc. to use for IRF plot    \n"
00186    "                                                                        \n"
00187    " Concat       Label     = Name to use as label for concatenation        \n"
00188    "              File      = File containing list of volume indices of     \n"
00189    "                            starting points for individual runs         \n"
00190    "              Col #     = Column of file which contains concat run list \n"
00191    "                                                                        \n"
00192    " Censor       Label     = Name to use as label for censor function      \n"
00193    "              File      = Time series file indicating censored points   \n"
00194    "              Col #     = Column of file which contains censor function \n"
00195    "                                                                        \n"
00196    " StimFnc      Label     = Name to use as label for this input stimulus  \n"
00197    "              File      = Time series file representing input stimulus  \n"
00198    "              Col #     = Column of file which contains input stimulus  \n"
00199    "              MinLag    = Minimum time delay for impulse response       \n"
00200    "              MaxLag    = Maximum time delay for impulse response       \n"
00201    "              NPTR      = Number of stim fn. time points per TR         \n"
00202    "              Base      = Is this input stimulus part of baseline model?\n"
00203    "                                                                        \n"
00204    " GLT Mat      Label     = Name to use as label for this GLT matrix      \n"
00205    "              File      = File containing the GLT matrix                \n"
00206    "              # Rows    = Number of rows (linear constraints) in GLT    \n"
00207    "                                                                        \n"
00208    "                                                                        \n"
00209    " For more details, see:                                                 \n"
00210    " `Deconvolution Analysis of FMRI Time Series Data' by B. Douglas Ward,  \n"
00211    " which is contained in file 3dDeconvolve.ps of the AFNI distribution.   \n"
00212 ;
00213 
00214 
00215 /*------- Strings for baseline control ------*/
00216 
00217 #define NBASE 7
00218 static char * baseline_strings[NBASE] = {"None", "Const", "Linear", 
00219                             "Quadrtc", "Cubic", "Quartic", "Quintic" };
00220 
00221 static char * false_or_true[2] = {"False", "True"};
00222 
00223 /*--------------- prototypes for internal routines ---------------*/
00224 
00225 static char * DC_main( PLUGIN_interface * ) ;  /* the entry point */
00226 
00227 static void DC_Fit (int nt, double to, double dt, float * vec, char ** label) ;
00228 static void DC_Err (int nt, double to, double dt, float * vec, char ** label) ;
00229 static void DC_IRF (int nt, double to, double dt, float * vec, char ** label) ;
00230 static int calculate_results();
00231 
00232 
00233 
00234 /*---------------- global data -------------------*/
00235 
00236 static PLUGIN_interface * global_plint = NULL;
00237 
00238 static int plug_polort;     /* degree of polynomial for baseline model */
00239 static int plug_p;          /* total number of parameters in the full model */
00240 static int plug_q;          /* total number of parameters in the base model */
00241 static int plug_qp;         /* number of poly. trend baseline parameters */
00242 static int plug_NFirst;     /* first image from input 3d+time dataset to use */
00243 static int plug_NLast;      /* last image from input 3d+time dataset to use */
00244 static int plug_IRF;        /* which impulse response fuction to plot */
00245 static int initialize;      /* flag for perform initialization */
00246 static int prev_nt;         /* previous time series length */
00247 static char * IRF_label;    /* label of stimulus for IRF plot */
00248 
00249 static char * concat_label;      /* label for concatenation */
00250 static int concat_column;        /* column containing list of blocks (runs) */
00251 static int num_blocks;           /* number of blocks (runs) */
00252 static int * block_list;         /* list of block (run) starting points */
00253 
00254 static int num_censor;                /* flag for include censor function */
00255 static char * censor_label;           /* censor time series label */
00256 static int censor_column;             /* column containing censor array */
00257 static float * censor_array;          /* censor time series array */
00258 static int censor_length;             /* length of censor time series */
00259 static int * good_list;               /* list of usable time points */
00260 
00261 static int num_stimts;                     /* number of stimulus time series */
00262 static int stim_base[MAX_STIMTS];          /* flag for baseline stimulus */
00263 static int stim_column[MAX_STIMTS];        /* column containing stimulus */
00264 static float * stimulus[MAX_STIMTS];       /* stimulus time series arrays */
00265 static int stim_length[MAX_STIMTS];        /* length of stimulus time series */
00266 static int min_lag[MAX_STIMTS];   /* minimum time delay for impulse response */
00267 static int max_lag[MAX_STIMTS];   /* maximum time delay for impulse response */
00268 static int nptr[MAX_STIMTS];      /* number of stim fn. time points per TR */
00269 static char * stim_label[MAX_STIMTS];      /* stimulus time series labels */
00270 
00271 static matrix xdata;              /* independent variable matrix */
00272 static matrix x_full;             /* extracted X matrix   for full model */
00273 static matrix xtxinv_full;        /* matrix:  1/(X'X)     for full model */
00274 static matrix xtxinvxt_full;      /* matrix:  (1/(X'X))X' for full model */
00275 static matrix x_base;             /* extracted X matrix   for baseline model */
00276 static matrix xtxinvxt_base;      /* matrix:  (1/(X'X))X' for baseline model */
00277 static matrix x_rdcd[MAX_STIMTS]; /* extracted X matrices for reduced models */
00278 static matrix xtxinvxt_rdcd[MAX_STIMTS];     
00279                                   /* matrix:  (1/(X'X))X' for reduced models */
00280 
00281 static int glt_num;                 /* number of general linear tests */
00282 static char * glt_label[MAX_GLT];   /* general linear test labels */
00283 static int glt_rows[MAX_GLT];       /* number of linear constraints in glt */
00284 static char * glt_filename[MAX_GLT];/* file containing glt matrix */
00285 static matrix cxtxinvct[MAX_GLT];   /* matrices: C(1/(X'X))C' for GLT */
00286 static matrix glt_cmat[MAX_GLT];    /* general linear test matrices */
00287 static matrix glt_amat[MAX_GLT];    /* constant GLT matrices for later use */
00288 static vector glt_coef[MAX_GLT];    /* linear combinations from GLT matrices */
00289 static vector glt_tcoef[MAX_GLT];   /* t-stats for GLT linear combinations */
00290 
00291 
00292 /*---------------------------------------------------------------------------*/
00293 /*
00294   Include deconvolution analysis software.
00295 */
00296 
00297 #undef USE_BASIS          /* 11 Aug 2004: for Deconvolve.c */
00298 
00299 #include "Deconvolve.c"
00300 
00301 
00302 
00303 /*---------------------------------------------------------------------------*/
00304 /*
00305   Initialize control options and global data.
00306 */
00307 
00308 static void initialize_options ()
00309 {
00310   int is;                     /* input stimulus index */
00311   int iglt;                   /* index for general linear test */
00312 
00313 
00314   /*----- Initialize control parameters -----*/
00315   plug_polort = 1;        /* degree of polynomial for baseline model */
00316   plug_p      = 0;        /* total number of parameters in the full model */
00317   plug_q      = 0;        /* total number of parameters in the base model */
00318   plug_qp     = 0;        /* number of poly. trend baseline parameters */
00319   plug_NFirst = 0;        /* first image from input 3d+time dataset to use */
00320   plug_NLast  = 32767;    /* last image from input 3d+time dataset to use */
00321   plug_IRF    = -1;       /* which impulse response fuction to plot */
00322   initialize  = 0;        /* flag for perform initialization */
00323   prev_nt     = 0;        /* previous time series length */
00324   IRF_label   = malloc (sizeof(char)*MAX_NAME_LENGTH);   MTEST (IRF_label);
00325   strcpy (IRF_label, " ");      /* label of stimulus for IRF plot */
00326 
00327 
00328   /*----- Initialization for concatenated runs -----*/
00329   concat_label = malloc (sizeof(char)*MAX_NAME_LENGTH);  MTEST (concat_label);
00330   strcpy (concat_label, " ");    /* label for concatenation */
00331   concat_column = 0;             /* column containing list of blocks (runs) */
00332   num_blocks = 1;                /* number of blocks (runs) */
00333   block_list = (int *) malloc (sizeof(int) * 1);  MTEST (block_list);  
00334   block_list[0] = 0;             /* list of block (run) starting points */
00335 
00336 
00337   /*----- Initialize censorship function -----*/
00338   num_censor = 0;
00339   censor_label = malloc (sizeof(char)*MAX_NAME_LENGTH);  MTEST (censor_label);
00340   strcpy (censor_label, " ");        /* censor time series label */
00341   censor_column = 0;                 /* column containing censor array */
00342   censor_array = NULL;               /* censor time series array */
00343   censor_length = 0;                 /* length of censor time series */
00344   good_list = NULL;                  /* list of usable time points */
00345 
00346 
00347   /*----- Initialize stimulus functions -----*/
00348   num_stimts = 0;                          /* number of stimulus time series */
00349   for (is =0;  is < MAX_STIMTS;  is++)
00350     {
00351       stim_label[is] = malloc (sizeof(char)*MAX_NAME_LENGTH);
00352       MTEST (stim_label[is]);
00353       sprintf (stim_label[is], "Stim #%d ", is+1);
00354                                            /* stimulus time series labels */
00355       stim_base[is] = 0;                   /* flag for baseline stimulus */
00356       stim_column[is] = 0;                 /* column containing stimulus */
00357       stimulus[is] = NULL;                 /* stimulus time series arrays */
00358       stim_length[is] = 0;                 /* length of stimulus time series */
00359       min_lag[is] = 0;            /* minimum time delay for impulse response */
00360       max_lag[is] = 0;            /* maximum time delay for impulse response */
00361       nptr[is] = 1;               /* number of stim fn. time points per TR */
00362    }
00363 
00364 
00365   /*----- Initialize matrices -----*/
00366   matrix_initialize (&xdata);         /* independent variable matrix */
00367   matrix_initialize (&x_full);        /* extracted X matrix   for full model */
00368   matrix_initialize (&xtxinv_full);   /* matrix:  1/(X'X)     for full model */
00369   matrix_initialize (&xtxinvxt_full); /* matrix:  (1/(X'X))X' for full model */
00370   matrix_initialize (&x_base);    /* extracted X matrix   for baseline model */
00371   matrix_initialize (&xtxinvxt_base);
00372                                   /* matrix:  (1/(X'X))X' for baseline model */
00373   for (is =0;  is < MAX_STIMTS;  is++)
00374     {
00375       matrix_initialize (&x_rdcd[is]); 
00376                                   /* extracted X matrices for reduced models */
00377       matrix_initialize (&xtxinvxt_rdcd[is]);
00378                                   /* matrix:  (1/(X'X))X' for reduced models */
00379     }
00380 
00381 
00382   /*----- Initialize GLT matrices -----*/
00383   glt_num = 0;                             /* number of general linear tests */
00384   for (iglt =0;  iglt < MAX_GLT;  iglt++)
00385     {
00386       glt_label[iglt] = malloc (sizeof(char)*MAX_NAME_LENGTH);
00387       MTEST (glt_label[iglt]);
00388       sprintf (glt_label[iglt], "GLT #%d ", iglt+1);
00389                                                /* general linear test labels */
00390       glt_rows[iglt] = 0;             /* number of linear constraints in glt */
00391       glt_filename[iglt] = malloc (sizeof(char)*MAX_NAME_LENGTH);
00392       MTEST (glt_filename[iglt]);
00393       strcpy (glt_filename[iglt], " ");        /* file containing glt matrix */
00394 
00395       matrix_initialize (&cxtxinvct[iglt]);
00396                                            /* matrices: C(1/(X'X))C' for GLT */
00397       matrix_initialize (&glt_cmat[iglt]);
00398                                              /* general linear test matrices */
00399       matrix_initialize (&glt_amat[iglt]); 
00400                                       /* constant GLT matrices for later use */
00401       vector_initialize (&glt_coef[iglt]);
00402                                     /* linear combinations from GLT matrices */
00403       vector_initialize (&glt_tcoef[iglt]);
00404                                     /* t-stats for GLT linear combinations   */
00405     }
00406 
00407 }
00408 
00409 
00410 /*---------------------------------------------------------------------------*/
00411 /*
00412   Reset control options and global data.
00413 */
00414 
00415 static void reset_options ()
00416 {
00417   int is;                     /* input stimulus index */
00418   int iglt;                   /* index for general linear test */
00419 
00420 
00421   /*----- Reset control parameters -----*/
00422   plug_polort = 1;        /* degree of polynomial for baseline model */
00423   plug_p      = 0;        /* total number of parameters in the full model */
00424   plug_q      = 0;        /* total number of parameters in the base model */
00425   plug_qp     = 0;        /* number of poly. trend baseline parameters */
00426   plug_NFirst = 0;        /* first image from input 3d+time dataset to use */
00427   plug_NLast  = 32767;    /* last image from input 3d+time dataset to use */
00428   plug_IRF    = -1;       /* which impulse response fuction to plot */
00429   initialize  = 0;        /* flag for perform initialization */
00430   prev_nt     = 0;        /* previous time series length */
00431   strcpy (IRF_label, " ");       /* label of stimulus for IRF plot */
00432 
00433 
00434   /*----- Reset for concatenated runs -----*/
00435   strcpy (concat_label, " ");    /* label for concatenation */
00436   concat_column = 0;             /* column containing list of blocks (runs) */
00437   num_blocks = 1;                /* number of blocks (runs) */
00438   if (block_list != NULL)  free (block_list);       
00439   block_list = (int *) malloc (sizeof(int) * 1);  MTEST (block_list);  
00440   block_list[0] = 0;             /* list of block (run) starting points */
00441 
00442 
00443   /*----- Reset censorship function -----*/
00444   num_censor = 0;
00445   strcpy (censor_label, " ");        /* censor time series label */
00446   censor_column = 0;                 /* column containing censor array */
00447   if (censor_array != NULL)
00448     {  
00449       free (censor_array);   
00450       censor_array = NULL;           /* censor time series array */
00451     }
00452   censor_length = 0;                 /* length of censor time series */
00453   if (good_list != NULL)
00454     {
00455       free (good_list);
00456       good_list = NULL;              /* list of usable time points */
00457     }
00458 
00459 
00460   /*----- Reset stimulus functions -----*/
00461   num_stimts = 0;                          /* number of stimulus time series */
00462   for (is =0;  is < MAX_STIMTS;  is++)
00463     {
00464       sprintf (stim_label[is], "Stim #%d ", is+1);
00465                                            /* stimulus time series labels */
00466       stim_base[is] = 0;                   /* flag for baseline stimulus */
00467       stim_column[is] = 0;                 /* column containing stimulus */
00468       if (stimulus[is] != NULL)
00469         {
00470           free (stimulus[is]);
00471           stimulus[is] = NULL;             /* stimulus time series arrays */
00472         }
00473       stim_length[is] = 0;                 /* length of stimulus time series */
00474       min_lag[is] = 0;            /* minimum time delay for impulse response */
00475       max_lag[is] = 0;            /* maximum time delay for impulse response */
00476       nptr[is] = 1;               /* number of stim fn. time points per TR */
00477    }
00478 
00479 
00480   /*----- Destroy matrices -----*/
00481   matrix_destroy (&xdata);         /* independent variable matrix */
00482   matrix_destroy (&x_full);        /* extracted X matrix   for full model */
00483   matrix_destroy (&xtxinv_full);   /* matrix:  1/(X'X)     for full model */
00484   matrix_destroy (&xtxinvxt_full); /* matrix:  (1/(X'X))X' for full model */
00485   matrix_destroy (&x_base);       /* extracted X matrix   for baseline model */
00486   matrix_destroy (&xtxinvxt_base);
00487                                   /* matrix:  (1/(X'X))X' for baseline model */
00488   for (is =0;  is < MAX_STIMTS;  is++)
00489     {
00490       matrix_destroy (&x_rdcd[is]); 
00491                                   /* extracted X matrices for reduced models */
00492       matrix_destroy (&xtxinvxt_rdcd[is]);
00493                                   /* matrix:  (1/(X'X))X' for reduced models */
00494     }
00495 
00496 
00497   /*----- Destroy GLT matrices -----*/
00498   glt_num = 0;                             /* number of general linear tests */
00499   for (iglt =0;  iglt < MAX_GLT;  iglt++)
00500     {
00501       sprintf (glt_label[iglt], "GLT #%d ", iglt+1);
00502                                                /* general linear test labels */
00503       glt_rows[iglt] = 0;             /* number of linear constraints in glt */
00504       strcpy (glt_filename[iglt], " ");        /* file containing glt matrix */
00505 
00506       matrix_destroy (&cxtxinvct[iglt]);
00507                                            /* matrices: C(1/(X'X))C' for GLT */
00508       matrix_destroy (&glt_cmat[iglt]);
00509                                              /* general linear test matrices */
00510       matrix_destroy (&glt_amat[iglt]); 
00511                                       /* constant GLT matrices for later use */
00512       vector_destroy (&glt_coef[iglt]);
00513                                     /* linear combinations from GLT matrices */
00514       vector_destroy (&glt_tcoef[iglt]);
00515                                     /* t-stats for GLT linear combinations   */
00516     }
00517 
00518 }
00519 
00520 
00521 /*---------------------------------------------------------------------------*/
00522 /*
00523    Set up the interface to the user.
00524 */
00525 
00526 
00527 DEFINE_PLUGIN_PROTOTYPE
00528 
00529 PLUGIN_interface * PLUGIN_init( int ncall )
00530 {
00531    PLUGIN_interface * plint ;     /* will be the output of this routine */
00532    int is;                        /* input stimulus index */
00533    int iglt;                   /* index for general linear test */
00534 
00535 
00536    if( ncall > 0 ) return NULL ;  /* generate interface for ncall 0 */
00537 
00538 
00539    /***** do interface #0 *****/
00540 
00541    /*---------------- set titles and call point ----------------*/
00542 
00543    plint = PLUTO_new_interface ("Deconvolution" ,
00544            "Control DC_Fit, DC_Err, and DC_IRF Deconvolution Functions" ,
00545            helpstring, PLUGIN_CALL_VIA_MENU, DC_main);
00546 
00547    global_plint = plint ;  /* make global copy */
00548 
00549    PLUTO_short_choose(plint) ;  /* 29 Mar 2002 [RWCox]: */
00550    PLUTO_short_number(plint) ;  /* make 'Choose' and number widgets shorter */
00551 
00552    PLUTO_add_hint (plint, 
00553      "Control DC_Fit, DC_Err, and DC_IRF Deconvolution Functions");
00554 
00555    PLUTO_set_sequence( plint , "A:funcs:fitting" ) ;
00556 
00557    PLUTO_set_runlabels( plint , "Set+Keep" , "Set+Close" ) ;  /* 04 Nov 2003 */
00558 
00559    /*----- Parameters -----*/
00560    PLUTO_add_option (plint, "Control", "Control", TRUE);
00561    PLUTO_add_string (plint, "Base", NBASE, baseline_strings, 2);
00562    PLUTO_add_number (plint, "NFirst", -1, 32767, 0, -1, TRUE);
00563    PLUTO_add_number (plint, "NLast",  0, 32767, 0, 32767,  TRUE);
00564    PLUTO_add_string( plint, "IRF ",    0, NULL, 1);
00565 
00566 
00567    /*----- Concatenation Function -----*/
00568    PLUTO_add_option (plint, "Concat", "Concat", FALSE);
00569    PLUTO_add_string( plint, "Label", 0, NULL, 1);
00570    PLUTO_add_timeseries (plint, "File");
00571    PLUTO_add_number (plint, "Col #", 0, 100, 0, 0, TRUE);
00572 
00573 
00574    /*----- Censor Function -----*/
00575    PLUTO_add_option (plint, "Censor", "Censor", FALSE);
00576    PLUTO_add_string( plint, "Label", 0, NULL, 1);
00577    PLUTO_add_timeseries (plint, "File");
00578    PLUTO_add_number (plint, "Col #", 0, 100, 0, 0, TRUE);
00579 
00580 
00581    /*----- Input Stimulus -----*/
00582    for (is = 0;  is < MAX_STIMTS;  is++)
00583      {
00584        PLUTO_add_option (plint, "StimFnc", "StimFnc", FALSE);
00585        PLUTO_add_string( plint, "Label", 0, NULL, 1);
00586        PLUTO_add_timeseries (plint, "File");
00587        PLUTO_add_number (plint, "Col #", 0, 100, 0, 0, TRUE);
00588        PLUTO_add_number (plint, "MinLag", 0, 100, 0, 0, TRUE);
00589        PLUTO_add_number (plint, "MaxLag", 0, 100, 0, 0, TRUE);
00590        PLUTO_add_number (plint, "NPTR",    1, 100, 0, 0, TRUE);
00591        PLUTO_add_string (plint, "Base", 2, false_or_true, 0);
00592      }
00593 
00594    /*----- General Linear Test -----*/
00595    for (is = 0;  is < MAX_GLT;  is++)
00596      {
00597        PLUTO_add_option (plint, "GLT Mat", "GLT Mat", FALSE);
00598        PLUTO_add_string( plint, "Label", 0, NULL, 1);
00599        PLUTO_add_string( plint, "File", 0, NULL, 1);     
00600        PLUTO_add_number (plint, "# Rows", 1, MAX_CONSTR, 0, 0, TRUE);
00601      }
00602 
00603    /*--------- done with interface setup ---------*/
00604    PLUTO_register_1D_funcstr ("DC_Fit" , DC_Fit);
00605    PLUTO_register_1D_funcstr ("DC_Err" , DC_Err);
00606    PLUTO_register_1D_funcstr ("DC_IRF" , DC_IRF);
00607    
00608 
00609    /*----- Initialize options and global data -----*/
00610    initialize_options ();
00611 
00612 
00613    return plint ;
00614 }
00615 
00616 
00617 /*---------------------------------------------------------------------------*/
00618 
00619 static void show_options ()
00620 {
00621   int ib;                         /* block (run) index */
00622   int it;                         /* time index */
00623   int is;                         /* stimulus index */
00624   int iglt;                       /* general linear test index */
00625 
00626 
00627   /*----- Identify software -----*/
00628   printf ("\n\n");
00629   printf ("Program:          %s \n", PROGRAM_NAME);
00630   printf ("Author:           %s \n", PROGRAM_AUTHOR); 
00631   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
00632   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
00633   printf ("\n");
00634 
00635 
00636   /*----- Show current input options -----*/
00637   printf ("\nControls: \n");
00638   printf ("Baseline  = %10s \n", baseline_strings[plug_polort+1]);
00639   printf ("NFirst    = %10d \n", plug_NFirst);
00640   printf ("NLast     = %10d \n", plug_NLast);
00641   printf ("IRF label = %10s \n", IRF_label);
00642 
00643 
00644   /*----- Identify concatenation function -----*/
00645   if (num_blocks > 0)
00646     {
00647       printf ("\n");
00648       printf ("Concatenation:     Label = %8s ", concat_label);
00649       printf ("Column = %3d  \n", concat_column);
00650       for (ib = 0;  ib < num_blocks;  ib++)
00651         printf ("Run #%d  Initial Point = %d \n", ib+1, block_list[ib]); 
00652     }
00653 
00654 
00655   /*----- Identify censor function -----*/
00656   if (num_censor > 0)
00657     {
00658       printf ("\n");
00659       printf ("Censor Function:   Label = %8s ", censor_label);
00660       printf ("Column = %3d  \n", censor_column);
00661       printf ("Censored Points: ");
00662       for (it = 0;  it < censor_length;  it++)
00663         {
00664           if (censor_array[it] == 0)  printf (" %d", it);
00665         }
00666       printf ("\n");
00667     }
00668 
00669 
00670   /*----- List stimulus functions -----*/
00671   if (num_stimts > 0)
00672     {
00673       printf ("\n");
00674       for (is = 0;  is < num_stimts;  is++)
00675         {
00676           if (stim_base[is])
00677             printf ("Baseline:      Label = %8s ", stim_label[is]);
00678           else
00679             printf ("Stimulus:      Label = %8s ", stim_label[is]);
00680           printf ("Column = %3d   Min. Lag = %3d   Max. Lag = %3d   ", 
00681                   stim_column[is], min_lag[is], max_lag[is]);
00682           printf ("NPTR = %d \n", nptr[is]);
00683         }
00684     }
00685 
00686 
00687   /*----- List GLT matrices -----*/
00688   if (glt_num > 0)
00689     {
00690       printf ("\n");
00691       for (iglt = 0;  iglt < glt_num;  iglt++)
00692         {
00693           printf ("GLT:       Label = %8s   ", glt_label[iglt]);
00694           printf ("#Rows = %2d   Input File: %s \n", 
00695                   glt_rows[iglt], glt_filename[iglt]);
00696         }
00697     }
00698  
00699 }
00700 
00701 
00702 /*---------------------------------------------------------------------------*/
00703 /*
00704   Main routine for this plugin (will be called from AFNI).
00705   If the return string is not NULL, some error transpired, and
00706   AFNI will popup the return string in a message box.
00707 */
00708 
00709 static char * DC_main( PLUGIN_interface * plint )
00710 {
00711   char * str;                           /* input string */
00712   int is;                               /* stimulus index */
00713   int iglt;                             /* general linear test index */
00714   MRI_IMAGE * stim;     /* pointers to image structures 
00715                            -- used to read 1D ASCII */
00716   float * far;          /* pointer to MRI_IMAGE floating point data */
00717   int ipt;              /* time point index */
00718   
00719 
00720   /*----- Reset options and global data -----*/
00721   reset_options ();
00722 
00723 
00724   /*--------- go to Control input line ---------*/
00725   PLUTO_next_option (plint);
00726   str    = PLUTO_get_string (plint);
00727   plug_polort = PLUTO_string_index (str, NBASE, baseline_strings) - 1;
00728   plug_NFirst = PLUTO_get_number (plint);
00729   plug_NLast  = PLUTO_get_number (plint);
00730   strcpy (IRF_label, PLUTO_get_string (plint));
00731           
00732 
00733   /*--------- go to Concat Runs input line ---------*/
00734   str = PLUTO_peek_optiontag (plint);
00735   if (str != NULL) 
00736 
00737     /*----- Read Concat Runs input Line -----*/
00738     if (strcmp (str, "Concat") == 0)
00739       {
00740         str = PLUTO_get_optiontag (plint);
00741         strcpy (concat_label, PLUTO_get_string(plint));
00742 
00743         stim = PLUTO_get_timeseries(plint) ;
00744         
00745         if (stim == NULL || stim->nx < 1 
00746             ||  stim->kind != MRI_float)
00747           return "**************************\n"
00748             "Illegal Concat File Input!\n"
00749             "**************************"  ;
00750         
00751         /*----- Column in file which contains the concat function -----*/
00752         concat_column = PLUTO_get_number(plint);
00753         
00754         if ((concat_column < 0) 
00755             ||(concat_column > stim->ny - 1))
00756           return "**********************************\n"
00757             "Illegal Concat File Column Number!\n"
00758             "**********************************"  ;
00759 
00760         /*----- Extract concat run start list from MRI data structure -----*/
00761         far = MRI_FLOAT_PTR(stim);
00762         num_blocks = stim->nx;
00763         if (block_list != NULL)  free (block_list);
00764         block_list = (int *) malloc (sizeof(int) * num_blocks);
00765         MTEST (block_list);
00766         
00767         for (ipt = 0;  ipt < num_blocks;  ipt++)
00768           block_list[ipt] = floor (far[ipt+concat_column*(stim->nx)] + 0.5); 
00769       }
00770   
00771 
00772   /*--------- go to Censor input line ---------*/
00773   str = PLUTO_peek_optiontag (plint);
00774   if (str != NULL) 
00775 
00776     /*----- Read Censor input Line -----*/
00777     if (strcmp (str, "Censor") == 0)
00778       {
00779         str = PLUTO_get_optiontag (plint);
00780         strcpy (censor_label, PLUTO_get_string(plint));
00781 
00782         stim = PLUTO_get_timeseries(plint) ;
00783         
00784         if (stim == NULL || stim->nx < 3 
00785             ||  stim->kind != MRI_float)
00786           return "**************************\n"
00787             "Illegal Censor File Input!\n"
00788             "**************************"  ;
00789 
00790         
00791         /*----- Column in file which contains the censor function -----*/
00792         censor_column = PLUTO_get_number(plint);
00793         
00794         if ((censor_column < 0) 
00795             ||(censor_column > stim->ny - 1))
00796           return "**********************************\n"
00797             "Illegal Censor File Column Number!\n"
00798             "**********************************"  ;
00799 
00800 
00801         /*----- Extract censor time series from MRI data structure -----*/
00802         if (censor_array != NULL) 
00803           {
00804             free (censor_array);
00805             censor_array = NULL;
00806           }
00807         far = MRI_FLOAT_PTR(stim);
00808         censor_length = stim->nx;
00809         censor_array = (float *) malloc (sizeof(float) * (stim->nx));
00810         MTEST (censor_array);
00811         
00812         num_censor = 1;
00813         for (ipt = 0;  ipt < (stim->nx);  ipt++)
00814           censor_array[ipt] 
00815             = far[ipt + censor_column*(stim->nx)]; 
00816       }
00817   
00818 
00819   /*------ Loop over input line(s) -----*/
00820   do
00821     {
00822       str = PLUTO_get_optiontag(plint); 
00823       if (str == NULL)  break;
00824       if ((strcmp (str, "StimFnc") != 0) && (strcmp (str, "GLT Mat") != 0))
00825         return "************************\n"
00826                "Illegal optiontag found!\n"
00827                "************************";
00828      
00829 
00830       /*----- Read Input Stimulus Line -----*/
00831       if (strcmp (str, "StimFnc") == 0)
00832         {      
00833           str =  PLUTO_get_string(plint);
00834           if (strlen(str) != 0)  strcpy (stim_label[num_stimts], str);
00835 
00836           if (strcmp(stim_label[num_stimts], IRF_label) == 0)
00837             plug_IRF = num_stimts;
00838           
00839           stim = PLUTO_get_timeseries(plint) ;
00840           
00841           if (stim == NULL || stim->nx < 3 
00842               ||  stim->kind != MRI_float)
00843             return "*************************\n"
00844                    "Illegal Timeseries Input!\n"
00845                    "*************************"  ;
00846 
00847 
00848           /*----- Column in file which contains the stimulus function -----*/
00849           stim_column[num_stimts] = PLUTO_get_number(plint);
00850 
00851           if ((stim_column[num_stimts] < 0) 
00852             ||(stim_column[num_stimts] > stim->ny - 1))
00853             return "********************************\n"
00854                    "Illegal Stim File Column Number!\n"
00855                    "********************************"  ;
00856 
00857 
00858           /*----- Extract stimulus time series from MRI data structure -----*/
00859           if (stimulus[num_stimts] != NULL) 
00860             {
00861               free (stimulus[num_stimts]);
00862               stimulus[num_stimts] = NULL;
00863             }
00864           far = MRI_FLOAT_PTR(stim);
00865           stim_length[num_stimts] = stim->nx;
00866           stimulus[num_stimts] = (float *) malloc (sizeof(float) * (stim->nx));
00867           MTEST (stimulus[num_stimts]);
00868 
00869           for (ipt = 0;  ipt < (stim->nx);  ipt++)
00870             stimulus[num_stimts][ipt] 
00871               = far[ipt + stim_column[num_stimts]*(stim->nx)]; 
00872 
00873 
00874           /*----- Minimum and Maximum time lags for model -----*/
00875           min_lag[num_stimts] = PLUTO_get_number(plint);
00876           max_lag[num_stimts] = PLUTO_get_number(plint);
00877           nptr[num_stimts]    = PLUTO_get_number(plint);
00878           str    = PLUTO_get_string (plint);
00879           stim_base[num_stimts] = PLUTO_string_index (str, 2, false_or_true);
00880 
00881           
00882           if (min_lag[num_stimts] > max_lag[num_stimts])
00883             return "**************************\n"
00884                    "Require Min Lag <= Max Lag\n"
00885                    "**************************"  ;
00886           
00887           num_stimts++;
00888         }
00889 
00890 
00891       /*----- Read General Matrix Test Line -----*/
00892       if (strcmp (str, "GLT Mat") == 0)
00893         {      
00894           str =  PLUTO_get_string(plint);
00895           if (strlen(str) != 0)  strcpy (glt_label[glt_num], str);
00896 
00897           strcpy (glt_filename[glt_num], PLUTO_get_string(plint));
00898 
00899           glt_rows[glt_num] = PLUTO_get_number(plint);
00900       
00901           glt_num++;
00902         }
00903 
00904     }
00905 
00906   while (1);
00907 
00908 
00909   /*----- Determine total number of parameters in the model -----*/
00910   plug_qp = (plug_polort + 1) * num_blocks;
00911   plug_q = plug_qp;
00912   plug_p = plug_qp;
00913   for (is = 0;  is < num_stimts;  is++)
00914     {
00915       if (stim_base[is])  plug_q += max_lag[is] - min_lag[is] + 1;
00916       plug_p += max_lag[is] - min_lag[is] + 1;
00917       if (plug_p > MAX_XVARS)
00918         { 
00919           return "****************************\n"
00920             "Too many parameters in model \n"
00921             "****************************"  ;
00922         }
00923     }
00924  
00925 
00926   /*----- Read the general linear test matrices -----*/
00927   if (glt_num > 0)
00928     for (iglt = 0;  iglt < glt_num;  iglt++)
00929       {
00930         matrix_file_read (glt_filename[iglt],
00931                           glt_rows[iglt],
00932                           plug_p,
00933                           &(glt_cmat[iglt]), 0);
00934         if (glt_cmat[iglt].elts == NULL)
00935           { 
00936             return "**************************\n"
00937                    "Unable to read GLT matrix \n"
00938                    "**************************"  ;
00939           }
00940       } 
00941 
00942 
00943   /*----- Show the user input options -----*/
00944   show_options ();
00945 
00946 
00947   /*--- nothing left to do until data arrives ---*/
00948   initialize = 1 ;  /* successful initialization */
00949   prev_nt = 0;      /* previous time series length */
00950   
00951   return NULL ;
00952 }
00953 
00954 
00955 /*---------------------------------------------------------------------------*/
00956 /*
00957   Calculate the impulse response function and associated statistics.
00958 */
00959 
00960 static int calculate_results 
00961 (
00962   int nt,               /* number of time points */
00963   double dt,            /* delta time */
00964   float * vec,          /* input measured data vector */
00965   int * NN,             /* number of usable data points */
00966   int * nfit,           /* number of fit parameters */
00967   float * fit,          /* fit parameters (regression coefficients) */
00968   char ** label,        /* string containing statistical summary of results */
00969   float ** fitts,       /* full model fitted time series */
00970   float ** errts        /* full model residual error time series */
00971 
00972 )
00973   
00974 {
00975   float * ts_array;           /* array of measured data for one voxel */
00976 
00977   int N;                      /* number of usable data points */
00978   int p;                      /* number of parameters in the full model */
00979   int q;                      /* number of parameters in the baseline model */
00980   int qp;                     /* number of poly. trend baseline parameters */
00981   int m;                      /* parameter index */
00982   int n;                      /* data point index */
00983 
00984   vector coef;                /* regression parameters */
00985   vector scoef;               /* std. devs. for regression parameters */
00986   vector tcoef;               /* t-statistics for regression parameters */
00987   float fpart[MAX_STIMTS];    /* partial F-statistics for the stimuli */
00988   float rpart[MAX_STIMTS];    /* partial R^2 stats. for the stimuli */
00989   float ffull;                /* full model F-statistic */
00990   float rfull;                /* full model R^2 statistic */
00991   float mse;                  /* mean square error from full model */
00992 
00993   vector y;                   /* vector of measured data */       
00994 
00995   int NFirst;            /* first image from input 3d+time dataset to use */
00996   int NLast;             /* last image from input 3d+time dataset to use */
00997   int i, it;             /* data point index */
00998   int is;                /* stimulus index */
00999 
01000   float rms_min = 0.0;   /* minimum variation in data to fit full model */
01001   int ok;                /* flag for successful matrix calculation */
01002   int novar;             /* flag for insufficient variation in data */
01003   float fglt[MAX_GLT];   /* F-statistics for the general linear tests */
01004   float rglt[MAX_GLT];   /* R^2 statistics for the general linear tests */
01005 
01006   int ib;                /* block (run) index */
01007   int irb;               /* time index relative to start of block (run) */
01008 
01009 
01010   /*----- Check initialization flag -----*/
01011   if (! initialize)  return (0);
01012 
01013 
01014   /*----- Initialize vectors -----*/
01015   vector_initialize (&coef);
01016   vector_initialize (&scoef);
01017   vector_initialize (&tcoef);
01018   vector_initialize (&y);
01019   
01020 
01021   /*----- Initialize local variables -----*/
01022   qp = plug_qp;
01023   q = plug_q;
01024   p = plug_p;
01025   *nfit = p;
01026 
01027 
01028   /*----- Allocate memory for fitted time series and residuals -----*/
01029   *fitts    = (float *) malloc (sizeof(float) * nt);    MTEST (*fitts);
01030   *errts    = (float *) malloc (sizeof(float) * nt);    MTEST (*errts);
01031 
01032 
01033   /*----- Check length of censor array -----*/
01034   if ((num_censor != 0) && (censor_length < nt))
01035     {
01036       DC_error ("Input censor time series file is too short");
01037       return (0);
01038     }
01039 
01040 
01041   /*----- Check validity of concatenated runs list -----*/
01042   for (ib = 0;  ib < num_blocks;  ib++)
01043     if ((block_list[ib] < 0) || (block_list[ib] >= nt))
01044       {
01045         DC_error ("Invalid concatenated runs list");
01046         return (0);
01047       }
01048   if (num_blocks > 1)
01049     for (ib = 1;  ib < num_blocks;  ib++)
01050       if (block_list[ib] <= block_list[ib-1])
01051         {
01052           DC_error ("Invalid concatenated runs list");
01053           return (0);
01054         }
01055     
01056 
01057   /*----- Create list of good (usable) data points -----*/
01058   good_list = (int *) malloc (sizeof(int) * nt);  MTEST (good_list);
01059   NFirst = plug_NFirst;
01060   if (NFirst < 0)
01061     for (is = 0;  is < num_stimts;  is++)
01062       if (NFirst < (max_lag[is]+nptr[is]-1)/nptr[is])  
01063         NFirst = (max_lag[is]+nptr[is]-1)/nptr[is];
01064   NLast = plug_NLast;   
01065 
01066   N = 0;
01067   ib = 0;
01068   for (it = block_list[0];  it < nt;  it++)
01069     {
01070       if (ib+1 < num_blocks)
01071         if (it >= block_list[ib+1])  ib++;
01072       
01073       irb = it - block_list[ib];
01074           
01075       if ((irb >= NFirst) && (irb <= NLast))
01076         if ((num_censor == 0) || (censor_array[it]))
01077           {
01078             good_list[N] = it;
01079             N++;
01080           }
01081     }
01082 
01083   if (N == 0)
01084     {
01085       DC_error ("No usable time points?");
01086       return (0);
01087     }
01088   if (N <= p)  
01089     {
01090       DC_error ("Insufficient time series data for deconvolution fit");
01091       return (0);
01092     }
01093 
01094   *NN = N;
01095 
01096 
01097   /*----- Perform initialization only if something has changed -----*/
01098   if (nt == prev_nt)
01099     {
01100       ok = 1;
01101     }
01102   else
01103     {
01104       /*----- Initialize the independent variable matrix -----*/
01105       demean_base = (plug_polort >= 0) ;
01106       ok = init_indep_var_matrix (p, qp, plug_polort, nt, N, good_list, 
01107                                   block_list, num_blocks, num_stimts, stimulus,
01108                                   stim_length, min_lag, max_lag, nptr, stim_base, &xdata);
01109 
01110       
01111       /*----- Initialization for the regression analysis -----*/
01112       if (ok)
01113         ok = init_regression_analysis (p, qp, num_stimts, stim_base, min_lag, 
01114                         max_lag, xdata, &x_full, &xtxinv_full, &xtxinvxt_full,
01115                         &x_base, &xtxinvxt_base, x_rdcd, xtxinvxt_rdcd);
01116 
01117 
01118       /*----- Initialization for the general linear test analysis -----*/
01119       if (ok)
01120         if (glt_num > 0)
01121           ok = init_glt_analysis (xtxinv_full, glt_num, glt_cmat, glt_amat, 
01122                                   cxtxinvct);
01123     }
01124       
01125   if (ok)
01126     {
01127       /*----- Extract Y-data for this voxel -----*/
01128       vector_create (N, &y);
01129       ts_array = vec;
01130       for (i = 0;  i < N;  i++)
01131         y.elts[i] = ts_array[good_list[i]];
01132       
01133       
01134       /*----- Perform the regression analysis for this voxel-----*/
01135       regression_analysis (N, p, q, num_stimts, min_lag, max_lag,
01136                            x_full, xtxinv_full, xtxinvxt_full, x_base,
01137                            xtxinvxt_base, x_rdcd, xtxinvxt_rdcd, 
01138                            y, rms_min, &mse, &coef, &scoef, &tcoef, 
01139                            fpart, rpart, &ffull, &rfull, &novar, 
01140                            *fitts, *errts);
01141       
01142           
01143       /*----- Perform the general linear tests for this voxel -----*/
01144       if (glt_num > 0)
01145         glt_analysis (N, p, x_full, y, mse*(N-p), coef, novar, cxtxinvct,
01146                       glt_num, glt_rows, glt_cmat, glt_amat, 
01147                       glt_coef, glt_tcoef, fglt, rglt);
01148       
01149      
01150       /*----- Save the fit parameters -----*/
01151       vector_to_array (coef, fit);
01152   
01153       
01154       /*----- Report results for this voxel -----*/
01155       printf ("\nResults for Voxel: \n");
01156       report_results (N, qp, q, p, plug_polort, block_list, num_blocks, 
01157                       num_stimts, stim_label, stim_base, min_lag, max_lag,
01158                       coef, tcoef, fpart, rpart, ffull, rfull, mse, 
01159                       glt_num, glt_label, glt_rows, glt_coef, 
01160                       glt_tcoef, fglt, rglt, label);
01161       printf ("%s \n", *label);
01162 
01163       prev_nt = nt;
01164     }
01165 
01166   else
01167     {
01168       vector_create (p, &coef);
01169       vector_to_array (coef, fit);
01170       strcpy (lbuf, "");
01171       *label = lbuf;
01172       prev_nt = 0;
01173     }
01174 
01175   
01176   /*----- Dispose of vectors -----*/
01177   vector_destroy (&y);
01178   vector_destroy (&tcoef);
01179   vector_destroy (&scoef);
01180   vector_destroy (&coef);
01181 
01182 
01183   if (ok)  return (1);
01184   else     return (0);
01185 }
01186 
01187 
01188 /*---------------------------------------------------------------------------*/
01189 /*
01190   Calculate the multiple linear regression least squares fit.
01191 */
01192 
01193 static void DC_Fit (int nt, double to, double dt, float * vec, char ** label)
01194 
01195 {
01196   int N;              /* first image from input 3d+time dataset to use */
01197   int n;                   /* time index */
01198   int ifit;                /* parameter index */
01199   int nfit;                /* number of fit parameters */
01200   float fit[MAX_XVARS];    /* fit parameters (regression coefficients) */
01201   int ok;                  /* Boolean for successful calculation */
01202   float * fitts = NULL;    /* full model fitted time series */
01203   float * errts = NULL;    /* full model residual error time series */
01204 
01205 
01206   /*----- Calculate the multiple linear regression -----*/
01207   ok = calculate_results (nt, dt, vec, &N, &nfit, fit, label,
01208                           &fitts, &errts);
01209 
01210 
01211   /*----- If unable to complete the calculation, return all zeros -----*/
01212   if (!ok)
01213     for (n = 0;  n < nt;  n++)  vec[n] = 0.0;
01214 
01215 
01216   /*----- Use full model fit to the time series data -----*/
01217   else
01218     {
01219       for (n = 0;  n < N;  n++)
01220         vec[good_list[n]] = fitts[n];
01221     }
01222 
01223 
01224   /*----- Deallocate memory -----*/
01225   free (fitts);   fitts = NULL;
01226   free (errts);   errts = NULL;
01227 
01228   return;
01229 }
01230 
01231 
01232 /*---------------------------------------------------------------------------*/
01233 /*
01234   Calculate the multiple linear regression residuals.
01235 */
01236 
01237 static void DC_Err (int nt, double to, double dt, float * vec, char ** label)
01238 
01239 {
01240   int N;              /* first image from input 3d+time dataset to use */
01241   float val;               /* residual value at a time point */ 
01242   int n;                   /* time index */
01243   int ifit;                /* parameter index */
01244   int nfit;                /* number of fit parameters */
01245   float fit[MAX_XVARS];    /* fit parameters (regression coefficients) */
01246   int ok;                  /* Boolean for successful calculation */
01247   float * fitts = NULL;    /* full model fitted time series */
01248   float * errts = NULL;    /* full model residual error time series */
01249 
01250 
01251   /*----- Calculate the multiple linear regression -----*/
01252   ok = calculate_results (nt, dt, vec, &N, &nfit, fit, label,
01253                           &fitts, &errts);
01254 
01255 
01256   /*----- If unable to complete the calculation, return all zeros -----*/
01257   for (n = 0;  n < nt;  n++)  vec[n] = 0.0;
01258 
01259 
01260   /*----- Use residuals from full model fit to time series data -----*/
01261   if (ok)
01262     {
01263       for (n = 0;  n < N;  n++)
01264         vec[good_list[n]] = errts[n];
01265     }
01266 
01267 
01268   /*----- Deallocate memory -----*/
01269   free (fitts);   fitts = NULL;
01270   free (errts);   errts = NULL;
01271 
01272   return;
01273 }
01274 
01275 
01276 /*---------------------------------------------------------------------------*/
01277 /*
01278   Estimate the system impulse response function.
01279 */
01280  
01281 static void DC_IRF (int nt, double to, double dt, float * vec, char ** label)
01282 
01283 {
01284   int N;              /* first image from input 3d+time dataset to use */
01285   int nfit;                /* number of fit parameters */
01286   float fit[MAX_XVARS];    /* fit parameters (regression coefficients) */
01287   int np;                  /* length of estimated impulse reponse function */
01288   int ip;                  /* impulse response function parameter index */
01289   int q;                   /* number of parameters in the baseline model */
01290   int it;                  /* array index */
01291   int ntdnp;               /* number of array points per IRF parameter */  
01292   int ok;                  /* Boolean for successful calculation */
01293   float * fitts = NULL;    /* full model fitted time series */
01294   float * errts = NULL;    /* full model residual error time series */
01295 
01296 
01297   /*----- Calculate the multiple linear regression -----*/
01298   ok = calculate_results (nt, dt, vec, &N, &nfit, fit, label,
01299                           &fitts, &errts);
01300 
01301 
01302   /*----- If unable to complete the calculation, return all zeros -----*/
01303   if (!ok || (num_stimts < 1))
01304     for (it = 0;  it < nt;  it++)  vec[it] = 0.0;
01305 
01306 
01307   /*----- Plot the system impulse response function -----*/
01308   else
01309     {
01310       if ((num_stimts == 1) || (plug_IRF < 0) || (plug_IRF >= num_stimts))
01311         plug_IRF = 0;
01312       
01313       np = max_lag[plug_IRF] - min_lag[plug_IRF] + 1;
01314       
01315       q = plug_qp;
01316       for (ip = 0;  ip < plug_IRF;  ip++)
01317         q += max_lag[ip] - min_lag[ip] + 1;
01318 
01319       if (np == 1)
01320         {
01321           for (it = 0;  it < nt;  it++)
01322             vec[it] = fit[q];
01323         }
01324       else
01325         {
01326           float r;
01327 
01328           ntdnp = nt / (np-1);
01329       
01330           vec[0] = fit[q];
01331           for (it = 0;  it < nt;  it++)
01332             {
01333               ip = it / ntdnp + 1;
01334               if (ip < np)
01335                 {
01336                   r = (float) it / (float) ntdnp - (ip-1);
01337                   vec[it] = r * fit[q+ip] + (1.0-r) * fit[q+ip-1]; 
01338                 }
01339               else
01340                 vec[it] = fit[q+np-1];
01341             }
01342         }
01343     }
01344 
01345 
01346   /*----- Deallocate memory -----*/
01347   free (fitts);   fitts = NULL;
01348   free (errts);   errts = NULL;
01349   
01350   return;
01351 }
01352 
01353 
01354 /*---------------------------------------------------------------------------*/
01355 
01356 
01357 
01358 
01359 
01360 
01361 
01362 
01363 
01364 
01365 
01366 
01367 
01368