do{ fprintf(stderr,"** ERROR: Can't read GLT matrix from file %s\n",fname);  \
     exit(1) ; } while(0)

03067 {
03068   int ip;                  /* parameter index */
03069   int qp;                  /* number of polynomial trend baseline parameters */
03070   int q;                   /* number of baseline model parameters */
03071   int p;                   /* number of full model parameters */
03072   int nxyz;                /* total number of voxels */
03073   int ixyz;                /* voxel index */
03074   int is;                  /* stimulus index */
03075   int num_stimts;          /* number of stimulus time series */
03076   int num_glt;             /* number of general linear tests */
03077   int iglt;                /* general linear test index */
03078   int nlc;                 /* number of linear combinations in a GLT */
03079   int ilc;                 /* linear combination index */
03080   int it;                  /* time point index */
03081   int nt;                  /* number of images in input 3d+time dataset */
03082   int * min_lag;           /* minimum time delay for impulse response */
03083   int * max_lag;           /* maximum time delay for impulse response */
03084 
03085   int fout;             /* flag to output F-statistics */
03086   int rout;             /* flag to output R^2 statistics */
03087   int tout;             /* flag to output t-statistics */
03088   int vout;             /* flag to output variance map */
03089   int bout;             /* flag to output baseline coefficients */
03090   int cout;             /* flag to output fit coefficients */
03091   int ibot,itop ;
03092 
03093 
03094 ENTRY("allocate_memory") ;
03095 
03096   /*----- Initialize local variables -----*/
03097   nxyz = option_data->nxyz;
03098   num_stimts = option_data->num_stimts;
03099   num_glt = option_data->num_glt;
03100   qp = option_data->qp;
03101   q  = option_data->q;
03102   p  = option_data->p;
03103   nt = option_data->nt;
03104   min_lag = option_data->stim_minlag;
03105   max_lag = option_data->stim_maxlag;
03106 
03107   fout  = option_data->fout;
03108   rout  = option_data->rout;
03109   tout  = option_data->tout;
03110   vout  = option_data->vout;
03111   bout  = 1 - (option_data->nobout || option_data->nocout);
03112   cout  = 1 - option_data->nocout;
03113 
03114   if( CoefFilename != NULL ){ bout = cout = 1 ; }  /* 26 Jul 2004 */
03115 
03116   /*----- Allocate memory space for volume data -----*/
03117   *coef_vol  = (float **) malloc (sizeof(float *) * p);   MTEST(*coef_vol);
03118   *scoef_vol = (float **) malloc (sizeof(float *) * p);   MTEST(*scoef_vol);
03119   *tcoef_vol = (float **) malloc (sizeof(float *) * p);   MTEST(*tcoef_vol);
03120   for (ip = 0;  ip < p;  ip++)
03121     {
03122       (*coef_vol)[ip] = NULL;
03123       (*scoef_vol)[ip] = NULL;
03124       (*tcoef_vol)[ip] = NULL;
03125     }
03126 
03127   if (bout)
03128     for (ip = 0;  ip < qp;  ip++)
03129       {
03130       zero_fill_volume (&((*coef_vol)[ip]),  nxyz);
03131       if (tout) zero_fill_volume (&((*tcoef_vol)[ip]),  nxyz);
03132       }
03133 
03134   ip = qp - 1;
03135   for (is = 0;  is < num_stimts;  is++)
03136     {
03137       if( basis_stim[is] != NULL ){ ibot=0 ; itop=basis_stim[is]->nfunc-1 ; }
03138       else                        { ibot=min_lag[is] ; itop=max_lag[is] ;   }
03139       for (it = ibot;  it <= itop;  it++)
03140       {
03141         ip++;
03142         if (option_data->stim_base[is])
03143           {
03144             if (bout || (option_data->iresp_filename[is] != NULL))
03145             zero_fill_volume (&((*coef_vol)[ip]),  nxyz);
03146             if (bout && tout)
03147             zero_fill_volume (&((*tcoef_vol)[ip]), nxyz);
03148           }
03149         else
03150           {
03151             if (cout || (option_data->iresp_filename[is] != NULL))
03152             zero_fill_volume (&((*coef_vol)[ip]),  nxyz);
03153             if (cout && tout)
03154             zero_fill_volume (&((*tcoef_vol)[ip]), nxyz);
03155           }
03156         if (option_data->sresp_filename[is] != NULL)
03157           zero_fill_volume (&((*scoef_vol)[ip]), nxyz);
03158       }
03159     }
03160 
03161 
03162   if (fout)
03163     {
03164       *fpart_vol = (float **) malloc (sizeof(float *) * num_stimts);
03165       MTEST(*fpart_vol);
03166       for (is = 0;  is < num_stimts;  is++)
03167       if ((! option_data->stim_base[is]) || (! option_data->nobout))
03168         zero_fill_volume (&((*fpart_vol)[is]), nxyz);
03169       else
03170         (*fpart_vol)[is] = NULL;
03171     }
03172 
03173 
03174   if (rout)
03175     {
03176       *rpart_vol = (float **) malloc (sizeof(float *) * num_stimts);
03177       MTEST(*rpart_vol);
03178       for (is = 0;  is < num_stimts;  is++)
03179       if ((! option_data->stim_base[is]) || (! option_data->nobout))
03180         zero_fill_volume (&((*rpart_vol)[is]), nxyz);
03181       else
03182         (*rpart_vol)[is] = NULL;
03183     }
03184 
03185 
03186   if (vout || basis_need_mse)  zero_fill_volume (&(*mse_vol),   nxyz);
03187   if (fout)  zero_fill_volume (&(*ffull_vol), nxyz);
03188   if (rout)  zero_fill_volume (&(*rfull_vol), nxyz);
03189 
03190 
03191   /*----- Allocate memory space for GLT volume data -----*/
03192   if (num_glt > 0)
03193     {
03194       *glt_coef_vol  = (float ***) malloc (sizeof(float **) * num_glt);
03195        MTEST(*glt_coef_vol);
03196 
03197        if (tout)
03198        {
03199          *glt_tcoef_vol  = (float ***) malloc (sizeof(float **) * num_glt);
03200          MTEST(*glt_tcoef_vol);
03201        }
03202 
03203        if (fout)
03204        {
03205          *glt_fstat_vol = (float **)  malloc (sizeof(float *)  * num_glt);
03206          MTEST(*glt_fstat_vol);
03207        }
03208 
03209        if (rout)
03210        {
03211          *glt_rstat_vol = (float **)  malloc (sizeof(float *)  * num_glt);
03212          MTEST(*glt_rstat_vol);
03213        }
03214 
03215       for (iglt = 0;  iglt < num_glt;  iglt++)
03216       {
03217         nlc = option_data->glt_rows[iglt];
03218 
03219         (*glt_coef_vol)[iglt] = (float **) malloc (sizeof(float *) * nlc);
03220         MTEST((*glt_coef_vol)[iglt]);
03221 
03222         if (tout)
03223           {
03224             (*glt_tcoef_vol)[iglt]
03225             = (float **) malloc (sizeof(float *) * nlc);
03226             MTEST((*glt_tcoef_vol)[iglt]);
03227           }
03228 
03229         for (ilc = 0;  ilc < nlc;  ilc++)
03230           {
03231             zero_fill_volume (&((*glt_coef_vol)[iglt][ilc]),  nxyz);
03232             if (tout)
03233             zero_fill_volume (&((*glt_tcoef_vol)[iglt][ilc]),  nxyz);
03234           }
03235 
03236         if (fout)  zero_fill_volume (&((*glt_fstat_vol)[iglt]),  nxyz);
03237         if (rout)  zero_fill_volume (&((*glt_rstat_vol)[iglt]),  nxyz);
03238       }
03239     }
03240 
03241 
03242   /*----- Allocate memory for fitted time series -----*/
03243   if (option_data->fitts_filename != NULL)
03244     {
03245       *fitts_vol = (float **) malloc (sizeof(float **) * nt);
03246       MTEST (*fitts_vol);
03247       for (it = 0;  it < nt;  it++)
03248       {
03249         zero_fill_volume (&((*fitts_vol)[it]),  nxyz);
03250       }
03251     }
03252 
03253   /*----- Allocate memory for residual errors -----*/
03254   if (option_data->errts_filename != NULL)
03255     {
03256       *errts_vol = (float **) malloc (sizeof(float **) * nt);
03257       MTEST (*errts_vol);
03258       for (it = 0;  it < nt;  it++)
03259       {
03260         zero_fill_volume (&((*errts_vol)[it]),  nxyz);
03261       }
03262     }
03263 
03264 #ifdef PROC_MAX
03265   if( proc_numjob > 1 ) proc_finalize_shm_volumes() ;  /* RWCox */
03266 #endif
03267 
03268   EXRETURN ;
03269 }

04684 {
04685   const float EPSILON = 1.0e-10;
04686   float factor;             /* factor is new scale factor for sub-brick #ib */
04687   short *sbr ;
04688 
04689 ENTRY("attach_sub_brick") ;
04690 
04691 
04692   /*----- allocate memory for output sub-brick -----*/
04693   sbr = (short *) malloc (sizeof(short) * nxyz);
04694   MTEST (sbr);
04695   factor = EDIT_coerce_autoscale_new(nxyz, MRI_float,volume, MRI_short,sbr);
04696 
04697   if (factor < EPSILON)  factor = 0.0;
04698   else                   factor = 1.0 / factor;
04699 
04700   /*----- edit the sub-brick -----*/
04701   EDIT_BRICK_LABEL (new_dset, ibrick, brick_label);
04702   EDIT_BRICK_FACTOR(new_dset, ibrick, factor);
04703 
04704   if (brick_type == FUNC_TT_TYPE)
04705     EDIT_BRICK_TO_FITT (new_dset, ibrick, dof);
04706   else if (brick_type == FUNC_FT_TYPE)
04707     EDIT_BRICK_TO_FIFT (new_dset, ibrick, ndof, ddof);
04708 
04709 
04710   /*----- attach sbr to be sub-brick #ibrick -----*/
04711   EDIT_substitute_brick (new_dset, ibrick, MRI_short, sbr);
04712   if( bar != NULL ) bar[ibrick] = sbr ;
04713 
04714   EXRETURN ;
04715 }

07789 {
07790    register int jj ;
07791    register double sum ;
07792 
07793    sum = 0.0 ;
07794    for( jj=0 ; jj < nParam ; jj++ )
07795      if( Xcol_inbase[jj] == 2 ) sum += coef.elts[jj] * Xcol_mean[jj] ;
07796 
07797    return (float)sum ;
07798 }

07057 {
07058    float w , tp , pp ;
07059 
07060    w = basis_block_hrf4(t,T) ;
07061    if( w > 0.0f && peak > 0.0f ){
07062      tp = TPEAK4(T) ;  pp = basis_block_hrf4(tp,T) ;
07063      if( pp > 0.0f ) w *= peak / pp ;
07064    }
07065    return w ;
07066 }

07121 {
07122    float w , tp , pp ;
07123 
07124    w = basis_block_hrf5(t,T) ;
07125    if( w > 0.0f && peak > 0.0f ){
07126      tp = TPEAK5(T) ;  pp = basis_block_hrf5(tp,T) ;
07127      if( pp > 0.0f ) w *= peak / pp ;
07128    }
07129    return w ;
07130 }

07025 {
07026   register double t26, t2, t4, t1, t42, t12, t34, t35, t16, t46, t,L ;
07027   double w ;
07028 
07029   if( tt <= 0.0f || tt >= (TT+15.0f) ) return 0.0f ;
07030 
07031   t = tt ; L = TT ; t4 = exp(0.4e1 - t);
07032   if( t < L ){ L = t ; t16 = 54.5982 ; }
07033   else       { t16 = exp(4.0-t+L) ;    }
07034 
07035   t1 = t * t;
07036   t2 = t1 * t1;
07037   t4 = exp(4.0 - t);
07038   t12 = t1 * t;
07039   t26 = t16 * L;
07040   t34 = L * L;
07041   t35 = t16 * t34;
07042   t42 = t16 * t34 * L;
07043   t46 = t34 * t34;
07044 
07045   w = -t2 * t4 / 0.256e3 - 0.3e1 / 0.32e2 * t4 - 0.3e1 / 0.32e2 * t4 * t
07046       - 0.3e1 / 0.64e2 * t4 * t1 - t4 * t12 / 0.64e2 + t16 * t2 / 0.256e3
07047       + 0.3e1 / 0.32e2 * t16 + 0.3e1 / 0.32e2 * t16 * t
07048       + 0.3e1 / 0.64e2 * t1 * t16 + t16 * t12 / 0.64e2 - 0.3e1 / 0.32e2 * t26
07049       - 0.3e1 / 0.32e2 * t26 * t - 0.3e1 / 0.64e2 * t1 * t26
07050       - t26 * t12 / 0.64e2 + 0.3e1 / 0.64e2 * t35 + 0.3e1 / 0.64e2 * t35 * t
07051       + 0.3e1 / 0.128e3 * t1 * t35 - t42 / 0.64e2 - t42 * t / 0.64e2
07052       + t16 * t46 / 0.256e3 ;
07053   return (float)w ;
07054 }

07078 {
07079    register double t , T ;
07080    register double t2,t3,t4,t5,t6,t7,t9,t10,t11,t14,t20,t25,t28,t37,t57 ;
07081    double w ;
07082 
07083    if( tt <= 0.0f || tt >= (TT+15.0f) ) return 0.0f ;
07084 
07085    t = tt ; T = TT ;
07086 
07087 #if 1
07088    t2 = exp(-t) ;
07089    if( t <= T ){ t3 = t ; t4 = 1.0/t2 ; }
07090    else        { t3 = T ; t4 = exp(T)  ; }
07091    t2 *= 148.413 ;    /* 148.413 = exp(5) */
07092 #else
07093    t2 = exp(0.5e1 - t);
07094    t3 = (t <= T ? t : T);
07095    t4 = exp(t3);
07096 #endif
07097    t5 = t * t;
07098    t6 = t5 * t5;
07099    t7 = t6 * t;
07100    t9 = t3 * t3;
07101    t10 = t9 * t9;
07102    t11 = t4 * t10;
07103    t14 = t4 * t9 * t3;
07104    t20 = t4 * t3;
07105    t25 = t4 * t9;
07106    t28 = t5 * t;
07107    t37 = -0.120e3 + t4 * t7 + 0.5e1 * t11 - 0.20e2 * t14 - t4 * t10 * t3
07108          - 0.10e2 * t14 * t5 - 0.120e3 * t20 * t - 0.20e2 * t14 * t
07109          + 0.30e2 * t25 * t5 + 0.10e2 * t25 * t28 + 0.5e1 * t11 * t
07110          + 0.20e2 * t4 * t28 + 0.60e2 * t25 * t;
07111    t57 = -0.5e1 * t20 * t6 - 0.20e2 * t20 * t28 - 0.60e2 * t20 * t5
07112          - 0.5e1 * t6 - 0.20e2 * t28 + 0.120e3 * t4 - 0.120e3 * t
07113          - 0.120e3 * t20 + 0.60e2 * t25 - t7 - 0.60e2 * t5 + 0.120e3 * t4 * t
07114          + 0.60e2 * t4 * t5 + 0.5e1 * t4 * t6;
07115    w = t2 * (t37 + t57) / 0.3125e4;
07116 
07117    return (float)w ;
07118 }

06960 {
06961    if( x < bot || x > top ) return 0.0f ;
06962    return (float)cos(freq*(x-bot)) ;
06963 }

06918 {
06919    float sum=0.0 ;
06920    int ii ;
06921 
06922    if( x >= be->tbot && x <= be->ttop ){
06923      for( ii=0 ; ii < be->nfunc ; ii++ )
06924        sum += wt[ii] * basis_funceval( be->bfunc[ii] , x ) ;
06925    }
06926 
06927    return sum ;
06928 }

07181 {
07182    PARSER_code *pc = (PARSER_code *)q ;
07183    double atoz[26] , val ;
07184 
07185    if( x < bot || x > top ) return 0.0f ;
07186    atoz[ITT] = x ;                            /* t = true time from stim */
07187    atoz[IXX] = (x-bot)*dtinv ;                /* x = scaled to [0,1] */
07188    atoz[IZZ] = 2.0*atoz[IXX] - 1.0 ;          /* z = scaled to [-1,1] */
07189    val = PARSER_evaluate_one( pc , atoz ) ;
07190    return (float)val ;
07191 }

06984 {
06985    if( x <= 0.0f || x > top ) return 0.0f ;
06986    return (float)(pow(x/(b*c),b)*exp(b-x/c)) ;
06987 }

07139 {
07140    float xq ;
07141 
07142    x = 2.0f*(x-bot)/(top-bot) - 1.0f ;  /* now in range -1..1 */
07143 
07144    if( x < -1.0f || x > 1.0f ) return 0.0f ;
07145 
07146    xq = x*x ;
07147 
07148    switch( (int)n ){
07149     case 0: return 1.0f ;
07150     case 1: return x ;
07151     case 2: return (3.0f*xq-1.0f)/2.0f ;
07152     case 3: return (5.0f*xq-3.0f)*x/2.0f ;
07153     case 4: return ((35.0f*xq-30.0f)*xq+3.0f)/8.0f ;
07154     case 5: return ((63.0f*xq-70.0f)*xq+15.0f)*x/8.0f ;
07155     case 6: return (((231.0f*xq-315.0f)*xq+105.0f)*xq-5.0f)/16.0f ;
07156     case 7: return (((429.0f*xq-693.0f)*xq+315.0f)*xq-35.0f)*x/16.0f ;
07157 
07158     case 8: return ((((6435.0f*xq-12012.0f)*xq+6930.0f)*xq-1260.0f)*xq+35.0f)
07159                   /128.0f;
07160 
07161     case 9: return ((((12155.0f*xq-25740.0f)*xq+18018.0f)*xq-4620.0f)*xq+315.0f)
07162                   *x/128.0f ;
07163    }
07164 
07165    return 0.0f ;   /* should never be reached */
07166 }

06948 {
06949    if( x < bot || x > top ) return 0.0f ;
06950    return 1.0f ;
06951 }

07198 {
07199    basis_expansion *be ;
07200    char *cpt , *scp ;
07201    float bot=0.0f, top=0.0f ;
07202    int nn , nord=0 ;
07203 
07204    if( sym == NULL ) return NULL ;
07205 
07206    scp = strdup(sym) ;                        /* duplicate, for editing */
07207    cpt = strchr(scp,'(') ;                    /* find opening '(' */
07208    if( cpt != NULL ){ *cpt = '\0' ; cpt++ ; } /* cut string there */
07209 
07210    be = (basis_expansion *)malloc(sizeof(basis_expansion)) ;
07211    be->name = NULL ;   /* will be fixed later */
07212 
07213    /*--- GAM(b,c) ---*/
07214 
07215    if( strcmp(scp,"GAM") == 0 ){
07216 
07217      be->nfunc = 1 ;
07218      be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07219      be->bfunc[0].f = basis_gam ;
07220      if( cpt == NULL ){
07221        be->bfunc[0].a = 8.6f ;     /* Mark Cohen's parameters */
07222        be->bfunc[0].b = 0.547f ;   /* t_peak=4.7 FWHM=3.7 */
07223        be->bfunc[0].c = 11.1f ;    /* return to zero-ish */
07224      } else {
07225        sscanf(cpt,"%f,%f",&bot,&top) ;
07226        if( bot <= 0.0f || top <= 0.0f ){
07227          fprintf(stderr,"** ERROR: 'GAM(%s' is illegal\n",cpt) ;
07228          fprintf(stderr,
07229           "          Correct format: 'GAM(b,c)' with b > 0 and c > 0.\n");
07230          free((void *)be->bfunc); free((void *)be); free(scp); return NULL;
07231        }
07232        be->bfunc[0].a = bot ;    /* t_peak = bot*top */
07233        be->bfunc[0].b = top ;    /* FWHM   = 2.3*sqrt(bot)*top */
07234        be->bfunc[0].c = bot*top + 4.0f*sqrt(bot)*top ;  /* long enough */
07235      }
07236      be->tbot = 0.0f ; be->ttop = be->bfunc[0].c ;
07237 
07238    /*--- TENT(bot,top,order) ---*/
07239 
07240    } else if( strcmp(scp,"TENT") == 0 ){
07241      float dx ;
07242 
07243      if( cpt == NULL ){
07244        fprintf(stderr,"** ERROR: 'TENT' by itself is illegal\n") ;
07245        fprintf(stderr,
07246         "          Correct format: 'TENT(bot,top,n)' with bot < top and n > 0.\n") ;
07247        free((void *)be); free(scp); return NULL;
07248      }
07249      sscanf(cpt,"%f,%f,%d",&bot,&top,&nord) ;
07250      if( bot >= top || nord < 2 ){
07251        fprintf(stderr,"** ERROR: 'TENT(%s' is illegal\n",cpt) ;
07252        fprintf(stderr,
07253         "          Correct format: 'TENT(bot,top,n)' with bot < top and n > 1.\n") ;
07254        free((void *)be); free(scp); return NULL;
07255      }
07256      be->nfunc = nord ;
07257      be->tbot  = bot  ; be->ttop = top ;
07258      be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07259      dx        = (top-bot) / (nord-1) ;
07260 
07261      be->bfunc[0].f = basis_tent ;
07262      be->bfunc[0].a = bot-0.001*dx ;
07263      be->bfunc[0].b = bot ;
07264      be->bfunc[0].c = bot+dx ;
07265      for( nn=1 ; nn < nord-1 ; nn++ ){
07266        be->bfunc[nn].f = basis_tent ;
07267        be->bfunc[nn].a = bot + (nn-1)*dx ;
07268        be->bfunc[nn].b = bot +  nn   *dx ;
07269        be->bfunc[nn].c = bot + (nn+1)*dx ;
07270      }
07271      be->bfunc[nord-1].f = basis_tent ;
07272      be->bfunc[nord-1].a = bot + (nord-2)*dx ;
07273      be->bfunc[nord-1].b = top ;
07274      be->bfunc[nord-1].c = top + 0.001*dx ;
07275 
07276    /*--- TRIG(bot,top,order) ---*/
07277 
07278    } else if( strcmp(scp,"TRIG") == 0 ){
07279 
07280      if( cpt == NULL ){
07281        fprintf(stderr,"** ERROR: 'TRIG' by itself is illegal\n") ;
07282        fprintf(stderr,
07283         "          Correct format: 'TRIG(bot,top,n)' with bot < top and n > 2.\n") ;
07284        free((void *)be); free(scp); return NULL;
07285      }
07286      sscanf(cpt,"%f,%f,%d",&bot,&top,&nord) ;
07287      if( bot >= top || nord < 3 ){
07288        fprintf(stderr,"** ERROR: 'TRIG(%s' is illegal\n",cpt) ;
07289        fprintf(stderr,
07290         "          Correct format: 'TRIG(bot,top,n)' with bot < top and n > 2.\n") ;
07291        free((void *)be); free(scp); return NULL;
07292      }
07293      be->nfunc = nord ;
07294      be->tbot  = bot  ; be->ttop = top ;
07295      be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07296 
07297      be->bfunc[0].f = basis_one ;
07298      be->bfunc[0].a = bot ;
07299      be->bfunc[0].b = top ;
07300      be->bfunc[0].c = 0.0f ;
07301      for( nn=1 ; nn < nord ; nn++ ){
07302        be->bfunc[nn].f = basis_cos ;
07303        be->bfunc[nn].a = bot ;
07304        be->bfunc[nn].b = top ;
07305        be->bfunc[nn].c = (2.0*PI)*((nn+1)/2)/(top-bot) ;
07306        nn++ ; if( nn >= nord ) break ;
07307        be->bfunc[nn].f = basis_sin ;
07308        be->bfunc[nn].a = bot ;
07309        be->bfunc[nn].b = top ;
07310        be->bfunc[nn].c = (2.0*PI)*((nn+1)/2)/(top-bot) ;
07311      }
07312 
07313    /*--- SIN(bot,top,order) ---*/
07314 
07315    } else if( strcmp(scp,"SIN") == 0 ){
07316 
07317      if( cpt == NULL ){
07318        fprintf(stderr,"** ERROR: 'SIN' by itself is illegal\n") ;
07319        fprintf(stderr,
07320         "          Correct format: 'SIN(bot,top,n)' with bot < top and n > 0.\n") ;
07321        free((void *)be); free(scp); return NULL;
07322      }
07323      sscanf(cpt,"%f,%f,%d",&bot,&top,&nord) ;
07324      if( bot >= top || nord < 1 ){
07325        fprintf(stderr,"** ERROR: 'SIN(%s' is illegal\n",cpt) ;
07326        fprintf(stderr,
07327         "          Correct format: 'SIN(bot,top,n)' with bot < top and n > 0.\n") ;
07328        free((void *)be); free(scp); return NULL;
07329      }
07330      be->nfunc = nord ;
07331      be->tbot  = bot  ; be->ttop = top ;
07332      be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07333 
07334      for( nn=0 ; nn < nord ; nn++ ){
07335        be->bfunc[nn].f = basis_sin ;
07336        be->bfunc[nn].a = bot ;
07337        be->bfunc[nn].b = top ;
07338        be->bfunc[nn].c = PI*(nn+1)/(top-bot) ;
07339      }
07340 
07341    /*--- POLY(bot,top,order) ---*/
07342 
07343    } else if( strcmp(scp,"POLY") == 0 ){
07344 
07345      if( cpt == NULL ){
07346        fprintf(stderr,"** ERROR: 'POLY' by itself is illegal\n") ;
07347        fprintf(stderr,
07348         "          Correct format: 'POLY(bot,top,n)' with bot < top and n > 0.\n") ;
07349        free((void *)be); free(scp); return NULL;
07350      }
07351      sscanf(cpt,"%f,%f,%d",&bot,&top,&nord) ;
07352      if( bot >= top || nord < 1 || nord > POLY_MAX ){
07353        fprintf(stderr,"** ERROR: 'POLY(%s' is illegal\n",cpt) ;
07354        fprintf(stderr,
07355         "          Correct format: 'POLY(bot,top,n)' with bot < top and n > 0.\n") ;
07356        free((void *)be); free(scp); return NULL;
07357      }
07358      be->nfunc = nord ;
07359      be->tbot  = bot  ; be->ttop = top ;
07360      be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07361      be->bfunc[0].f = basis_one ;
07362      be->bfunc[0].a = bot ;
07363      be->bfunc[0].b = top ;
07364      be->bfunc[0].c = 0.0f ;
07365      for( nn=1 ; nn < nord ; nn++ ){
07366        be->bfunc[nn].f = basis_legendre ;
07367        be->bfunc[nn].a = bot ;
07368        be->bfunc[nn].b = top ;
07369        be->bfunc[nn].c = (float)nn ;
07370      }
07371 
07372    /*--- SPMG ---*/
07373 
07374    } else if( strncmp(scp,"SPMG",4) == 0 ){
07375 
07376      be->nfunc = 2 ;
07377      be->tbot  = 0.0f ; be->ttop = 25.0f ;
07378      be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07379      be->bfunc[0].f = basis_spmg1 ;
07380      be->bfunc[0].a = 0.0f ;
07381      be->bfunc[0].b = 0.0f ;
07382      be->bfunc[0].c = 0.0f ;
07383      be->bfunc[1].f = basis_spmg2 ;
07384      be->bfunc[1].a = 0.0f ;
07385      be->bfunc[1].b = 0.0f ;
07386      be->bfunc[1].c = 0.0f ;
07387 
07388    /*--- BLOCKn(duration,peak) for n=4 or 5 ---*/
07389 
07390    } else if( strncmp(scp,"BLOCK",5) == 0 || strncmp(scp,"IGFUN",5) == 0 ){
07391      int nb=4 ;
07392 
07393      if( scp[5] != '\0' ){
07394        nb = strtol( scp+5 , NULL , 10 ) ;
07395        if( nb != 4 && nb != 5 ){
07396          fprintf(stderr,
07397          "** ERROR: '%s' has illegal power: only 4 or 5 allowed\n",scp) ;
07398          free((void *)be); free(scp); return NULL;
07399        }
07400      }
07401 
07402      if( cpt == NULL ){
07403        fprintf(stderr,"** ERROR: '%s' by itself is illegal\n",scp) ;
07404        fprintf(stderr,
07405         "          Correct format: 'BLOCKn(dur)' with dur > 0, n=4 or 5\n"
07406         "                     *OR* 'BLOCKn(dur,peak)' with peak > 0 too.\n") ;
07407        free((void *)be); free(scp); return NULL;
07408      }
07409      sscanf(cpt,"%f,%f",&top,&bot) ;
07410      if( top <= 0.0f ){
07411        fprintf(stderr,"** ERROR: '%s(%s' is illegal\n",scp,cpt) ;
07412        fprintf(stderr,
07413         "          Correct format: 'BLOCKn(dur)' with dur > 0, n=4 or 4\n"
07414         "                      or: 'BLOCKn(dur,peak)' with peak > 0 too.\n") ;
07415        free((void *)be); free(scp); return NULL;
07416      }
07417      if( bot < 0.0f ) bot = 0.0f ;
07418 
07419      be->nfunc = 1 ;
07420      be->tbot  = 0.0f ; be->ttop = top+15.0f ;
07421      be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07422      be->bfunc[0].f = (nb==5) ? basis_block5 : basis_block4 ;
07423      be->bfunc[0].a = top ;
07424      be->bfunc[0].b = bot ;
07425      be->bfunc[0].c = 0.0f ;
07426 
07427    /*--- EXPR(bot,top) exp1 exp2 ... ---*/
07428 
07429    } else if( strcmp(scp,"EXPR") == 0 ){   /* 28 Aug 2004 */
07430      char *ept ;
07431      NI_str_array *sar ;
07432      int nexpr , ie ;
07433      PARSER_code *pc ;
07434 
07435      if( cpt == NULL ){
07436        fprintf(stderr,"** ERROR: 'EXPR' by itself is illegal\n") ;
07437        free((void *)be); free(scp); return NULL;
07438      }
07439      sscanf(cpt,"%f,%f",&bot,&top) ;
07440      if( top <= bot ){
07441        fprintf(stderr,"** ERROR: 'EXPR(%f,%f)' has illegal time range\n",bot,top) ;
07442        free((void *)be); free(scp); return NULL;
07443      }
07444      ept = strchr( cpt , ')' ) ;
07445      if( ept == NULL ){
07446        fprintf(stderr,"** ERROR: 'EXPR(%f,%f)' has no expressions!?\n",bot,top);
07447        free((void *)be); free(scp); return NULL;
07448      }
07449      sar = NI_decode_string_list( ept+1 , "~" ) ;
07450      if( sar == NULL || sar->num == 0 ){
07451        fprintf(stderr,"** ERROR: 'EXPR(%f,%f)' has no expressions!?\n",bot,top);
07452        free((void *)be); free(scp); return NULL;
07453      }
07454      be->nfunc = nexpr = sar->num ;
07455      be->tbot  = bot  ; be->ttop = top ;
07456      be->bfunc = (basis_func *)calloc(sizeof(basis_func),be->nfunc) ;
07457      PARSER_set_printout(1) ;
07458      for( ie=0 ; ie < nexpr ; ie++ ){
07459        pc = PARSER_generate_code( sar->str[ie] ) ;
07460        if( pc == NULL ){
07461          fprintf(stderr,"** ERROR: unparsable EXPRession '%s'\n",sar->str[ie]) ;
07462          free((void *)be); free(scp); return NULL;
07463        }
07464        if( strcmp(sar->str[ie],"1") != 0 &&    /* must either be "1" */
07465            !PARSER_has_symbol("t",pc)    &&    /* or contain symbol  */
07466            !PARSER_has_symbol("x",pc)    &&    /* 't', 'x', or 'z'   */
07467            !PARSER_has_symbol("z",pc)      ){
07468          fprintf(stderr,
07469            "** ERROR: illegal EXPRession '%s' lacks symbol t, x, or z\n",
07470            sar->str[ie]) ;
07471          free((void *)be); free(scp); return NULL;
07472        }
07473        be->bfunc[ie].f = basis_expr ;
07474        be->bfunc[ie].a = bot ;
07475        be->bfunc[ie].b = top ;
07476        be->bfunc[ie].c = 1.0f/(top-bot) ;  /* for ease of scaling */
07477        be->bfunc[ie].q = (void *)pc ;
07478      }
07479      PARSER_set_printout(0) ;
07480 
07481      NI_delete_str_array(sar) ;
07482 
07483    /*--- NO MORE BASIS FUNCTION CHOICES ---*/
07484 
07485    } else {
07486      fprintf(stderr,"** ERROR: '%s' is unknown response function type\n",scp) ;
07487      free((void *)be); free(scp); return NULL;
07488    }
07489 
07490    /* 28 Apr 2005: set scaling factors */
07491 
07492    for( nn=0 ; nn < be->nfunc ; nn++ )  /* initialize factors to 1.0 */
07493      be->bfunc[nn].ffac = 1.0 ;
07494 
07495    /* for each basis function, find its peak value, then
07496       set ffac so that the peak value becomes basis_normall */
07497 
07498 #undef  BNSUB
07499 #define BNSUB 999
07500    if( basis_normall > 0.0f ){
07501      int jj ; float dt , ftop , val ;
07502      bot = be->tbot ; top = be->ttop ; dt = (top-bot)/BNSUB ;
07503      for( nn=0 ; nn < be->nfunc ; nn++ ){
07504        ftop = 0.0f ;
07505        for( jj=0 ; jj <= BNSUB ; jj++ ){
07506          val = basis_funceval( be->bfunc[nn] , bot+jj*dt ) ;
07507          val = fabs(val) ; if( val > ftop ) ftop = val ;
07508        }
07509        if( ftop > 0.0f ) be->bfunc[nn].ffac = basis_normall / ftop ;
07510      }
07511    }
07512 
07513    free(scp); return be;
07514 }

06972 {
06973    if( x <= bot || x >= top ) return 0.0f ;
06974    return (float)sin(freq*(x-bot)) ;
06975 }

07002 {
07003    if( x <= 0.0f || x >= 25.0f ) return 0.0f ;
07004    return (float)(exp(-x)*( SPM_A1*pow(x,SPM_P1)
07005                            -SPM_A2*pow(x,SPM_P2) )) ;
07006 }

07009 {
07010    if( x <= 0.0f || x >= 25.0f ) return 0.0f ;
07011    return (float)(exp(-x)*( SPM_A1*pow(x,SPM_P1-1.0)*(SPM_P1-x)
07012                            -SPM_A2*pow(x,SPM_P2-1.0)*(SPM_P2-x) )) ;
07013 }

06937 {
06938    if( x <= bot || x >= top ) return 0.0f ;
06939    if( x <= mid )             return (x-bot)/(mid-bot) ;
06940                               return (top-x)/(top-mid) ;
06941 }

07525 {
07526    int nvox, ii, nf, allz, ts_length ;
07527    register int pp, ib ;
07528    float *wt , *tt , **hout , factor , **bb ;
07529    short *bar ;
07530    THD_3dim_dataset *in_dset = NULL;
07531    THD_3dim_dataset *out_dset = NULL;
07532    char *commandline , label[512] ;
07533    const float EPSILON = 1.0e-10 ;
07534 
07535    /* open input 3D+time dataset to get some parameters */
07536 
07537    in_dset = THD_open_dataset(option_data->input_filename);
07538    nvox    = in_dset->daxes->nxx * in_dset->daxes->nyy * in_dset->daxes->nzz;
07539 
07540    if( dt <= 0.0f ) dt = DSET_TR(in_dset) ;
07541    if( dt <= 0.0f ) dt = 1.0f ;
07542 
07543    /* create output dataset on the input as a model, then close the input */
07544 
07545    out_dset = EDIT_empty_copy( in_dset ) ;
07546    tross_Copy_History( in_dset , out_dset ) ;
07547    DSET_delete( in_dset ) ;
07548 
07549    /* historicize the output */
07550 
07551    commandline = tross_commandline( PROGRAM_NAME , argc , argv ) ;
07552    sprintf( label , "Impulse response: %s" , output_filename ) ;
07553    if( commandline != NULL )
07554      tross_multi_Append_History( out_dset , commandline,label,NULL ) ;
07555    else
07556      tross_Append_History ( out_dset, label);
07557    free((void *)commandline) ;
07558 
07559    ts_length = 1 + (int)ceil( (be->ttop - be->tbot)/dt ) ; /* 13 Apr 2005: +1 */
07560 
07561    /* modify the output dataset appropriately */
07562 
07563    (void ) EDIT_dset_items( out_dset,
07564                              ADN_prefix,      output_filename,
07565                              ADN_label1,      output_filename,
07566                              ADN_self_name,   output_filename,
07567                              ADN_malloc_type, DATABLOCK_MEM_MALLOC,
07568                              ADN_datum_all,   MRI_short,
07569                              ADN_nvals,       ts_length,
07570                              ADN_ntt,         ts_length,
07571                              ADN_ttdel,       dt ,
07572                              ADN_ttorg,       be->tbot,
07573                             ADN_none ) ;
07574 
07575    if( THD_is_file(out_dset->dblk->diskptr->header_name) ){
07576      fprintf(stderr,
07577              "** ERROR: Output dataset file %s already exists - won't overwrite\n",
07578              out_dset->dblk->diskptr->header_name ) ;
07579      DSET_delete(out_dset) ; return ;
07580    }
07581 
07582    /* create output bricks (float for now, will scale to shorts later) */
07583 
07584    hout = (float **) malloc( sizeof(float *) * ts_length ) ;
07585    for( ib=0 ; ib < ts_length ; ib++ )
07586      hout[ib] = (float *)calloc(sizeof(float),nvox) ;
07587 
07588    /* create basis vectors for the expansion on the dt time grid */
07589 
07590    nf = be->nfunc ;
07591    wt = (float *) malloc( sizeof(float) * nf ) ;
07592    tt = (float *) malloc( sizeof(float) * ts_length ) ;
07593 
07594    for( ib=0 ; ib < ts_length ; ib++ )     /* output time grid */
07595      tt[ib] = be->tbot + ib*dt ;
07596 
07597    bb = (float **) malloc( sizeof(float *) * nf ) ;
07598    for( pp=0 ; pp < nf ; pp++ ){
07599      bb[pp] = (float *) malloc( sizeof(float) * ts_length ) ;
07600      for( ib=0 ; ib < ts_length ; ib++ )
07601        bb[pp][ib] = basis_funceval( be->bfunc[pp] , tt[ib] ) ;
07602    }
07603 
07604    /* loop over voxels:
07605         extract coefficient (weights) for each basis function into wt
07606         sum up basis vectors times wt to get result, save into output arrays */
07607 
07608    for( ii=0 ; ii < nvox ; ii++ ){
07609      allz = 1 ;
07610      for( pp=0 ; pp < nf ; pp++ ){
07611        wt[pp] = wtar[pp][ii] ; allz = ( allz && (wt[pp] == 0.0f) );
07612      }
07613 
07614      if( allz ){
07615        for( ib=0 ; ib < ts_length ; ib++ ) hout[ib][ii] = 0.0f ;
07616      } else {
07617        register float sum ;
07618        for( ib=0 ; ib < ts_length ; ib++ ){
07619          sum = 0.0f ;
07620          for( pp=0 ; pp < nf ; pp++ ) sum += wt[pp] * bb[pp][ib] ;
07621          hout[ib][ii] = sum ;
07622        }
07623      }
07624    }
07625 
07626    /* toss some trash */
07627 
07628    for( pp=0 ; pp < nf ; pp++ ) free((void *)bb[pp]) ;
07629    free((void *)bb) ; free((void *)tt) ; free((void *)wt) ;
07630 
07631    /* scale floating point bricks to shorts and insert into output dataset */
07632 
07633    for( ib=0 ; ib < ts_length ; ib++ ){
07634      bar = (short *) malloc( sizeof(short) * nvox ) ;
07635      factor = EDIT_coerce_autoscale_new(nvox, MRI_float,hout[ib], MRI_short,bar) ;
07636      if( factor < EPSILON ) factor = 0.0f ;          /* if brick is all zero */
07637      else                   factor = 1.0f / factor ;
07638      EDIT_BRICK_FACTOR( out_dset , ib , factor ) ;
07639      EDIT_substitute_brick( out_dset , ib , MRI_short , bar ) ;
07640      free((void *)hout[ib]) ;
07641    }
07642    free((void *)hout) ;
07643 
07644    /* and save the results to disk! */
07645 
07646    DSET_write( out_dset ) ;
07647    if( verb )
07648     fprintf(stderr,"++ Wrote iresp 3D+time dataset into %s\n",DSET_BRIKNAME(out_dset)) ;
07649 
07650    DSET_delete( out_dset ) ;
07651    return ;
07652 }

07661 {
07662    int nvox, ii, nf, allz, ts_length ;
07663    register int pp,qq, ib ;
07664    register float sum ;
07665    float *vv , *tt , **hout , factor , **bb ;
07666    short *bar ;
07667    THD_3dim_dataset *in_dset = NULL;
07668    THD_3dim_dataset *out_dset = NULL;
07669    char *commandline , label[512] ;
07670    const float EPSILON = 1.0e-10 ;
07671 
07672    /* open input 3D+time dataset to get some parameters */
07673 
07674    in_dset = THD_open_dataset(option_data->input_filename);
07675    nvox    = in_dset->daxes->nxx * in_dset->daxes->nyy * in_dset->daxes->nzz;
07676 
07677    if( dt <= 0.0f ) dt = DSET_TR(in_dset) ;
07678    if( dt <= 0.0f ) dt = 1.0f ;
07679 
07680    /* create output dataset on the input as a model, then close the input */
07681 
07682    out_dset = EDIT_empty_copy( in_dset ) ;
07683    tross_Copy_History( in_dset , out_dset ) ;
07684    DSET_delete( in_dset ) ;
07685 
07686    /* historicize the output */
07687 
07688    commandline = tross_commandline( PROGRAM_NAME , argc , argv ) ;
07689    sprintf( label , "Sigma response: %s" , output_filename ) ;
07690    if( commandline != NULL )
07691      tross_multi_Append_History( out_dset , commandline,label,NULL ) ;
07692    else
07693      tross_Append_History ( out_dset, label);
07694    free((void *)commandline) ;
07695 
07696    ts_length = 1 + (int)ceil( (be->ttop - be->tbot)/dt ) ;
07697 
07698    /* modify the output dataset appropriately */
07699 
07700    (void ) EDIT_dset_items( out_dset,
07701                              ADN_prefix,      output_filename,
07702                              ADN_label1,      output_filename,
07703                              ADN_self_name,   output_filename,
07704                              ADN_malloc_type, DATABLOCK_MEM_MALLOC,
07705                              ADN_datum_all,   MRI_short,
07706                              ADN_nvals,       ts_length,
07707                              ADN_ntt,         ts_length,
07708                              ADN_ttdel,       dt ,
07709                              ADN_ttorg,       be->tbot,
07710                             ADN_none ) ;
07711 
07712    if( THD_is_file(out_dset->dblk->diskptr->header_name) ){
07713      fprintf(stderr,
07714              "** ERROR: Output dataset file %s already exists - won't overwrite\n",
07715              out_dset->dblk->diskptr->header_name ) ;
07716      DSET_delete(out_dset) ; return ;
07717    }
07718 
07719    /* create output bricks (float for now, will scale to shorts later) */
07720 
07721    hout = (float **) malloc( sizeof(float *) * ts_length ) ;
07722    for( ib=0 ; ib < ts_length ; ib++ )
07723      hout[ib] = (float *)calloc(sizeof(float),nvox) ;
07724 
07725    nf = be->nfunc ;
07726    tt = (float *) malloc( sizeof(float) * ts_length ) ;
07727 
07728    for( ib=0 ; ib < ts_length ; ib++ )     /* output time grid */
07729      tt[ib] = be->tbot + ib*dt ;
07730 
07731    /* evaluate basis vectors for output on dt time grid */
07732 
07733    bb = (float ** ) malloc( sizeof(float *) * nf ) ;
07734    for( pp=0 ; pp < nf ; pp++ ){
07735      bb[pp] = (float * ) malloc( sizeof(float  ) * ts_length ) ;
07736      for( ib=0 ; ib < ts_length ; ib++ )
07737        bb[pp][ib] = basis_funceval( be->bfunc[pp] , tt[ib] ) ;
07738    }
07739    free((void *)tt) ;
07740 
07741    /* evaluate unscaled variance on dt time grid */
07742 
07743    vv = (float *) malloc( sizeof(float) * ts_length ) ;
07744    for( ib=0 ; ib < ts_length ; ib++ ){
07745      sum = 0.0f ;
07746      for( pp=0 ; pp < nf ; pp++ ){
07747        for( qq=0 ; qq < nf ; qq++ )
07748          sum += cvar.elts[pbot+pp][pbot+qq] * bb[pp][ib] * bb[qq][ib] ;
07749      }
07750      vv[ib] = (sum >= 0.0f) ? sum : 0.0f ;
07751    }
07752    for( pp=0 ; pp < nf ; pp++ ) free((void *)bb[pp]) ;
07753    free((void *)bb) ;
07754 
07755    /* loop over voxels, scale by mse to get variance */
07756 
07757    for( ii=0 ; ii < nvox ; ii++ ){
07758      for( ib=0 ; ib < ts_length ; ib++ )
07759        hout[ib][ii] = sqrt( vv[ib] * mse[ii] ) ;
07760    }
07761    free((void *)vv) ;
07762 
07763    /* scale floating point bricks to shorts and insert into output dataset */
07764 
07765    for( ib=0 ; ib < ts_length ; ib++ ){
07766      bar = (short *) malloc( sizeof(short) * nvox ) ;
07767      factor = EDIT_coerce_autoscale_new(nvox, MRI_float,hout[ib], MRI_short,bar) ;
07768      if( factor < EPSILON ) factor = 0.0f ;          /* if brick is all zero */
07769      else                   factor = 1.0f / factor ;
07770      EDIT_BRICK_FACTOR( out_dset , ib , factor ) ;
07771      EDIT_substitute_brick( out_dset , ib , MRI_short , bar ) ;
07772      free((void *)hout[ib]) ;
07773    }
07774    free((void *)hout) ;
07775 
07776    /* and save the results to disk! */
07777 
07778    DSET_write( out_dset ) ;
07779    if( verb )
07780     fprintf(stderr,"++ Wrote sresp 3D+time dataset into %s\n",DSET_BRIKNAME(out_dset)) ;
07781 
07782    DSET_delete( out_dset ) ;
07783    return ;
07784 }

03698 {
03699   float * ts_array = NULL;    /* array of measured data for one voxel */
03700 
03701 #ifdef USE_GET
03702   int do_get = 0 ;            /* flag to use multi-gets */
03703 #endif
03704 
03705   int qp;                     /* number of poly. trend baseline parameters */
03706   int q;                      /* number of baseline model parameters */
03707   int p;                      /* number of full model parameters */
03708   int polort;                 /* degree of polynomial for baseline model */
03709   int m;                      /* parameter index */
03710   int n;                      /* data point index */
03711 
03712   vector coef;                /* regression parameters */
03713   vector scoef;               /* std. devs. for regression parameters */
03714   vector tcoef;               /* t-statistics for regression parameters */
03715   float * fpart = NULL;       /* partial F-statistics for the stimuli */
03716   float * rpart = NULL;       /* partial R^2 stats. for the stimuli */
03717   float ffull;                /* full model F-statistic */
03718   float rfull;                /* full model R^2 stat. */
03719   float mse;                  /* mean square error from full model */
03720 
03721   matrix xdata;               /* independent variable matrix */
03722   matrix x_full;              /* extracted X matrix    for full model */
03723   matrix xtxinv_full;         /* matrix:  1/(X'X)      for full model */
03724   matrix xtxinvxt_full;       /* matrix:  (1/(X'X))X'  for full model */
03725   matrix x_base;              /* extracted X matrix    for baseline model */
03726   matrix xtxinvxt_base;       /* matrix:  (1/(X'X))X'  for baseline model */
03727   matrix * x_rdcd = NULL;     /* extracted X matrices  for reduced models */
03728   matrix * xtxinvxt_rdcd = NULL;
03729                               /* matrix:  (1/(X'X))X'  for reduced models */
03730   vector y;                   /* vector of measured data */
03731 
03732   int ixyz;                   /* voxel index */
03733   int nxyz;                   /* number of voxels per dataset */
03734 
03735   int nt;                  /* number of images in input 3d+time dataset */
03736   int N;                   /* number of usable data points */
03737 
03738   int num_stimts;          /* number of stimulus time series */
03739   int * baseline;          /* flag for stim function in baseline model */
03740   int * min_lag;           /* minimum time delay for impulse response */
03741   int * max_lag;           /* maximum time delay for impulse response */
03742   int * nptr;              /* number of stim fn. time points per TR */
03743   char ** stim_label;      /* label for stimulus time series */
03744 
03745   int i;                   /* data point index */
03746   int is;                  /* stimulus index */
03747   int ilag;                /* time lag index */
03748   float stddev;            /* normalized parameter standard deviation */
03749   float rms_min;           /* minimum variation in data to fit full model */
03750   char * label;            /* string containing stat. summary of results */
03751   int nodata;              /* flag for 'no data' option */
03752   int novar;               /* flag for insufficient variation in data */
03753 
03754   int iglt;                    /* general linear test index */
03755   int ilc;                     /* linear combination index */
03756   int num_glt;                 /* number of general linear tests */
03757   char ** glt_label;           /* label for general linear test */
03758   int * glt_rows;              /* number of linear constraints in glt */
03759   matrix * cxtxinvct = NULL;   /* matrices: C(1/(X'X))C' for GLT */
03760   matrix * glt_amat = NULL;    /* constant GLT matrices for later use */
03761   vector * glt_coef = NULL;    /* linear combinations from GLT matrices */
03762   vector * glt_tcoef = NULL;   /* t-statistics for GLT linear combinations */
03763   float * fglt = NULL;         /* F-statistics for the general linear tests */
03764   float * rglt = NULL;         /* R^2 stats. for the general linear tests */
03765 
03766   float * fitts = NULL;        /* full model fitted time series */
03767   float * errts = NULL;        /* full model residual error time series */
03768 
03769   int vstep ;                  /* interval progress meter dots */
03770   double ct ;                  /* clock time */
03771 
03772 ENTRY("calculate_results") ;
03773 
03774   /*----- Initialize local variables -----*/
03775   nodata = option_data->nodata;
03776   nxyz = option_data->nxyz;
03777   nt = option_data->nt;
03778   rms_min = option_data->rms_min;
03779   num_stimts = option_data->num_stimts;
03780   stim_label = option_data->stim_label;
03781   baseline = option_data->stim_base;
03782   min_lag = option_data->stim_minlag;
03783   max_lag = option_data->stim_maxlag;
03784   nptr    = option_data->stim_nptr;
03785   num_glt = option_data->num_glt;
03786   glt_label = option_data->glt_label;
03787   glt_rows = option_data->glt_rows;
03788 
03789   polort = option_data->polort;
03790   qp = option_data->qp;
03791   q  = option_data->q;
03792   p  = option_data->p;
03793 
03794   N = option_data->N;
03795 
03796 
03797   /*----- Initialize matrices and vectors -----*/
03798   matrix_initialize (&xdata);
03799   matrix_initialize (&x_full);
03800   matrix_initialize (&xtxinv_full);
03801   matrix_initialize (&xtxinvxt_full);
03802   matrix_initialize (&x_base);
03803   matrix_initialize (&xtxinvxt_base);
03804   vector_initialize (&coef) ;
03805   vector_initialize (&scoef);
03806   vector_initialize (&tcoef);
03807   vector_initialize (&y);
03808 
03809   if (num_stimts > 0)
03810     {
03811       x_rdcd = (matrix *) malloc (sizeof(matrix) * num_stimts);
03812       MTEST (x_rdcd);
03813       xtxinvxt_rdcd = (matrix *) malloc (sizeof(matrix) * num_stimts);
03814       MTEST (xtxinvxt_rdcd);
03815 
03816       for (is =0;  is < num_stimts;  is++)
03817       {
03818         matrix_initialize (&x_rdcd[is]);
03819         matrix_initialize (&xtxinvxt_rdcd[is]);
03820       }
03821     }
03822 
03823   if (num_glt > 0)
03824     {
03825       cxtxinvct = (matrix *) malloc (sizeof(matrix) * num_glt);
03826       glt_amat  = (matrix *) malloc (sizeof(matrix) * num_glt);
03827       glt_coef  = (vector *) malloc (sizeof(vector) * num_glt);
03828       glt_tcoef = (vector *) malloc (sizeof(vector) * num_glt);
03829 
03830       for (iglt =0;  iglt < num_glt;  iglt++)
03831       {
03832         matrix_initialize (&cxtxinvct[iglt]);
03833         matrix_initialize (&glt_amat[iglt]);
03834         vector_initialize (&glt_coef[iglt]);
03835         vector_initialize (&glt_tcoef[iglt]);
03836       }
03837     }
03838 
03839 
03840   /*----- Allocate memory -----*/
03841   if (num_stimts > 0)
03842     {
03843       fpart  = (float *) malloc (sizeof(float) * num_stimts);   MTEST (fpart);
03844       rpart  = (float *) malloc (sizeof(float) * num_stimts);   MTEST (rpart);
03845     }
03846   if (num_glt > 0)
03847     {
03848       fglt  = (float *) malloc (sizeof(float) * num_glt);   MTEST (fglt);
03849       rglt  = (float *) malloc (sizeof(float) * num_glt);   MTEST (rglt);
03850     }
03851 
03852   if (option_data->input1D_filename == NULL)
03853     {
03854 #ifdef USE_GET
03855       do_get = 1 ;    /* don't need a pre-malloc-ed array for multi-gets */
03856 #else
03857       ts_array = (float *) malloc (sizeof(float) * nt);   MTEST (ts_array);
03858 #endif
03859     }
03860   fitts  = (float *) malloc (sizeof(float) * nt);   MTEST (fitts);
03861   errts  = (float *) malloc (sizeof(float) * nt);   MTEST (errts);
03862 
03863 
03864   /*----- Initialize the independent variable matrix -----*/
03865 
03866   ct = COX_clock_time() ;
03867 
03868   init_indep_var_matrix (p, qp, polort, nt, N, good_list, block_list,
03869                    num_blocks, num_stimts, stimulus, stim_length,
03870                    min_lag, max_lag, nptr, option_data->stim_base , &xdata);
03871   if (option_data->xout)  matrix_sprint ("X matrix:", xdata);
03872 
03873   if( option_data->xjpeg_filename != NULL )    /* 21 Jul 2004 */
03874     JPEG_matrix_gray( xdata , option_data->xjpeg_filename ) ;
03875 
03876 
03877   /*-- 14 Jul 2004: check matrix for bad columns - RWCox --*/
03878 
03879   { int *iar , k , nerr=0 ;
03880     iar = matrix_check_columns( xdata , 1.e-7 ) ;
03881     if( iar != NULL ){
03882       fprintf(stderr,"** WARNING: Problems with the X matrix columns:\n") ;
03883       for( k=0 ; iar[2*k] >= 0 ; k++ ){
03884         if( iar[2*k+1] >= 0 ){
03885           fprintf(stderr," * Columns %d and %d are nearly collinear!\n",
03886                   iar[2*k],iar[2*k+1] ) ; nerr++ ;
03887         } else {
03888           fprintf(stderr," * Column %d is all zeros: %s\n",
03889                   iar[2*k] ,
03890                   use_psinv ? "SVD on => will be kept"
03891                             : "SVD off => will be excised" ) ;
03892         }
03893       }
03894       if( nerr > 0 && AFNI_yesenv("AFNI_3dDeconvolve_nodup") ){
03895         fprintf(stderr,"** ERROR: can't continue after above warnings!\n");
03896         exit(1) ;
03897       }
03898       free(iar) ;
03899     }
03900   }
03901 
03902   /*-- 14 Jul 2004: calculate matrix condition number - RWCox --*/
03903 
03904 #ifndef PSINV_EPS
03905 #define PSINV_EPS 1.e-12
03906 #endif
03907   if( !option_data->nocond ){
03908     double *ev , ebot,emin,emax ; int i,nsmall=0 ;
03909     ev = matrix_singvals( xdata ) ;
03910     emax = 1.e-38 ;
03911     for( i=0 ; i < xdata.cols ; i++ ){
03912       if( ev[i] > emax ) emax = ev[i] ;
03913     }
03914     ebot = sqrt(PSINV_EPS)*emax ; emin = 1.e+38 ;
03915     for( i=0 ; i < xdata.cols ; i++ ){
03916       if( ev[i] >= ebot && ev[i] < emin ) emin = ev[i] ;
03917       if( ev[i] <  ebot ) nsmall++ ;
03918     }
03919     if( nsmall > 0 ){
03920       fprintf(stderr,
03921         "** WARNING: Largest singular value=%g;"
03922                   " %d %s less than cutoff=%g\n"
03923         "            Implies strong collinearity in the input regressors\n",
03924               emax , nsmall , (nsmall==1)?"is":"are" , ebot ) ;
03925       show_singvals = 1 ;
03926     }
03927     if( show_singvals ){
03928       fprintf(stderr,"++ Matrix singular values:") ;
03929       for( i=0; i < xdata.cols ; i++ ) fprintf(stderr," %g",ev[i]) ;
03930       fprintf(stderr,"\n") ;
03931     }
03932     free((void *)ev) ;
03933     if( emin <= 0.0 || emax <= 0.0 ){
03934       fprintf(stderr,"** Matrix condition:  UNDEFINED: "
03935                      "min ev=%g  max ev=%g  ** VERY BAD **\n",emin,emax ) ;
03936     } else {
03937       double cond = sqrt(emax/emin) ;
03938       if( !use_psinv ) cond = cond*cond ;  /* Gaussian elim is twice as bad */
03939       fprintf(stderr,"++ (%dx%d) Matrix condition [%s]:  %g",
03940               xdata.rows,xdata.cols , (use_psinv) ? "X" : "XtX" , cond  ) ;
03941 #ifdef FLOATIZE
03942       if( cond > 100.0 ) fprintf(stderr,"  ** BEWARE **") ;
03943 #else
03944       if( cond > 1.e6  ) fprintf(stderr,"  ** BEWARE **") ;
03945 #endif
03946       fprintf(stderr,"\n") ;
03947     }
03948   }
03949 
03950   /*----- Initialization for the regression analysis -----*/
03951   init_regression_analysis (p, qp, num_stimts, baseline, min_lag, max_lag,
03952                       xdata, &x_full, &xtxinv_full, &xtxinvxt_full,
03953                       &x_base, &xtxinvxt_base, x_rdcd, xtxinvxt_rdcd);
03954   if ((option_data->xout) || nodata)
03955       matrix_sprint ("(X'X) inverse matrix:", xtxinv_full);
03956 
03957   /*----- Save some of this stuff for later, dude -----*/
03958   X        = x_full        ;  /* 25 Jul 2004 (RWCox) */
03959   XtXinv   = xtxinv_full   ;
03960   XtXinvXt = xtxinvxt_full ;
03961 
03962   { int m , npar , j ;        /* 31 Aug 2004 (RWCox) */
03963     register double sum ;
03964     float mmax ;
03965 
03966     Xcol_inbase = (int *)  calloc(sizeof(int)  ,p) ;
03967     Xcol_mean   = (float *)calloc(sizeof(float),p) ;
03968 
03969     for( is=0 ; is < qp ; is++ ) Xcol_inbase[is] = 1 ; /* mark baseline columns */
03970     m = qp ;
03971     for( is=0 ; is < num_stimts ; is++ ){
03972       npar = (basis_stim[is] != NULL)
03973             ? basis_stim[is]->nfunc
03974             : option_data->stim_maxlag[is] - option_data->stim_minlag[is] + 1 ;
03975 
03976       if( baseline[is] ) for( j=0 ; j < npar ; j++ ) Xcol_inbase[m+j] = 1 ;
03977       m += npar ;
03978     }
03979 
03980     mmax = 0.0f ;
03981     for( j=0 ; j < p ; j++ ){   /* compute mean of each column */
03982       sum = 0.0 ;
03983       for( i=0 ; i < X.rows ; i++ ) sum += X.elts[i][j] ;
03984       Xcol_mean[j] = (float)(sum/X.rows) ;
03985       if( Xcol_inbase[j] && fabs(Xcol_mean[j]) > mmax )  /* find largest */
03986         mmax = fabs(Xcol_mean[j]) ;             /* mean of baseline cols */
03987     }
03988 
03989     if( mmax > 0.0f ){    /* mark baseline cols that have nontrivial means */
03990       mmax *= 9.99e-6 ;
03991       for( j=0 ; j < p ; j++ )
03992         if( Xcol_inbase[j] && fabs(Xcol_mean[j]) > mmax ) Xcol_inbase[j] = 2 ;
03993     }
03994   }
03995 
03996   /*--- Compute abs sum of matrix [xtxinvxt][xdata]-I [19 Aug 2004]---*/
03997   if( !option_data->nocond ){
03998     double esum , sum ;
03999     int nn=xdata.rows , mm=xdata.cols , ii,jj,kk ;
04000     char *www = "\0" ;
04001     fprintf(stderr,"++ Matrix inverse average error = ") ;
04002     esum = 0.0 ;
04003     for( ii=0 ; ii < mm ; ii++ ){
04004       for( jj=0 ; jj < mm ; jj++ ){
04005         sum = (ii==jj) ? -1.0 : 0.0 ;
04006         for( kk=0 ; kk < nn ; kk++ )
04007           sum += xtxinvxt_full.elts[ii][kk]*xdata.elts[kk][jj] ;
04008         esum += fabs(sum) ;
04009       }
04010     }
04011     esum /= (mm*mm) ;
04012     if( esum > 1.e-3 ) www = " ** WARNING!!! **" ;
04013     fprintf(stderr,"%g %s\n",esum,www) ;
04014   }
04015 
04016   /*-- 19 Aug 2004: plot matrix pseudoinverse as well --*/
04017   if( option_data->xjpeg_filename != NULL ){
04018     char *jpt , *jsuf=".jpg" ;
04019     char *fn = calloc( sizeof(char) , strlen(option_data->xjpeg_filename)+16 ) ;
04020     matrix xpsinv ;
04021 
04022     matrix_initialize( &xpsinv ) ;
04023     matrix_transpose( xtxinvxt_full , &xpsinv ) ;
04024 
04025     strcpy(fn,option_data->xjpeg_filename) ;
04026                        jpt = strstr(fn,".jpg") ;
04027     if( jpt == NULL ){ jpt = strstr(fn,".JPG") ; jsuf = ".JPG" ; }
04028     if( jpt == NULL )  jpt = fn + strlen(fn) ;
04029     strcpy(jpt,"_psinv") ; strcat(fn,jsuf) ;
04030 
04031     JPEG_matrix_gray( xpsinv , fn ) ;
04032     free((void *)fn) ; matrix_destroy( &xpsinv ) ;
04033   }
04034 
04035   /*----- Initialization for the general linear test analysis -----*/
04036   if (num_glt > 0)
04037     init_glt_analysis (xtxinv_full, num_glt, glt_cmat, glt_amat, cxtxinvct);
04038 
04039   ct = COX_clock_time() - ct ;
04040   fprintf(stderr,"++ Matrix setup time = %.2f s\n",ct) ;  /* 25 Apr 2005 */
04041 
04042 
04043   vector_create (N, &y);
04044 
04045 
04046   if (nodata)
04047     {
04048       report_evaluation(qp, num_stimts, stim_label, min_lag, max_lag,
04049                         x_full, xtxinv_full,
04050                         num_glt, glt_label, glt_rows, cxtxinvct);
04051     }
04052 
04053   else
04054     {     /*============== actually process data ====================*/
04055 
04056       int ixyz_bot=0 , ixyz_top=nxyz ;  /* voxel indexes to process */
04057 
04058 #ifdef USE_GET
04059 #define NGET 128              /* number to get at one time */
04060       int nget=0 ,            /* number of time series current gotten */
04061           cget=0 ,            /* index of next timeseries in iget & imget */
04062           jget   ,            /* loop index for iget */
04063           iget[NGET] ;        /* voxel index of timeseries */
04064       MRI_IMARR *imget=NULL ; /* array of timeseries */
04065 #endif
04066 
04067 #ifdef PROC_MAX
04068       if( proc_numjob > 1 ){    /*---- set up multiple processes ----*/
04069         int vv , nvox=nxyz , nper , pp , nv ;
04070         pid_t newpid ;
04071 
04072         /* count number of voxels to compute with into nvox */
04073 
04074         if( mask_vol != NULL ){
04075           for( vv=nvox=0 ; vv < nxyz ; vv++ )
04076             if( mask_vol[vv] != 0 ) nvox++ ;
04077         }
04078 
04079         if( nvox < proc_numjob ){  /* too few voxels for multiple jobs? */
04080 
04081           proc_numjob = 1 ;
04082 
04083         } else {                   /* prepare jobs */
04084 
04085           /* split voxels between jobs evenly */
04086 
04087           nper = nvox / proc_numjob ;  /* # voxels per job */
04088           if( mask_vol == NULL ){
04089             proc_vox_bot[0] = 0 ;
04090             for( pp=0 ; pp < proc_numjob ; pp++ ){
04091               proc_vox_top[pp] = proc_vox_bot[pp] + nper ;
04092               if( pp < proc_numjob-1 ) proc_vox_bot[pp+1] = proc_vox_top[pp] ;
04093             }
04094             proc_vox_top[proc_numjob-1] = nxyz ;
04095           } else {
04096             proc_vox_bot[0] = 0 ;
04097             for( pp=0 ; pp < proc_numjob ; pp++ ){
04098               for( nv=0,vv=proc_vox_bot[pp] ;         /* count ahead until */
04099                    nv < nper && vv < nxyz  ; vv++ ){  /* find nper voxels */
04100                 if( mask_vol[vv] != 0 ) nv++ ;        /* inside the mask */
04101               }
04102               proc_vox_top[pp] = vv ;
04103               if( pp < proc_numjob-1 ) proc_vox_bot[pp+1] = proc_vox_top[pp] ;
04104             }
04105             proc_vox_top[proc_numjob-1] = nxyz ;
04106           }
04107 
04108           /* make sure dataset is in memory before forks */
04109 
04110           DSET_load(dset) ;  /* so dataset will be common */
04111 
04112           /* start processes */
04113 
04114           if( !option_data->quiet ) fprintf(stderr,"++ Voxels in dataset: %d\n",nxyz) ;
04115           if( nvox < nxyz )
04116           if( !option_data->quiet ) fprintf(stderr,"++ Voxels in mask:    %d\n",nvox) ;
04117           if( !option_data->quiet ) fprintf(stderr,"++ Voxels per job:    %d\n",nper) ;
04118 
04119           for( pp=1 ; pp < proc_numjob ; pp++ ){
04120             ixyz_bot = proc_vox_bot[pp] ;   /* these 3 variables   */
04121             ixyz_top = proc_vox_top[pp] ;   /* are for the process */
04122             proc_ind = pp ;                 /* we're about to fork */
04123             newpid   = fork() ;
04124             if( newpid == -1 ){
04125               fprintf(stderr,"** Can't fork job #%d! Error exit!\n",pp);
04126               exit(1) ;
04127             }
04128             if( newpid == 0 ) break ;   /* I'm the child */
04129             proc_pid[pp] = newpid ;     /* I'm the parent */
04130             iochan_sleep(10) ;          /* 10 ms, to let the child get going */
04131           }
04132           if( pp == proc_numjob ){       /* only in the parent */
04133             ixyz_bot = proc_vox_bot[0] ; /* set the 3 control */
04134             ixyz_top = proc_vox_top[0] ; /* variables needed */
04135             proc_ind = 0 ;               /* below           */
04136           }
04137           if( !option_data->quiet )
04138             fprintf(stderr,"++ Job #%d: processing voxels %d to %d; elapsed time=%.3f\n",
04139                     proc_ind,ixyz_bot,ixyz_top-1,COX_clock_time()) ;
04140         }
04141       }
04142 #endif /* PROC_MAX */
04143 
04144       if( proc_numjob == 1 && !option_data->quiet )
04145         fprintf(stderr,"++ Calculations starting; elapsed time=%.3f\n",COX_clock_time()) ;
04146 
04147       vstep = nxyz / 50 ;
04148       if( option_data->quiet        ||
04149           option_data->fdisp >= 0.0 ||
04150           option_data->progress > 0 ||
04151           proc_numjob > 1             ) vstep = 0 ;
04152 
04153       if( vstep > 0 ) fprintf(stderr,"++ voxel loop:") ;
04154 
04155       /*----- Loop over all voxels -----*/
04156       for (ixyz = ixyz_bot;  ixyz < ixyz_top;  ixyz++)
04157       {
04158 
04159           if( vstep > 0 && ixyz%vstep==vstep-1 ) vstep_print() ;
04160 
04161         /*----- Apply mask? -----*/
04162         if (mask_vol != NULL)
04163           if (mask_vol[ixyz] == 0)  continue;
04164 
04165 #ifdef USE_GET
04166           /*** race ahead and extract a bunch of voxel time series at once ***/
04167 
04168           if( do_get && cget == nget ){
04169             if( imget != NULL ) DESTROY_IMARR(imget) ;
04170             iget[0] = ixyz ; nget = 1 ;
04171             for( jget=ixyz+1 ; jget < nxyz && nget < NGET ; jget++ ){
04172               if( mask_vol == NULL || mask_vol[jget] != 0 )
04173                 iget[nget++] = jget ;
04174             }
04175             imget = THD_extract_many_series( nget, iget, dset ) ;
04176             cget  = 0 ;  /* the next one to take out of imget */
04177           }
04178 #endif
04179 
04180         /*----- Extract Y-data for this voxel -----*/
04181         if (option_data->input1D_filename != NULL)
04182           ts_array = fmri_data;
04183         else {
04184 #ifdef USE_GET
04185             ts_array = MRI_FLOAT_PTR(IMARR_SUBIM(imget,cget));  /* the GET way */
04186             cget++ ;  /* take this one next time */
04187 #else
04188             extract_ts_array (dset, ixyz, ts_array);            /* the OLD way */
04189 #endif
04190           }
04191 
04192         for (i = 0;  i < N;  i++)
04193           y.elts[i] = ts_array[good_list[i]];
04194 
04195 
04196         /*----- Perform the regression analysis for this voxel-----*/
04197         regression_analysis (N, p, q, num_stimts, min_lag, max_lag,
04198                          x_full, xtxinv_full, xtxinvxt_full, x_base,
04199                          xtxinvxt_base, x_rdcd, xtxinvxt_rdcd,
04200                          y, rms_min, &mse, &coef, &scoef, &tcoef,
04201                          fpart, rpart, &ffull, &rfull, &novar,
04202                          fitts, errts);
04203 
04204           if( voxel_base != NULL )
04205             voxel_base[ixyz] = baseline_mean( coef ) ;   /* 31 Aug 2004 */
04206 
04207 
04208         /*----- Perform the general linear tests for this voxel -----*/
04209         if (num_glt > 0)
04210           glt_analysis (N, p, x_full, y, mse*(N-p), coef, novar, cxtxinvct,
04211                     num_glt, glt_rows, glt_cmat, glt_amat,
04212                     glt_coef, glt_tcoef, fglt, rglt);
04213 
04214 
04215         /*----- Save results for this voxel into arrays -----*/
04216         save_voxel (option_data, ixyz, coef, scoef, tcoef, fpart, rpart, mse,
04217                   ffull, rfull, glt_coef, glt_tcoef, fglt, rglt,
04218                   nt, ts_array, good_list, fitts, errts,
04219                   coef_vol, scoef_vol, tcoef_vol, fpart_vol, rpart_vol,
04220                   mse_vol, ffull_vol, rfull_vol, glt_coef_vol,
04221                   glt_tcoef_vol, glt_fstat_vol, glt_rstat_vol,
04222                   fitts_vol, errts_vol);
04223 
04224 
04225         /*----- Report results for this voxel -----*/
04226         if ( ((ffull > option_data->fdisp) && (option_data->fdisp >= 0.0))
04227              || ((option_data->progress > 0)
04228                && (ixyz % option_data->progress == 0))
04229              || (option_data->input1D_filename != NULL) )
04230           {
04231 
04232               if( proc_ind == 0 ){
04233               printf ("\n\nResults for Voxel #%d: \n", ixyz);
04234               report_results (N, qp, q, p, polort, block_list, num_blocks,
04235                        num_stimts, stim_label, baseline, min_lag, max_lag,
04236                        coef, tcoef, fpart, rpart, ffull, rfull, mse,
04237                          num_glt, glt_label, glt_rows, glt_coef,
04238                          glt_tcoef, fglt, rglt, &label);
04239               printf ("%s \n", label); fflush(stdout);
04240               }
04241           }
04242 
04243       }  /*----- Loop over voxels -----*/
04244 
04245         if( vstep > 0 ) fprintf(stderr,"\n") ;
04246 
04247 #ifdef USE_GET
04248         if( do_get ){
04249           if( imget != NULL ) DESTROY_IMARR(imget) ;
04250           ts_array = NULL ;
04251         }
04252 #endif
04253 
04254         /*-- if this is a child process, we're done.
04255              if this is the parent process, wait for the children --*/
04256 
04257 #ifdef PROC_MAX
04258         if( proc_numjob > 1 ){
04259           if( proc_ind > 0 ){                          /* death of child */
04260             if( !option_data->quiet )
04261               fprintf(stderr,"++ Job #%d finished; elapsed time=%.3f\n",proc_ind,COX_clock_time()) ;
04262             _exit(0) ;
04263 
04264           } else {                      /* parent waits for children */
04265             int pp ;
04266             if( !option_data->quiet )
04267               fprintf(stderr,"++ Job #0 waiting for children to finish; elapsed time=%.3f\n",COX_clock_time()) ;
04268             for( pp=1 ; pp < proc_numjob ; pp++ )
04269               waitpid( proc_pid[pp] , NULL , 0 ) ;
04270             if( !option_data->quiet )
04271               fprintf(stderr,"++ Job #0 now finishing up; elapsed time=%.3f\n",COX_clock_time()) ;
04272           }
04273 
04274           /* when get to here, only parent process is left alive,
04275              and all the results are in the shared memory segment arrays */
04276         }
04277 #endif
04278         if( proc_numjob == 1 && !option_data->quiet )
04279           fprintf(stderr,"++ Calculations finished; elapsed time=%.3f\n",COX_clock_time()) ;
04280 
04281         if( option_data->input1D_filename == NULL)  /* don't need data anymore */
04282           DSET_unload(dset) ;
04283 
04284     }  /*----- NOT nodata -----*/
04285 
04286 
04287   /*----- Dispose of matrices and vectors -----*/
04288   vector_destroy (&y);
04289   vector_destroy (&tcoef);
04290   vector_destroy (&scoef);
04291   vector_destroy (&coef);
04292 
04293   if (num_stimts > 0)
04294     {
04295       for (is = 0;  is < num_stimts;  is++)
04296       {
04297        matrix_destroy (&x_rdcd[is]);
04298        matrix_destroy (&xtxinvxt_rdcd[is]);
04299       }
04300       free (x_rdcd);         x_rdcd = NULL;
04301       free (xtxinvxt_rdcd);  xtxinvxt_rdcd = NULL;
04302     }
04303 
04304   matrix_destroy (&xtxinvxt_base);
04305   matrix_destroy (&x_base);
04306   matrix_destroy (&xdata);
04307 
04308   if (num_glt > 0)
04309     {
04310       for (iglt = 0;  iglt < num_glt;  iglt++)
04311       {
04312         matrix_destroy (&cxtxinvct[iglt]);
04313         matrix_destroy (&glt_amat[iglt]);
04314         vector_destroy (&glt_coef[iglt]);
04315         vector_destroy (&glt_tcoef[iglt]);
04316       }
04317       free (cxtxinvct);   cxtxinvct = NULL;
04318       free (glt_amat);    glt_amat  = NULL;
04319       free (glt_coef);    glt_coef  = NULL;
04320       free (glt_tcoef);   glt_tcoef = NULL;
04321     }
04322 
04323 
04324   if (fpart != NULL)     { free (fpart);     fpart    = NULL; }
04325   if (rpart != NULL)     { free (rpart);     rpart    = NULL; }
04326   if (fglt  != NULL)     { free (fglt);      fglt     = NULL; }
04327   if (rglt  != NULL)     { free (rglt);      rglt     = NULL; }
04328   if (ts_array != NULL)  { free (ts_array);  ts_array = NULL; }
04329   if (fitts != NULL)     { free (fitts);     fitts    = NULL; }
04330   if (errts != NULL)     { free (errts);     errts    = NULL; }
04331 
04332   EXRETURN ;
04333 }

02599                                                 : stimulus time series arrays */
02600   int ** good_list                /* list of usable time points */
02601 )
02602 
02603 {
02604   char message[THD_MAX_NAME];  /* error message */
02605   int is;                  /* stimulus index */
02606   int num_stimts;          /* number of stimulus time series */
02607   int * min_lag;           /* minimum time delay for impulse response */
02608   int * max_lag;           /* maximum time delay for impulse response */
02609   int * nptr;              /* number of stim fn. time points per TR */
02610   int m;                   /* number of time delays for impulse response */
02611   int qp;                  /* number of polynomial trend baseline parameters */
02612   int q;                   /* number of baseline model parameters */
02613   int p;                   /* number of full model parameters */
02614   int nbricks;             /* number of sub-bricks in bucket dataset output */
02615   int it;                  /* time point index */
02616   int nt;                  /* number of images in input 3d+time dataset */
02617   int NFirst;              /* first image from input 3d+time dataset to use */
02618   int NLast;               /* last image from input 3d+time dataset to use */
02619   int N;                   /* number of usable time points */
02620   int ib;                  /* block (run) index */
02621   int irb;                 /* time index relative to start of block (run) */
02622   int num_glt;             /* number of general linear tests */
02623   int * glt_rows;          /* number of linear constraints in glt */
02624   int iglt;                /* general linear test index */
02625   int nerr=0 ;             /* 22 Oct 2003 */
02626 
02627 
02628   /*----- Initialize local variables -----*/
02629   nt = option_data->nt;
02630   num_stimts = option_data->num_stimts;
02631   min_lag = option_data->stim_minlag;
02632   max_lag = option_data->stim_maxlag;
02633   num_glt = option_data->num_glt;
02634   glt_rows = option_data->glt_rows;
02635   nptr = option_data->stim_nptr;
02636   p  = option_data->p;
02637   q  = option_data->q;
02638   qp = option_data->qp;
02639 
02640   /*----- Check if -xsave was given without -bucket ------*/
02641   if( xsave && option_data->bucket_filename == NULL ){
02642     fprintf(stderr,"** WARNING: -xsave given without -bucket; -xsave is disabled!\n") ;
02643     xsave = 0 ;
02644   }
02645 
02646   /*----- Check length of censor array -----*/
02647   if (censor_length < nt)
02648     {
02649       sprintf (message, "Input censor time series file %s is too short",
02650              option_data->censor_filename);
02651       DC_error (message);
02652     }
02653 
02654 
02655   /*----- Check validity of concatenated runs list -----*/
02656   for (ib = 0;  ib < num_blocks;  ib++)
02657     if ((block_list[ib] < 0) || (block_list[ib] >= nt))
02658       {
02659       sprintf (message, "Invalid -concat input: %d ", block_list[ib]);
02660       DC_error (message);
02661       }
02662   if (num_blocks > 1)
02663     for (ib = 1;  ib < num_blocks;  ib++)
02664       if (block_list[ib] <= block_list[ib-1])
02665       DC_error ("Invalid concatenated runs list");
02666 
02667 
02668   /*----- Create list of good (usable) time points -----*/
02669   *good_list = (int *) malloc (sizeof(int) * nt);  MTEST (*good_list);
02670   NFirst = option_data->NFirst;
02671   if (NFirst < 0){
02672     for (is = 0;  is < num_stimts;  is++)
02673       if( basis_stim[is] == NULL && NFirst < (max_lag[is]+nptr[is]-1)/nptr[is] )
02674         NFirst = (max_lag[is]+nptr[is]-1)/nptr[is];
02675   }
02676   NLast = option_data->NLast;
02677   if (NLast < 0)  NLast = nt;
02678 
02679   N = 0;
02680   ib = 0;
02681   for (it = block_list[0];  it < nt;  it++)
02682     {
02683       if (ib+1 < num_blocks)
02684       if (it >= block_list[ib+1])  ib++;
02685 
02686       irb = it - block_list[ib];
02687 
02688       if ((irb >= NFirst) && (irb <= NLast) && (censor_array[it]))
02689       {
02690         (*good_list)[N] = it;
02691         N++;
02692       }
02693     }
02694 
02695 
02696   /*----- Check for sufficient data -----*/
02697   if (N == 0)  DC_error ("No usable time points?");
02698   if (N <= p)
02699     {
02700        sprintf (message,  "Insufficient data for estimating %d parameters", p);
02701        DC_error (message);
02702    }
02703   option_data->N = N;
02704 
02705    GoodList = *good_list ;
02706   nGoodList = N ;
02707   nParam    = p ;
02708 
02709   /*----- Check number of stimulus time series -----*/
02710   if (num_stimts < 0)
02711     {
02712       DC_error ("Require: 0 <= num_stimts (DUH! --or-- D'oh!)");
02713     }
02714 
02715 
02716   /*----- Check lengths of stimulus time series -----*/
02717 
02718 #define ALLOW_EXTEND
02719   for (is = 0;  is < num_stimts;  is++)
02720     {
02721       if( basis_stim[is] != NULL ) continue ;  /* 12 Aug 2004 */
02722 
02723       if (stim_length[is] < nt*nptr[is])
02724       {
02725 #ifndef ALLOW_EXTEND
02726         sprintf (message, "Input stimulus time series file %s is too short",
02727                option_data->stim_filename[is]);
02728         DC_error (message);
02729 #else
02730           int nlen=nt*nptr[is], qq ;
02731           fprintf(stderr,
02732                   "** WARNING: input stimulus time series file %s is too short:\n"
02733                   "            length = %d, but should be at least %d.\n" ,
02734             option_data->stim_filename[is] , stim_length[is] , nlen ) ;
02735           stimulus[is] = (float *) realloc( stimulus[is] , sizeof(float)*nlen ) ;
02736           for( qq=stim_length[is] ; qq < nlen ; qq++ ) stimulus[is][qq] = 0.0 ;
02737           stim_length[is] = nlen ; nerr++ ;
02738 #endif
02739       }
02740     }
02741 #ifdef ALLOW_EXTEND
02742     if( nerr > 0 ){
02743       char *eee = getenv("AFNI_3dDeconvolve_extend") ;
02744       if( eee != NULL && (*eee=='n' || *eee=='N') ){
02745         fprintf(stderr,"** ERROR: Can't continue with too short files!\n"
02746                        "          AFNI_3dDeconvolve_extend = %s\n",eee ) ;
02747         exit(1) ;
02748       }
02749       fprintf(stderr,"++ EXTENDING short files with zero values\n"
02750                      "   (to stop this behavior, setenv AFNI_3dDeconvolve_extend NO)\n") ;
02751     }
02752 #endif
02753 
02754 
02755   /*----- Check whether time lags are reasonable -----*/
02756   for (is = 0;  is < num_stimts;  is++)
02757     {
02758       if( basis_stim[is] != NULL ) continue ; /* 12 Aug 2004 */
02759 
02760       m = max_lag[is] - min_lag[is] + 1;
02761       if (m < 2)
02762       {
02763         if (option_data->iresp_filename[is] != NULL)
02764           {
02765             sprintf (message, "Only %d time point for output dataset %s ",
02766                    m, option_data->iresp_filename[is]);
02767             DC_error (message);
02768           }
02769 
02770         if (option_data->sresp_filename[is] != NULL)
02771           {
02772             sprintf (message, "Only %d time point for output dataset %s",
02773                    m, option_data->sresp_filename[is]);
02774             DC_error (message);
02775           }
02776       }
02777       if ((m < 4) && (option_data->tshift))
02778       {
02779         if (option_data->iresp_filename[is] != NULL)
02780           {
02781             sprintf (message, "Only %d time points for 3d+time dataset %s\n",
02782                     m, option_data->iresp_filename[is]);
02783             strcat (message, "Require >= 4 data points for -tshift option");
02784             DC_error (message);
02785           }
02786       }
02787     }
02788 
02789 
02790   /*----- Calculate number of sub-bricks in the bucket dataset,
02791           and check for illegal number of sub-bricks -----*/
02792   nbricks = 0;
02793   if (option_data->bucket_filename != NULL)
02794     {
02795       if (! option_data->nocout)
02796       {
02797         if (! option_data->nobout)
02798           nbricks += qp * (1 + option_data->tout);
02799 
02800         for (is = 0;  is < num_stimts;  is++)
02801           {
02802             if ((!option_data->stim_base[is]) || (!option_data->nobout))
02803             {
02804                   if( basis_stim[is] != NULL ) m = basis_stim[is]->nfunc ;
02805                   else                         m = max_lag[is] - min_lag[is] + 1;
02806               nbricks += m * (1 + option_data->tout);
02807               nbricks += option_data->rout + option_data->fout;
02808             }
02809           }
02810       }
02811 
02812       nbricks += option_data->rout + option_data->fout + option_data->vout;
02813 
02814       if (num_glt > 0)
02815       for (iglt = 0;  iglt < num_glt;  iglt++)
02816         {
02817           nbricks += glt_rows[iglt] * (1 + option_data->tout);
02818           nbricks += option_data->rout + option_data->fout;
02819         }
02820 
02821       if (nbricks <= 0)
02822       {
02823         sprintf (message,
02824                "User requested bucket dataset with only %d sub-bricks",
02825                nbricks);
02826         DC_error (message);
02827       }
02828 
02829     }
02830   option_data->nbricks = nbricks;
02831 
02832 
02833   /*----- Check for zero slice offsets with nptr option -----*/
02834   if (dset_time != NULL)
02835     if (dset_time->taxis->nsl > 0)
02836       for (is = 0;  is < num_stimts;  is++)
02837       if (nptr[is] > 1)
02838         {
02839           sprintf (message, "Must align all slices to 0 offset time, \n ");
02840           strcat  (message, "before using -stim_nptr option.  ");
02841           strcat  (message, "See program 3dTshift. ");
02842           DC_error (message);
02843         }
02844 
02845 
02846   /*----- Check for -tshift and -input1D option -----*/
02847   if ((option_data->tshift) && (option_data->input1D_filename != NULL))
02848     DC_error ("-tshift option is not compatible with -input1D option");
02849 
02850 
02851   /*----- Check whether any of the output files already exist -----*/
02852   if (option_data->input_filename != NULL)
02853     check_output_files (option_data, dset_time);
02854 
02855 }

02508 {
02509   char message[THD_MAX_NAME];      /* error message */
02510   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
02511   int ierror;                         /* number of errors in editing data */
02512 
02513 
02514   /*----- make an empty copy of input dataset -----*/
02515   new_dset = EDIT_empty_copy( dset_time ) ;
02516 
02517 
02518   ierror = EDIT_dset_items( new_dset ,
02519                       ADN_prefix , filename ,
02520                       ADN_label1 , filename ,
02521                       ADN_self_name , filename ,
02522                       ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE :
02523                                                             GEN_FUNC_TYPE ,
02524                       ADN_none ) ;
02525 
02526   if( ierror > 0 )
02527     {
02528       sprintf (message,
02529              "*** %d errors in attempting to create output dataset!\n",
02530              ierror);
02531       DC_error (message);
02532     }
02533 
02534   if( THD_is_file(new_dset->dblk->diskptr->header_name) )
02535     {
02536       sprintf (message,
02537              "Output dataset file %s already exists "
02538              " -- cannot continue! ",
02539              new_dset->dblk->diskptr->header_name);
02540       DC_error (message);
02541     }
02542 
02543   /*----- deallocate memory -----*/
02544   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
02545 
02546 }

02560 {
02561   int is;                         /* stimulus time series index */
02562 
02563 
02564   if (option_data->bucket_filename != NULL)
02565     check_one_output_file (dset_time, option_data->bucket_filename);
02566 
02567   if (option_data->fitts_filename != NULL)
02568     check_one_output_file (dset_time, option_data->fitts_filename);
02569 
02570   if (option_data->errts_filename != NULL)
02571     check_one_output_file (dset_time, option_data->errts_filename);
02572 
02573   for (is = 0;  is < option_data->num_stimts;  is++)
02574     {
02575       if (option_data->iresp_filename[is] != NULL)
02576       check_one_output_file (dset_time, option_data->iresp_filename[is]);
02577 
02578       if (option_data->sresp_filename[is] != NULL)
02579       check_one_output_file (dset_time, option_data->sresp_filename[is]);
02580     }
02581 }

06326 {
06327    int nerr = 0 ;
06328 
06329    if( X.rows < 1 || X.cols < 1 ){
06330      fprintf(stderr,
06331              "** ERROR: -xrestore %s has bad X matrix\n",xrestore_filename) ;
06332      nerr++ ;
06333    }
06334    if( XtXinv.rows < 1 || XtXinv.cols < 1 ){
06335      fprintf(stderr,
06336              "** ERROR: -xrestore %s has bad XtXinv matrix\n",xrestore_filename) ;
06337      nerr++ ;
06338    }
06339    if( XtXinvXt.rows < 1 || XtXinvXt.cols < 1 ){
06340      fprintf(stderr,
06341              "** ERROR: -xrestore %s has bad XtXinvXt matrix\n",xrestore_filename) ;
06342      nerr++ ;
06343    }
06344    if( nParam > 0 && X.cols != nParam ){
06345      fprintf(stderr,
06346              "** ERROR: -xrestore %s X matrix cols mismatch: %d != %d\n",
06347              xrestore_filename, X.cols , nParam ) ;
06348      nerr++ ;
06349    }
06350    if( GoodList == NULL ){
06351      fprintf(stderr,
06352              "** ERROR: -xrestore %s missing GoodList field\n",xrestore_filename) ;
06353      nerr++ ;
06354    } else if( nGoodList != X.rows ){
06355      fprintf(stderr,
06356              "** ERROR: -xrestore %s X matrix rows mismatch: %d != %d\n",
06357              xrestore_filename, X.cols , nGoodList ) ;
06358      nerr++ ;
06359    }
06360 #if 0
06361    if( ParamStim == NULL ){
06362      fprintf(stderr,
06363              "** ERROR: -xrestore %s missing ParamStim field\n",xrestore_filename) ;
06364      nerr++ ;
06365    }
06366    if( ParamLabel == NULL ){
06367      fprintf(stderr,
06368              "** ERROR: -xrestore %s missing ParamLabel field\n",xrestore_filename) ;
06369      nerr++ ;
06370    }
06371 #endif
06372    if( InputFilename == NULL ){
06373      fprintf(stderr,
06374              "** ERROR: -xrestore %s missing InputFilename field\n",xrestore_filename) ;
06375      nerr++ ;
06376    }
06377 
06378    if( nerr > 0 ) exit(1) ;   /** bad bad bad **/
06379    return ;
06380 }

04379 {
04380   THD_3dim_dataset * dset = NULL;        /* input afni data set pointer */
04381   int nx, ny, nz, nxyz;  /* dataset dimensions in voxels */
04382   int ixyz;              /* voxel index */
04383   int isl;               /* slice index */
04384   int i;                 /* data point index */
04385   float * yarray = NULL; /* impulse response function for a single voxel */
04386   float * sarray = NULL; /* second derivative of the cubic in each interval */
04387   matrix m, minv;        /* matrices for cubic spline interpolation */
04388   vector v, sv;          /* vectors for cubic spline interpolation */
04389   int n;                 /* number of intervals = ts_length-1 */
04390   float * a = NULL,
04391         * b = NULL,
04392         * c = NULL,
04393         * d = NULL;      /* cubic spline interpolation polynomial coefs. */
04394   float tslice;          /* slice acquisition time offset */
04395   float tdelta;          /* time between same slice acquisitons */
04396   float frac;            /* fraction of interval for slice acq. time offset */
04397   int k;                 /* interval to use for interpolation */
04398   float t;               /* time in fractions of TR */
04399   float delt;            /* time offset relative to interpolation interval */
04400   float y;               /* interpolated value */
04401 
04402 
04403   /*----- Initialize matrices and vectors -----*/
04404   matrix_initialize (&m);
04405   matrix_initialize (&minv);
04406   vector_initialize (&v);
04407   vector_initialize (&sv);
04408 
04409 
04410   /*----- Initialize local variables -----*/
04411   dset = THD_open_dataset (option_data->input_filename);
04412   n = ts_length - 1;
04413   tdelta = dset->taxis->ttdel;
04414   nx = dset->daxes->nxx;   ny = dset->daxes->nyy;   nz = dset->daxes->nzz;
04415   nxyz = nx * ny * nz;
04416 
04417 
04418   /*----- Allocate space for data and interpolation polynomials -----*/
04419   yarray = (float *) malloc (sizeof(float) * ts_length);
04420   sarray = (float *) malloc (sizeof(float) * (n+1));
04421   a = (float *) malloc (sizeof(float) * n);
04422   b = (float *) malloc (sizeof(float) * n);
04423   c = (float *) malloc (sizeof(float) * n);
04424   d = (float *) malloc (sizeof(float) * n);
04425 
04426 
04427   /*----- Calculate matrix for cubic spline interpolation -----*/
04428   matrix_create (n-1, n-1, &m);
04429   m.elts[0][0] = 4.0;
04430   m.elts[0][1] = 1.0;
04431   m.elts[n-2][n-3] = 1.0;
04432   m.elts[n-2][n-2] = 4.0;
04433   for (i = 1;  i < n-2;  i++)
04434     {
04435       m.elts[i][i] = 4.0;
04436       m.elts[i][i-1] = 1.0;
04437       m.elts[i][i+1] = 1.0;
04438     }
04439   matrix_inverse (m, &minv);
04440 
04441 
04442   vector_create (n-1, &v);
04443 
04444 
04445   /*----- Loop over all voxels -----*/
04446   for (ixyz = 0;  ixyz < nxyz;  ixyz++)
04447     {
04448 
04449       /*----- Get time offset for this slice -----*/
04450       isl = ixyz / (nx*ny);
04451       tslice = THD_timeof_slice (0, isl, dset);
04452       frac = -((tslice/tdelta) - 0.5);
04453 
04454       /*----- Get impulse response function for this voxel -----*/
04455       for (i = 0;  i < ts_length;  i++)
04456       yarray[i] = vol_array[i][ixyz];
04457 
04458 
04459       /*----- Calculate vector for cubic spline interpolation -----*/
04460       for (i = 1;  i < n;  i++)
04461            v.elts[i-1] = 6.0 * (yarray[i+1] - 2.0 * yarray[i] + yarray[i-1]);
04462       vector_multiply (minv, v, &sv);
04463 
04464 
04465       /*----- Set array of second derivatives -----*/
04466       for (i = 1;  i < n;  i++)
04467       {
04468         sarray[i] = sv.elts[i-1];
04469       }
04470       sarray[0] = 0.0;
04471       sarray[n] = 0.0;
04472 
04473 
04474       /*----- Calculate cubic spline polynomial coefficients -----*/
04475       for (i = 0;  i < n;  i++)
04476       {
04477         a[i] = (sarray[i+1] - sarray[i]) / 6.0;
04478         b[i] = sarray[i] / 2;
04479         c[i] = (yarray[i+1]-yarray[i]) - (2.0*sarray[i]+sarray[i+1]) / 6.0;
04480         d[i] = yarray[i];
04481       }
04482 
04483 
04484       /*----- Apply time shift to impulse response function -----*/
04485       for (i = 0;  i < ts_length;  i++)
04486       {
04487         t = i + frac;
04488 
04489         if (frac < 0.0)  k = i-1;
04490         else             k = i;
04491 
04492         if (k < 0)    k = 0;
04493         if (k > n-1)  k = n-1;
04494 
04495         delt = t - k;
04496 
04497         yarray[i] = a[k]*delt*delt*delt + b[k]*delt*delt + c[k]*delt + d[k];
04498       }
04499 
04500 
04501       /*----- Save interpolated impulse response function -----*/
04502       for (i = 0;  i < ts_length;  i++)
04503       vol_array[i][ixyz] = yarray[i];
04504 
04505 
04506 
04507     }  /* Loop over voxels */
04508 
04509 
04510   /*----- Deallocate memory -----*/
04511   THD_delete_3dim_dataset (dset, False);  dset = NULL ;
04512 
04513   matrix_destroy (&m);
04514   matrix_destroy (&minv);
04515   vector_destroy (&v);
04516   vector_destroy (&sv);
04517 
04518   free (sarray);   sarray = NULL;
04519   free (yarray);   yarray = NULL;
04520   free (a);        a = NULL;
04521   free (b);        b = NULL;
04522   free (c);        c = NULL;
04523   free (d);        d = NULL;
04524 }

00543 {
00544   fprintf (stderr, "%s Error: %s \a\n\n", PROGRAM_NAME, message);
00545   exit(1);
00546 }

00580 {
00581   identify_software();
00582 
00583   printf (
00584 "Program to calculate the deconvolution of a measurement 3d+time dataset    \n"
00585 "with a specified input stimulus time series.  This program will also       \n"
00586 "perform multiple linear regression using multiple input stimulus time      \n"
00587 "series. Output consists of an AFNI 'bucket' type dataset containing the    \n"
00588 "least squares estimates of the linear regression coefficients, t-statistics\n"
00589 "for significance of the coefficients, partial F-statistics for significance\n"
00590 "of the individual input stimuli, and the F-statistic for significance of   \n"
00591 "the overall regression.  Additional output consists of a 3d+time dataset   \n"
00592 "containing the estimated system impulse response function.                 \n"
00593     "                                                                       \n"
00594     "Usage:                                                                 \n"
00595     PROGRAM_NAME "\n"
00596     "                                                                       \n"
00597     "**** Input data and control options:                                   \n"
00598     "-input fname         fname = filename of 3d+time input dataset         \n"
00599     "                       (more than  one filename  can be  given)        \n"
00600     "                       (here,   and  these  datasets  will  be)        \n"
00601     "                       (catenated  in time;   if you do this, )        \n"
00602     "                       ('-concat' is not needed and is ignored)        \n"
00603     "[-input1D dname]     dname = filename of single (fMRI) .1D time series \n"
00604     "[-nodata [NT [TR]]   Evaluate experimental design only (no input data) \n"
00605     "[-mask mname]        mname = filename of 3d mask dataset               \n"
00606     "[-automask]          build a mask automatically from input data        \n"
00607     "                      (will be slow for long time series datasets)     \n"
00608     "[-censor cname]      cname = filename of censor .1D time series        \n"
00609     "[-concat rname]      rname = filename for list of concatenated runs    \n"
00610     "[-nfirst fnum]       fnum = number of first dataset image to use in the\n"
00611     "                       deconvolution procedure. (default = max maxlag) \n"
00612     "[-nlast  lnum]       lnum = number of last dataset image to use in the \n"
00613     "                       deconvolution procedure. (default = last point) \n"
00614     "[-polort pnum]       pnum = degree of polynomial corresponding to the  \n"
00615     "                       null hypothesis  (default: pnum = 1)            \n"
00616     "[-legendre]          use Legendre polynomials for null hypothesis      \n"
00617     "[-nolegendre]        use power polynomials for null hypotheses         \n"
00618     "                       (default is -legendre)                          \n"
00619     "[-nodmbase]          don't de-mean baseline time series                \n"
00620     "                       (i.e., polort>1 and -stim_base inputs)          \n"
00621     "[-dmbase]            de-mean baseline time series (default if polort>0)\n"
00622     "[-nocond]            don't calculate matrix condition number           \n"
00623     "[-svd]               Use SVD instead of Gaussian elimination (default) \n"
00624     "[-nosvd]             Use Gaussian elimination instead of SVD           \n"
00625     "[-rmsmin r]          r = minimum rms error to reject reduced model     \n"
00626     "                                                                       \n"
00627     "**** Input stimulus options:                                           \n"
00628     "-num_stimts num      num = number of input stimulus time series        \n"
00629     "                       (0 <= num)   (default: num = 0)                 \n"
00630     "-stim_file k sname   sname = filename of kth time series input stimulus\n"
00631     "[-stim_label k slabel] slabel = label for kth input stimulus           \n"
00632     "[-stim_base k]       kth input stimulus is part of the baseline model  \n"
00633     "[-stim_minlag k m]   m = minimum time lag for kth input stimulus       \n"
00634     "                       (default: m = 0)                                \n"
00635     "[-stim_maxlag k n]   n = maximum time lag for kth input stimulus       \n"
00636     "                       (default: n = 0)                                \n"
00637     "[-stim_nptr k p]     p = number of stimulus function points per TR     \n"
00638     "                       Note: This option requires 0 slice offset times \n"
00639     "                       (default: p = 1)                                \n"
00640     "                                                                       \n"
00641     "[-stim_times k tname Rmodel]                                           \n"
00642     "   Generate the k-th response model from a set of stimulus times       \n"
00643     "   given in file 'tname'.  The response model is specified by the      \n"
00644     "   'Rmodel' argument, which can be one of                              \n"
00645     "     'GAM(p,q)'    = 1 parameter gamma variate                         \n"
00646     "     'SPMG'        = 2 parameter SPM gamma variate + derivative        \n"
00647     "     'POLY(b,c,n)' = n parameter polynomial expansion                  \n"
00648     "     'SIN(b,c,n)'  = n parameter sine series expansion                 \n"
00649     "     'TENT(b,c,n)' = n parameter tent function expansion               \n"
00650     "     'BLOCK(d,p)'  = 1 parameter block stimulus of duration 'd'        \n"
00651     "                     (can also be called 'IGFUN' which stands)         \n"
00652     "                     (for 'incomplete gamma function'        )         \n"
00653     "     'EXPR(b,c) exp1 ... expn' = n parameter; arbitrary expressions    \n"
00654     "                                                                       \n"
00655     "[-basis_normall a]                                                     \n"
00656     "   Normalize all basis functions for '-stim_times' to have             \n"
00657     "   amplitude 'a' (must have a > 0).  The peak absolute value           \n"
00658     "   of each basis function will be scaled to be 'a'.                    \n"
00659     "   NOTE: -basis_normall only affect -stim_times options that           \n"
00660     "         appear LATER on the command line                              \n"
00661     "                                                                       \n"
00662     "**** General linear test (GLT) options:                                \n"
00663     "-num_glt num         num = number of general linear tests (GLTs)       \n"
00664     "                       (0 <= num)   (default: num = 0)                 \n"
00665     "[-glt s gltname]     Perform s simultaneous linear tests, as specified \n"
00666     "                       by the matrix contained in file gltname         \n"
00667     "[-glt_label k glabel]  glabel = label for kth general linear test      \n"
00668     "[-gltsym gltname]    Read the GLT with symbolic names from the file    \n"
00669     "                                                                       \n"
00670     "[-TR_irc dt]                                                           \n"
00671     "   Use 'dt' as the stepsize for computation of integrals in -IRC_times \n"
00672     "   options.  Default is to use value given in '-TR_times'.             \n"
00673     "                                                                       \n"
00674     "**** Options for output 3d+time datasets:                              \n"
00675     "[-iresp k iprefix]   iprefix = prefix of 3d+time output dataset which  \n"
00676     "                       will contain the kth estimated impulse response \n"
00677     "[-tshift]            Use cubic spline interpolation to time shift the  \n"
00678     "                       estimated impulse response function, in order to\n"
00679     "                       correct for differences in slice acquisition    \n"
00680     "                       times. Note that this effects only the 3d+time  \n"
00681     "                       output dataset generated by the -iresp option.  \n"
00682     "[-sresp k sprefix]   sprefix = prefix of 3d+time output dataset which  \n"
00683     "                       will contain the standard deviations of the     \n"
00684     "                       kth impulse response function parameters        \n"
00685     "[-fitts  fprefix]    fprefix = prefix of 3d+time output dataset which  \n"
00686     "                       will contain the (full model) time series fit   \n"
00687     "                       to the input data                               \n"
00688     "[-errts  eprefix]    eprefix = prefix of 3d+time output dataset which  \n"
00689     "                       will contain the residual error time series     \n"
00690     "                       from the full model fit to the input data       \n"
00691     "[-TR_times dt]                                                         \n"
00692     "   Use 'dt' as the stepsize for output of -iresp and -sresp file       \n"
00693     "   for response models generated by '-stim_times' options.             \n"
00694     "   Default is same as time spacing in the '-input' 3D+time dataset.    \n"
00695     "                                                                       \n"
00696     "**** Options to control the contents of the output bucket dataset:     \n"
00697     "[-fout]            Flag to output the F-statistics                     \n"
00698     "[-rout]            Flag to output the R^2 statistics                   \n"
00699     "[-tout]            Flag to output the t-statistics                     \n"
00700     "[-vout]            Flag to output the sample variance (MSE) map        \n"
00701     "[-nobout]          Flag to suppress output of baseline coefficients    \n"
00702     "                     (and associated statistics)                       \n"
00703     "[-nocout]          Flag to suppress output of regression coefficients  \n"
00704     "                     (and associated statistics)                       \n"
00705     "[-full_first]      Flag to specify that the full model statistics will \n"
00706     "                     appear first in the bucket dataset output         \n"
00707     "[-bucket bprefix]  Create one AFNI 'bucket' dataset containing various \n"
00708     "                     parameters of interest, such as the estimated IRF \n"
00709     "                     coefficients, and full model fit statistics.      \n"
00710     "                     Output 'bucket' dataset is written to bprefix.    \n"
00711     "                                                                       \n"
00712     "[-xsave]           Flag to save X matrix into file bprefix.xsave       \n"
00713     "                     (only works if -bucket option is also given)      \n"
00714     "[-noxsave]         Don't save X matrix (this is the default)           \n"
00715     "[-cbucket cprefix] Save the regression coefficients (no statistics)    \n"
00716     "                     into a dataset named 'cprefix'.  This dataset     \n"
00717     "                     will be used in a -xrestore run instead of the    \n"
00718     "                     bucket dataset, if possible.                      \n"
00719     "                                                                       \n"
00720     "[-xrestore f.xsave] Restore the X matrix, etc. from a previous run     \n"
00721     "                     that was saved into file 'f.xsave'.  You can      \n"
00722     "                     then carry out new -glt tests.  When -xrestore    \n"
00723     "                     is used, most other command line options are      \n"
00724     "                     ignored.                                          \n"
00725     "                                                                       \n"
00726     "**** The following options control the screen output only:             \n"
00727     "[-quiet]             Flag to suppress most screen output               \n"
00728     "[-xout]              Flag to write X and inv(X'X) matrices to screen   \n"
00729     "[-xjpeg filename]    Write a JPEG file graphing the X matrix           \n"
00730     "[-progress n]        Write statistical results for every nth voxel     \n"
00731     "[-fdisp fval]        Write statistical results for those voxels        \n"
00732     "                       whose full model F-statistic is > fval          \n"
00733     );
00734 
00735 #ifdef PROC_MAX
00736     printf( "\n"
00737             " -jobs J   Run the program with 'J' jobs (sub-processes).\n"
00738             "             On a multi-CPU machine, this can speed the\n"
00739             "             program up considerably.  On a single CPU\n"
00740             "             machine, using this option is silly.\n"
00741             "             J should be a number from 1 up to the\n"
00742             "             number of CPU sharing memory on the system.\n"
00743             "             J=1 is normal (single process) operation.\n"
00744             "             The maximum allowed value of J is %d.\n"
00745             "         * For more information on parallelizing, see\n"
00746             "             http://afni.nimh.nih.gov/afni/doc/misc/parallize.html \n"
00747             "         * Use -mask to get more speed; cf. 3dAutomask.\n"
00748           , PROC_MAX ) ;
00749 #endif
00750 
00751     printf("\n"
00752            "** NOTE **\n"
00753            "This version of the program has been compiled to use\n"
00754 #ifdef FLOATIZE
00755            "single precision arithmetic for most internal calculations.\n"
00756 #else
00757            "double precision arithmetic for most internal calculations.\n"
00758 #endif
00759           ) ;
00760 
00761   exit(0);
00762 }

06385 {
06386    THD_3dim_dataset *dset_time , *dset_coef, *dset_buck ;
06387    int nt , np , ii, nxyz, tout,rout,fout, ixyz,novar,view ;
06388    char *buck_prefix , *buck_name ;
06389    char brick_label[THD_MAX_NAME] ;
06390 
06391    float *ts_array = NULL; /* array of measured data for one voxel */
06392 #undef  NGET
06393 #define NGET 32            /* number to get at one time */
06394    int nget=0 ,            /* number of time series current gotten */
06395        cget=0 ,            /* index of next timeseries in iget & imget */
06396        jget   ,            /* loop index for iget */
06397        iget[NGET] ;        /* voxel index of timeseries */
06398    MRI_IMARR *imget=NULL ; /* array of timeseries */
06399 
06400    int num_glt , iglt , ilc,nlc ;
06401    matrix *glt_cmat , *glt_amat , *cxtxinvct ;
06402    vector *glt_coef , *glt_tcoef, y , coef ;
06403    float *fglt , *rglt ;
06404    char **glt_label ;
06405    int *glt_rows ;
06406    float ***glt_coef_vol       ,
06407          ***glt_tcoef_vol=NULL ,
06408          ** glt_fstat_vol=NULL ,
06409          ** glt_rstat_vol=NULL  ;
06410    float *cdar=NULL , ssef , *volume ;
06411    int ivol , nvol , nbuck , vstep ;
06412 
06413    /*----- Check for GLT options -----*/
06414 
06415    num_glt   = option_data->num_glt   ;
06416    glt_label = option_data->glt_label ;
06417    glt_rows  = option_data->glt_rows  ;
06418 
06419    tout = option_data->tout ;
06420    rout = option_data->rout ;
06421    fout = option_data->fout ;
06422 
06423    if( num_glt < 1 ){
06424      fprintf(stderr,"** ERROR: -xrestore with no new GLTs???\n"); exit(1);
06425    }
06426 
06427    /*----- initialize values to be read from xsave file -----*/
06428 
06429    matrix_initialize( &X )        ;
06430    matrix_initialize( &XtXinv )   ;
06431    matrix_initialize( &XtXinvXt ) ;
06432    nGoodList = 0 ; GoodList   = NULL ;
06433    nParam    = 0 ; ParamIndex = NULL ; ParamStim = NULL ; ParamLabel = NULL ;
06434 
06435    /*----- read xsave file -----*/
06436 
06437    if( verb ) fprintf(stderr,"++ Starting -xrestore %s\n",xrestore_filename) ;
06438 
06439    XSAVE_input( xrestore_filename ) ;
06440    check_xrestore_data() ;
06441 
06442    nt = X.rows ; np = X.cols ;  /* number of time points and parameters */
06443 
06444    /*----- read input time series dataset -----*/
06445 
06446    if( verb ) fprintf(stderr,"++ loading time series dataset '%s'\n",InputFilename);
06447    dset_time = THD_open_dataset( InputFilename ) ;
06448    if( dset_time == NULL ){
06449      fprintf(stderr,
06450              "** ERROR: -xrestore can't open time series dataset '%s'\n" ,
06451              InputFilename ) ;
06452      exit(1) ;
06453    }
06454    DSET_load( dset_time ) ;
06455    if( !DSET_LOADED(dset_time) ){
06456      fprintf(stderr,
06457              "** ERROR: -xrestore can't load time series dataset '%s'\n" ,
06458              InputFilename ) ;
06459      exit(1) ;
06460    }
06461    nxyz = DSET_NVOX(dset_time) ;  /* number of voxels */
06462    view = dset_time->view_type ;  /* +orig, +acpc, +tlrc ? */
06463 
06464    /*----- read coefficient dataset (if possible) -----*/
06465 
06466    dset_coef = NULL ;
06467    if( CoefFilename != NULL ){
06468      if( verb) fprintf(stderr,"++ loading coefficient dataset %s\n",CoefFilename);
06469      dset_coef = THD_open_dataset( CoefFilename ) ;
06470      if( dset_coef == NULL ){
06471        fprintf(stderr,
06472                "** WARNING: -xrestore can't open coefficient dataset %s\n",
06473                CoefFilename);
06474      } else {
06475        DSET_load(dset_coef) ;
06476        if( !DSET_LOADED(dset_coef) ){
06477          fprintf(stderr,
06478                  "** WARNING: -xrestore can't load coefficient dataset %s\n",
06479                  CoefFilename);
06480          DSET_delete(dset_coef) ; dset_coef = NULL ;
06481        }
06482      }
06483      if( dset_coef != NULL && DSET_NVALS(dset_coef) < np ){
06484        fprintf(stderr,
06485                "** WARNING: -xrestore coefficient dataset %s too short\n",
06486                CoefFilename);
06487        DSET_delete(dset_coef) ; dset_coef = NULL ;
06488      }
06489      if( dset_coef != NULL ){
06490        if( ParamIndex != NULL ) free((void *)ParamIndex) ;
06491        ParamIndex = (int *)malloc(sizeof(int)*np) ;
06492        for( ii=0 ; ii < np ; ii++ ) ParamIndex[ii] = ii ;
06493      }
06494    }
06495 
06496    /*----- if above failed, try the old bucket dataset -----*/
06497 
06498    if( dset_coef == NULL && BucketFilename != NULL && ParamIndex != NULL ){
06499      if( verb ) fprintf(stderr,"++ loading original bucket dataset %s\n",BucketFilename) ;
06500      dset_coef = THD_open_dataset( BucketFilename ) ;
06501      if( dset_coef == NULL ){
06502        fprintf(stderr,
06503                "** WARNING: -xrestore can't open old bucket dataset %s\n",
06504                BucketFilename);
06505      } else {
06506        DSET_load(dset_coef) ;
06507        if( !DSET_LOADED(dset_coef) ){
06508          fprintf(stderr,
06509                  "** WARNING: -xrestore can't load old bucket dataset %s\n",
06510                  BucketFilename);
06511          DSET_delete(dset_coef) ; dset_coef = NULL ;
06512        }
06513      }
06514      if( dset_coef != NULL && DSET_NVALS(dset_coef) < ParamIndex[np-1] ){
06515        fprintf(stderr,
06516                "** WARNING: -xrestore old bucket dataset %s too short\n",
06517                BucketFilename);
06518        DSET_delete(dset_coef) ; dset_coef = NULL ;
06519      }
06520    }
06521 
06522    if( ISVALID_DSET(dset_coef) && DSET_NVOX(dset_coef) != nxyz ){
06523      fprintf(stderr,
06524              "** ERROR: dataset mismatch between time series and coef!\n") ;
06525      exit(1) ;
06526    }
06527 
06528    /*----- neither worked ==> must recompute from input data time series -----*/
06529 
06530    if( dset_coef == NULL && verb )
06531     fprintf(stderr,
06532            "++ -xrestore recomputing coefficients from time series\n");
06533 
06534    /*----- read new GLT matrices -----*/
06535 
06536    glt_cmat = (matrix *) malloc( sizeof(matrix) * num_glt ) ;
06537 
06538    for( iglt=0; iglt < num_glt ; iglt++){
06539      matrix_initialize( glt_cmat + iglt ) ;
06540 #if 1
06541      read_glt_matrix( option_data->glt_filename[iglt] ,
06542                       option_data->glt_rows + iglt ,
06543                       np , glt_cmat + iglt             ) ;
06544 #else
06545      matrix_file_read ( option_data->glt_filename[iglt] ,
06546                             option_data->glt_rows[iglt] ,
06547                         np , glt_cmat+iglt , 1           ) ;
06548      if( glt_cmat[iglt].elts == NULL ){
06549        fprintf(stderr,
06550                "** ERROR: Can't read GLT matrix from file %s\n",
06551                option_data->glt_filename[iglt] ) ;
06552        exit(1) ;
06553      }
06554 #endif
06555    }
06556 
06557    /*----- initialize new GLT calculations:
06558             - setup space for matrices and vectors
06559             - calculate matrices
06560             - malloc space for output bricks -----*/
06561 
06562    cxtxinvct = (matrix *) malloc( sizeof(matrix) * num_glt ) ;
06563    glt_amat  = (matrix *) malloc( sizeof(matrix) * num_glt ) ;
06564    glt_coef  = (vector *) malloc( sizeof(vector) * num_glt ) ;
06565    glt_tcoef = (vector *) malloc( sizeof(vector) * num_glt ) ;
06566 
06567    for( iglt=0 ; iglt < num_glt ; iglt++ ){
06568      matrix_initialize( cxtxinvct+iglt ) ;  /* will be loaded in   */
06569      matrix_initialize( glt_amat +iglt ) ;  /* init_glt_analysis() */
06570      vector_initialize( glt_coef +iglt ) ;
06571      vector_initialize( glt_tcoef+iglt ) ;
06572    }
06573 
06574    fglt = (float *) malloc( sizeof(float) * num_glt ) ;
06575    rglt = (float *) malloc( sizeof(float) * num_glt ) ;
06576 
06577    vector_initialize( &y )   ; vector_create( nt, &y    );  /* time series */
06578    vector_initialize( &coef ); vector_create( np, &coef );  /* parameters */
06579 
06580    if( dset_coef != NULL )
06581      cdar = (float *)malloc( sizeof(float) * DSET_NVALS(dset_coef) ) ;
06582 
06583    init_glt_analysis( XtXinv , num_glt , glt_cmat , glt_amat , cxtxinvct ) ;
06584 
06585    /*-- malloc output bricks --*/
06586 
06587    glt_coef_vol = (float ***) malloc( sizeof(float **) * num_glt ) ;
06588    if( tout )
06589      glt_tcoef_vol  = (float ***) malloc( sizeof(float **) * num_glt ) ;
06590    if( fout )
06591      glt_fstat_vol = (float **) malloc( sizeof(float *)  * num_glt ) ;
06592    if( rout )
06593      glt_rstat_vol = (float **) malloc( sizeof(float *)  * num_glt ) ;
06594 
06595    nvol = 0 ;
06596    for( iglt=0 ; iglt < num_glt ; iglt++ ){
06597      nlc = glt_rows[iglt];
06598      glt_coef_vol[iglt] = (float **) malloc( sizeof(float *) * nlc ) ;
06599      if( tout )
06600        glt_tcoef_vol[iglt] = (float **) malloc( sizeof(float *) * nlc ) ;
06601      for( ilc=0 ; ilc < nlc ; ilc++ ){
06602        glt_coef_vol[iglt][ilc] = (float *)calloc(sizeof(float),nxyz); nvol++;
06603        if( tout ){
06604          glt_tcoef_vol[iglt][ilc] = (float *)calloc(sizeof(float),nxyz); nvol++;
06605        }
06606      }
06607      if( fout ){
06608        glt_fstat_vol[iglt] = (float *)calloc( sizeof(float), nxyz ); nvol++;
06609      }
06610      if( rout ){
06611        glt_rstat_vol[iglt] = (float *)calloc( sizeof(float), nxyz ); nvol++;
06612      }
06613    }
06614 
06615    /*----- loop over voxels:
06616             - fetch coefficients (check for all zero), or recompute them
06617             - fetch time series
06618             - compute SSE of full model
06619             - compute and store new GLT results in arrays -----*/
06620 
06621    vstep = nxyz / 50 ; if( !verb ) vstep = 0 ;
06622    if( vstep > 0 ) fprintf(stderr,"++ voxel loop:") ;
06623    for( ixyz=0 ; ixyz < nxyz ; ixyz++ ){
06624 
06625      if( vstep > 0 && ixyz%vstep == vstep-1 ) vstep_print() ;
06626 
06627      /*** race ahead and extract a bunch of voxel time series at once ***/
06628 
06629      if( cget == nget ){
06630        if( imget != NULL ) DESTROY_IMARR(imget) ;
06631        iget[0] = ixyz ; nget = 1 ;
06632        for( jget=ixyz+1 ; jget < nxyz && nget < NGET ; jget++ )
06633          iget[nget++] = jget ;
06634        imget = THD_extract_many_series( nget, iget, dset_time ) ;
06635        cget  = 0 ;  /* the next one to take out of imget */
06636      }
06637 
06638      /*** Extract Y-data for this voxel ***/
06639 
06640      ts_array = MRI_FLOAT_PTR(IMARR_SUBIM(imget,cget)) ; cget++ ;
06641      for( ii=0 ; ii < nt ; ii++ ) y.elts[ii] = ts_array[GoodList[ii]];
06642 
06643      /*** Extract or recompute coefficients for this voxel ***/
06644 
06645      if( dset_coef != NULL ){
06646        (void) THD_extract_array( ixyz , dset_coef , 0 , cdar ) ;
06647        for( ii=0 ; ii < np ; ii++ ) coef.elts[ii] = cdar[ParamIndex[ii]] ;
06648      } else {
06649        vector_multiply( XtXinvXt , y , &coef ) ;
06650      }
06651 
06652      /*** if coef is all zero, skip this voxel ***/
06653 
06654      for( ii=0 ; ii < np ; ii++ ) if( coef.elts[ii] != 0.0 ) break ;
06655      novar = (ii==np) ;
06656 
06657      if( !novar ) ssef = calc_sse( X , coef , y ) ;
06658 
06659      /************************************/
06660      /*** Do the work we came here for ***/
06661      /************************************/
06662 
06663      glt_analysis( nt , np ,
06664                    X , y , ssef , coef , novar ,
06665                    cxtxinvct , num_glt , glt_rows , glt_cmat , glt_amat ,
06666                    glt_coef , glt_tcoef , fglt , rglt ) ;
06667 
06668      /*** save results into output bricks ***/
06669 
06670      if (glt_coef_vol != NULL)
06671       for (iglt = 0;  iglt < num_glt;  iglt++)
06672        if (glt_coef_vol[iglt] != NULL)
06673         for (ilc = 0;  ilc < glt_rows[iglt];  ilc++)
06674          glt_coef_vol[iglt][ilc][ixyz] = glt_coef[iglt].elts[ilc];
06675 
06676      if (glt_tcoef_vol != NULL)
06677       for (iglt = 0;  iglt < num_glt;  iglt++)
06678        if (glt_tcoef_vol[iglt] != NULL)
06679         for (ilc = 0;  ilc < glt_rows[iglt];  ilc++)
06680          glt_tcoef_vol[iglt][ilc][ixyz] = glt_tcoef[iglt].elts[ilc];
06681 
06682      if (glt_fstat_vol != NULL)
06683       for (iglt = 0;  iglt < num_glt;  iglt++)
06684        if (glt_fstat_vol[iglt] != NULL)
06685         glt_fstat_vol[iglt][ixyz] = fglt[iglt];
06686 
06687      if (glt_rstat_vol != NULL)
06688       for (iglt = 0;  iglt < num_glt;  iglt++)
06689        if (glt_rstat_vol[iglt] != NULL)
06690         glt_rstat_vol[iglt][ixyz] = rglt[iglt];
06691 
06692    } /*** end of loop over voxels */
06693 
06694    if( vstep > 0 ) fprintf(stderr,"\n") ;  /* end of progress meter */
06695 
06696    /*** unload input datasets to save memory ***/
06697 
06698                            DSET_unload( dset_time ) ;
06699    if( dset_coef != NULL ) DSET_delete( dset_coef ) ;
06700 
06701    /*----- open old dataset for output if
06702              (a) -bucket was given for an existing dataset, or
06703              (b) no -bucket option was given;
06704             otherwise, open a new dataset for output of the GLT results -----*/
06705 
06706    if( option_data->bucket_filename != NULL ){
06707      buck_prefix = strdup(option_data->bucket_filename) ;
06708    } else if( BucketFilename != NULL ){
06709      buck_prefix = (char *)malloc(strlen(BucketFilename)+4) ;
06710      FILENAME_TO_PREFIX( BucketFilename , buck_prefix ) ;
06711    } else {
06712      buck_prefix = strdup("X") ;   /* bad user, bad bad bad */
06713    }
06714 
06715    /*** try to open dataset as an old one **/
06716 
06717    buck_name = (char *)malloc(strlen(buck_prefix)+32) ;
06718    strcpy(buck_name,buck_prefix) ;
06719    strcat(buck_name,"+") ; strcat(buck_name,VIEW_codestr[view]) ;
06720 
06721    dset_buck = THD_open_dataset( buck_name ) ;
06722 
06723    if( dset_buck != NULL ){
06724 
06725      if( verb) fprintf(stderr,"++ -xrestore appending to dataset %s\n",buck_name) ;
06726      DSET_mallocize( dset_buck ) ;
06727      if( DSET_NVOX(dset_buck) != nxyz ){
06728        fprintf(stderr,
06729                "** ERROR: dataset %s mismatch with time series '%s'\n" ,
06730                buck_name , DSET_HEADNAME(dset_time)       );
06731        exit(0) ;
06732      }
06733      DSET_load( dset_buck ) ;
06734      if( !DSET_LOADED(dset_buck) ){
06735        fprintf(stderr,
06736                "** ERROR: can't load dataset %s from disk\n" , buck_name ) ;
06737        exit(1) ;
06738      }
06739 
06740      /* add nvol empty bricks at the end */
06741 
06742      ivol = DSET_NVALS(dset_buck) ;   /* where to save first new brick */
06743 
06744      EDIT_add_bricklist( dset_buck , nvol, NULL, NULL, NULL ) ;
06745 
06746    } else {  /*** create a new dataset ***/
06747 
06748      if( verb ) fprintf(stderr,"++ -xrestore creating new dataset %s\n",buck_name) ;
06749 
06750      dset_buck = EDIT_empty_copy( dset_time ) ;
06751      (void) EDIT_dset_items( dset_buck ,
06752                              ADN_prefix,          buck_prefix ,
06753                              ADN_type,            HEAD_FUNC_TYPE,
06754                              ADN_func_type,       FUNC_BUCK_TYPE,
06755                              ADN_datum_all,       MRI_short ,
06756                              ADN_ntt,             0,               /* no time */
06757                              ADN_nvals,           nvol ,
06758                              ADN_malloc_type,     DATABLOCK_MEM_MALLOC ,
06759                              ADN_none ) ;
06760      ivol = 0 ;
06761    }
06762 
06763    tross_Make_History( PROGRAM_NAME , argc , argv , dset_buck ) ;
06764 
06765    nbuck = DSET_NVALS(dset_buck) ;
06766 
06767    /*** attach sub-bricks to output ***/
06768 
06769    for( iglt=0 ; iglt < num_glt ; iglt++ ){
06770 
06771      for( ilc=0 ; ilc < glt_rows[iglt] ; ilc++ ){
06772        sprintf( brick_label , "%s LC[%d] coef" , glt_label[iglt], ilc ) ;
06773        volume = glt_coef_vol[iglt][ilc];
06774        attach_sub_brick( dset_buck, ivol, volume, nxyz,
06775                          FUNC_FIM_TYPE, brick_label, 0, 0, 0, NULL);
06776        free((void *)volume) ; ivol++ ;
06777 
06778        if( tout ){
06779          sprintf( brick_label , "%s LC[%d] t-st" , glt_label[iglt], ilc ) ;
06780          volume = glt_tcoef_vol[iglt][ilc];
06781          attach_sub_brick( dset_buck, ivol, volume, nxyz,
06782                            FUNC_TT_TYPE, brick_label, nt-np, 0, 0, NULL);
06783          free((void *)volume) ; ivol++ ;
06784        }
06785      }
06786 
06787      if( rout ){
06788        sprintf( brick_label , "%s R^2" , glt_label[iglt] ) ;
06789        volume = glt_rstat_vol[iglt];
06790        attach_sub_brick( dset_buck, ivol, volume, nxyz,
06791                          FUNC_THR_TYPE, brick_label, 0, 0, 0, NULL);
06792        free((void *)volume) ; ivol++ ;
06793      }
06794 
06795      if( fout ){
06796        sprintf( brick_label , "%s F-stat" , glt_label[iglt] ) ;
06797        volume = glt_fstat_vol[iglt];
06798        attach_sub_brick( dset_buck, ivol, volume, nxyz,
06799                          FUNC_FT_TYPE, brick_label, 0, glt_rows[iglt], nt-np, NULL);
06800        free((void *)volume) ; ivol++ ;
06801      }
06802    }  /** End loop over general linear tests **/
06803 
06804    /*----- Write dataset out! -----*/
06805 
06806    DSET_write( dset_buck ) ;
06807    return ;
06808 }