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 }

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 }

00163 {
00164   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00165 
00166 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

 {"None", "Const", "Linear", 
                            "Quadrtc", "Cubic", "Quartic", "Quintic" }