00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*
00008   This file contains routines used by programs 3dNLfim, plug_nlfit, and
00009   3dTSgen for performing non-linear regression analysis of AFNI 3d+time 
00010   data seta, and investigating various statistical properties of the data. 
00011 
00012   File:     NLfit.c
00013   Author:   B. Douglas Ward
00014   Date:     19 June 1997
00015 
00016   Mod:      Added option for absolute noise parameter constraints.
00017   Date:     21 August 1997
00018 
00019   Mod:      Changed random number initialization from srand to srand48.
00020   Date:     29 August 1997
00021 
00022   Mod:      Added options for percent signal change above baseline, and
00023             calculation of area under the signal above baseline.
00024   Date:     26 November 1997
00025 
00026   Mod:      Continuation of previous modification.
00027   Date:     8 January 1998
00028 
00029   Mod:      Added novar flag to eliminate unnecessary calculations.
00030   Date:     13 July 1999
00031 
00032   Mod:      Adjust F-statistics if parameter constraints force a parameter
00033             to be a constant.
00034   Date:     08 February 2000
00035 
00036   Mod:      Changes for output of R^2 (coefficient of multiple determination),
00037             and standard deviation of residuals from full model fit.
00038             Added global variable calc_tstats.
00039             Also, added screen display of p-values.
00040   Date:     10 May 2000
00041 
00042   Mod:      Add stuff for longjmp() return from NLfit_error().
00043   Date:     01 May 2003 - RWCox
00044 
00045 */
00046 
00047 /*---------------------------------------------------------------------------*/
00048 
00049 #include "NLfit_model.c"
00050 
00051 /*---------------------------------------------------------------------------*/
00052 /*
00053    Routine to print error message and stop.
00054 */
00055 
00056 #include <setjmp.h>                    /* 01 May 2003 */
00057 static int jump_on_NLfit_error = 0 ;
00058 static jmp_buf NLfit_error_jmpbuf ;
00059 
00060 void NLfit_error
00061 (
00062   char * message         /* message to be displayed */
00063 )
00064 
00065 {
00066    fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00067    if( jump_on_NLfit_error ) longjmp( NLfit_error_jmpbuf , 1 ) ;  /* 01 May 2003 */
00068    exit(1);
00069 }
00070 
00071 
00072 /*---------------------------------------------------------------------------*/
00073 /*
00074   Routine to initialize the signal model by defining the number of parameters
00075   in the signal model, the default values for the minimum and maximum 
00076   parameter constraints, and setting the pointer to the function which
00077   implements the signal model.
00078 */
00079 
00080 void initialize_signal_model 
00081 (
00082   NLFIT_MODEL_array * model_array,       /* the array of SO models */
00083   char * sname,            /* name of signal model (user input) */
00084   vfp * smodel,            /* pointer to signal model */
00085   int * p,                 /* number of parameters in the signal model */
00086   char ** spname,          /* signal parameter names */
00087   float * min_sconstr,     /* minimum parameter constraints for signal model */
00088   float * max_sconstr      /* maximum parameter constraints for signal model */
00089 )
00090 
00091 {
00092   int im, ip, index;
00093   char message[MAX_NAME_LENGTH];    /* error message */
00094 
00095 
00096   index = -1;
00097   for (im = 0;  im < model_array->num;  im++)
00098     if (strncmp(model_array->modar[im]->interface->label, 
00099                 sname, MAX_NAME_LENGTH) == 0)  index = im;
00100   if (index < 0)
00101     {
00102       sprintf (message, "Unable to locate signal model %s", sname);
00103       NLfit_error (message);
00104     }
00105 
00106   if (model_array->modar[index]->interface->model_type != MODEL_SIGNAL_TYPE)
00107     {
00108       sprintf (message, "%s has not been declared a signal model", sname);
00109       NLfit_error (message);
00110     }
00111 
00112   *smodel = model_array->modar[index]->interface->call_func;
00113   if (*smodel == NULL) 
00114     {
00115       sprintf (message, "Signal model %s not properly implemented", sname);
00116       NLfit_error (message);
00117     }
00118 
00119   *p = model_array->modar[index]->interface->params;
00120   if ((*p < 0) || (*p > MAX_PARAMETERS)) 
00121     {
00122       sprintf (message, "Illegal number of parameters for signal model %s", 
00123                sname);
00124       NLfit_error (message);
00125     }
00126 
00127   for (ip = 0;  ip < *p;  ip++)
00128     {
00129       strncpy (spname[ip], model_array->modar[index]->interface->plabel[ip],
00130                MAX_NAME_LENGTH);
00131       min_sconstr[ip] = model_array->modar[index]->interface->min_constr[ip]; 
00132       max_sconstr[ip] = model_array->modar[index]->interface->max_constr[ip]; 
00133       if (min_sconstr[ip] > max_sconstr[ip])
00134         NLfit_error("Must have signal parameter min cnstrnts <= max cnstrnts");
00135     }
00136 
00137 }
00138 
00139 
00140 /*---------------------------------------------------------------------------*/
00141 /*
00142   Routine to initialize the noise model by defining the number of parameters
00143   in the noise model, the default values for the minimum and maximum 
00144   parameter constraints, and setting the pointer to the function which
00145   implements the noise model.
00146 */
00147 
00148 void initialize_noise_model 
00149 (
00150   NLFIT_MODEL_array * model_array,       /* the array of SO models */
00151   char * nname,            /* name of noise model (user input) */
00152   vfp * nmodel,            /* pointer to noise model */
00153   int * r,                 /* number of parameters in the noise model */
00154   char ** npname,          /* noise parameter names */
00155   float * min_nconstr,     /* minimum parameter constraints for noise model */
00156   float * max_nconstr      /* maximum parameter constraints for noise model */
00157 )
00158 
00159 {
00160   int im, ip, index;
00161   char message[MAX_NAME_LENGTH];    /* error message */
00162 
00163 
00164   index = -1;
00165   for (im = 0;  im < model_array->num;  im++)
00166     if (strncmp(model_array->modar[im]->interface->label, 
00167                 nname, MAX_NAME_LENGTH) == 0)  index = im;
00168   if (index < 0)
00169     {
00170       sprintf (message, "Unable to locate noise model %s", nname);
00171       NLfit_error (message);
00172     }
00173 
00174   if (model_array->modar[index]->interface->model_type != MODEL_NOISE_TYPE)
00175     {
00176       printf ("type = %d \n", 
00177               model_array->modar[index]->interface->model_type);
00178       sprintf (message, "%s has not been declared a noise model", nname);
00179       NLfit_error (message);
00180     }
00181 
00182   *nmodel = model_array->modar[index]->interface->call_func;
00183   if (*nmodel == NULL) 
00184     {
00185       sprintf (message, "Noise model %s not properly implemented", nname);
00186       NLfit_error (message);
00187     }
00188 
00189   *r = model_array->modar[index]->interface->params;
00190   if ((*r < 0) || (*r > MAX_PARAMETERS)) 
00191     {
00192       sprintf (message, "Illegal number of parameters for noise model %s", 
00193                nname);
00194       NLfit_error (message);
00195     }
00196 
00197   for (ip = 0;  ip < *r;  ip++)
00198     {
00199       strncpy (npname[ip], model_array->modar[index]->interface->plabel[ip],
00200                MAX_NAME_LENGTH);
00201       min_nconstr[ip] = model_array->modar[index]->interface->min_constr[ip]; 
00202       max_nconstr[ip] = model_array->modar[index]->interface->max_constr[ip]; 
00203       if (min_nconstr[ip] > max_nconstr[ip])
00204         NLfit_error("Must have noise parameter min cnstrnts <= max cnstrnts");
00205     }
00206 
00207 }
00208 
00209 
00210 /*---------------------------------------------------------------------------*/
00211 /*
00212   Routine to calculate the least squares linear regression.
00213 */
00214 
00215 void calc_linear_regression 
00216 ( 
00217   matrix x,              /* matrix of constants */
00218   vector y,              /* observation vector */
00219   vector * b,            /* estimated regression coefficients */
00220   float * sse            /* error sum of squares */
00221 )
00222   
00223 {
00224   matrix xt, xtx, xtxinv, xtxinvxt;     /* intermediate matrix results */
00225   vector yhat;                          /* estimate of observation vector */
00226   vector e;                             /* error in fit */
00227   int ok;                               /* successful matrix inversion? */
00228 
00229 
00230   /*----- initialize matrices -----*/
00231   matrix_initialize (&xt);
00232   matrix_initialize (&xtx);
00233   matrix_initialize (&xtxinv);
00234   matrix_initialize (&xtxinvxt);
00235   vector_initialize (&yhat);
00236   vector_initialize (&e);
00237 
00238   /*----- calculate least squares regression coefficients -----*/
00239   matrix_transpose (x, &xt);
00240   matrix_multiply (xt, x, &xtx);
00241   ok = matrix_inverse (xtx, &xtxinv);
00242   if (!ok)  NLfit_error ("Unable to invert matrix");
00243   matrix_multiply (xtxinv, xt, &xtxinvxt);
00244   vector_multiply (xtxinvxt, y, b);
00245 
00246   /*----- calculate least squares fit -----*/
00247   vector_multiply (x, *b, &yhat);
00248 
00249   /*----- calculate error sum of squares -----*/
00250   vector_subtract (y, yhat, &e);
00251   *sse = vector_dot (e, e);
00252 
00253   /*----- dispose of matrices -----*/
00254   vector_destroy (&e);
00255   vector_destroy (&yhat);
00256   matrix_destroy (&xtxinvxt);
00257   matrix_destroy (&xtxinv);
00258   matrix_destroy (&xtx);
00259   matrix_destroy (&xt);
00260 
00261 }
00262 
00263 
00264 /*---------------------------------------------------------------------------*/
00265 /*
00266   Routine to calculate the least squares fit corresponding the the (linear)
00267   reduced model.
00268 */
00269 
00270 void calc_reduced_model 
00271 ( 
00272   int n,                 /* number of observations */        
00273   int r,                 /* number of parameters in reduced model */
00274   float ** x_array,      /* data matrix */
00275   float * y_array,       /* observed time series */
00276   float * b_array,       /* estimated parameters for reduced model */
00277   float * sse            /* error sum of squares */
00278 )
00279 
00280 {
00281   matrix x;              /* independent variable matrix */
00282   vector y;              /* observation vector */
00283   vector b;              /* regression coefficient vector */
00284 
00285 
00286   /*----- initialize matrices and vectors -----*/
00287   matrix_initialize (&x);
00288   vector_initialize (&y);
00289   vector_initialize (&b);
00290 
00291   /*----- convert data to matrix and vector types -----*/
00292   array_to_matrix (n, r, x_array, &x);
00293   array_to_vector (n, y_array, &y);
00294 
00295   /*----- calculate linear regression -----*/
00296   calc_linear_regression (x, y, &b, sse);
00297 
00298   /*----- convert vector to array -----*/
00299   vector_to_array (b, b_array);
00300 
00301   /*----- dispose of matrices and vectors -----*/
00302   vector_destroy (&b);
00303   vector_destroy (&y);
00304   matrix_destroy (&x);
00305 }
00306 
00307 
00308 /*---------------------------------------------------------------------------*/
00309 /*
00310   Routine to allocate memory required for calculation of the full model.
00311 */
00312 
00313 void initialize_full_model 
00314 (
00315   int dimension,            /* number of parameters in full model */
00316   int nbest,                /* keep nbest vectors from random search */
00317   float *** parameters,     /* array of parameter vectors */
00318   float ** sse              /* array of error sum of squares */
00319 )
00320 
00321 {
00322   int i;                    /* parameter vector index */
00323 
00324   *sse = (float *) malloc (sizeof(float) * nbest);   
00325   *parameters = (float **) malloc (sizeof(float *) * nbest);
00326   for (i = 0;  i < nbest;  i++)
00327     (*parameters)[i] = (float *) malloc (sizeof(float) * dimension);
00328        
00329 }
00330 
00331 
00332 /*---------------------------------------------------------------------------*/
00333 /*
00334   Routine to determine if the parameter vector violates the constraints.
00335 */
00336 
00337 int calc_constraints 
00338 (
00339   int r,                  /* number of parameters in the noise model */
00340   int p,                  /* number of parameters in the signal model */
00341   int nabs,               /* use absolute constraints for noise parameters */
00342   float * b_array,        /* estimated parameters for the reduced model */
00343   float * min_nconstr,    /* minimum parameter constraints for noise model */
00344   float * max_nconstr,    /* maximum parameter constraints for noise model */
00345   float * min_sconstr,    /* minimum parameter constraints for signal model */
00346   float * max_sconstr,    /* maximum parameter constraints for signal model */
00347   float * vertex          /* test parameter vector */
00348 )
00349 
00350 {
00351   int i;                  /* parameter index */
00352 
00353 
00354   /*----- check noise model constraints -----*/
00355   if (nabs)   /*--- absolute noise parameter constraints ---*/
00356     for (i = 0;  i < r;  i++)
00357       {
00358         if (vertex[i] < min_nconstr[i])  return (1);
00359         if (vertex[i] > max_nconstr[i])  return (1);
00360       }
00361   else        /*--- relative noise parameter constraints ---*/
00362     for (i = 0;  i < r;  i++)
00363       {
00364         if (vertex[i] < min_nconstr[i] + b_array[i])  return (1);
00365         if (vertex[i] > max_nconstr[i] + b_array[i])  return (1);
00366       }
00367 
00368   /*----- check signal model constraints -----*/
00369   for (i = 0;  i < p;  i++)
00370     {
00371       if (vertex[i+r] < min_sconstr[i])  return (1);
00372       if (vertex[i+r] > max_sconstr[i])  return (1);
00373     }
00374 
00375   return (0);
00376 }
00377 
00378 #define SAVE_RAN  /* 27 July 1998 -- RWCox */
00379 #ifdef  SAVE_RAN
00380 /*===========================================================================*/
00381 /* Data to save and reuse the random search parameters for the signal model. */
00382 /* The purpose is to avoid recomputation of the signal model over and over.  */
00383 
00384 /* Parameters of the noise model (not used at this time) */
00385 
00386 static vfp     OLD_nmodel      = NULL ;
00387 static int     OLD_r           = -1 ;
00388 static float * OLD_min_nconstr = NULL ;  /* [r] */
00389 static float * OLD_max_nconstr = NULL ;  /* [r] */
00390 static float * RAN_npar = NULL ;         /* [nrand*r]        */
00391 static float * RAN_nts  = NULL ;         /* [nrand*ts_length */
00392 
00393 /* Parameters of the signal model */
00394 
00395 static vfp     OLD_smodel      = NULL ;
00396 static int     OLD_p           = -1 ;
00397 static float * OLD_min_sconstr = NULL ;  /* [p] */
00398 static float * OLD_max_sconstr = NULL ;  /* [p] */
00399 static float * RAN_spar = NULL ;         /* [nrand*p]        */
00400 static float * RAN_sts  = NULL ;         /* [nrand*ts_length */
00401 
00402 static int     OLD_ts_length   = -1 ;
00403 static int     OLD_nrand       = -1 ;
00404 static float   OLD_x0          = -666.0 ;
00405 static float   OLD_x1          = -777.0 ;
00406 
00407 static int     RAN_sind = -1 ;  /* if >= 0, index of random signal to use */
00408 
00409 /*---------------------------------------------------------------------------*/
00410 /* Function to compare two floating point vectors of the same length         */
00411 
00412 int RAN_compare_vect( int n , float * a , float * b )
00413 {
00414    int ii ;
00415    if( n < 1 || a == NULL || b == NULL ) return 1 ;
00416    for( ii=0 ; ii < n ; ii++ ) if( a[ii] != b[ii] ) return 1 ;
00417    return 0 ;
00418 }
00419 
00420 /*---------------------------------------------------------------------------*/
00421 /* Function to initialize the random parameters and their time series.       */
00422 /* At present, the noise model is not used, due to its 'relative' nature.    */
00423 
00424 void RAN_setup
00425 (
00426   vfp nmodel,             /* pointer to noise model */
00427   vfp smodel,             /* pointer to signal model */
00428   int r,                  /* number of parameters in the noise model */
00429   int p,                  /* number of parameters in the signal model */
00430   int nabs,               /* use absolute constraints for noise parameters */
00431   float * min_nconstr,    /* minimum parameter constraints for noise model */
00432   float * max_nconstr,    /* maximum parameter constraints for noise model */
00433   float * min_sconstr,    /* minimum parameter constraints for signal model */
00434   float * max_sconstr,    /* maximum parameter constraints for signal model */
00435   int ts_length,          /* length of time series array */
00436   float ** x_array,       /* independent variable matrix */
00437   int nrand               /* number of random vectors to generate */
00438 )
00439 {
00440    int ip , ipt ;
00441    float * par ;
00442    float * ts ;
00443 
00444    /*** check if the signal model is dead on arrival ***/
00445 
00446    if( smodel == NULL ){
00447       OLD_smodel = NULL ;
00448       OLD_p      = -1 ;
00449       if( OLD_min_sconstr != NULL ){ free(OLD_min_sconstr); OLD_min_sconstr = NULL; }
00450       if( OLD_max_sconstr != NULL ){ free(OLD_max_sconstr); OLD_max_sconstr = NULL; }
00451       if( RAN_spar        != NULL ){ free(RAN_spar)       ; RAN_spar        = NULL; }
00452       if( RAN_sts         != NULL ){ free(RAN_sts )       ; RAN_sts         = NULL; }
00453       return ;
00454    }
00455 
00456    /*** check if the signal model has changed ***/
00457 
00458    if( smodel        != OLD_smodel                     ||
00459        p             != OLD_p                          ||
00460        ts_length     != OLD_ts_length                  ||
00461        nrand         != OLD_nrand                      ||
00462        x_array[0][1] != OLD_x0                         ||
00463        x_array[1][1] != OLD_x1                         ||
00464        RAN_compare_vect(p,min_sconstr,OLD_min_sconstr) ||
00465        RAN_compare_vect(p,max_sconstr,OLD_max_sconstr)   ){
00466 
00467       /* save parameters of new signal model */
00468 
00469 #if 0
00470       fprintf(stderr,"++ NLfit: initializing random signal models") ;
00471 #endif
00472 
00473       OLD_smodel    = smodel ;
00474       OLD_p         = p ;
00475       OLD_ts_length = ts_length ;
00476       OLD_nrand     = nrand ;
00477       OLD_x0        = x_array[0][1] ;
00478       OLD_x1        = x_array[1][1] ;
00479       if( OLD_min_sconstr != NULL ) free(OLD_min_sconstr) ;
00480       if( OLD_max_sconstr != NULL ) free(OLD_max_sconstr) ;
00481       OLD_min_sconstr = (float *) malloc(sizeof(float)*p) ;
00482       OLD_max_sconstr = (float *) malloc(sizeof(float)*p) ;
00483       memcpy( OLD_min_sconstr , min_sconstr , sizeof(float)*p ) ;
00484       memcpy( OLD_max_sconstr , max_sconstr , sizeof(float)*p ) ;
00485 
00486       /* create space for new random signal model vectors */
00487 
00488       if( RAN_spar != NULL ) free(RAN_spar) ;
00489       if( RAN_sts  != NULL ) free(RAN_sts ) ;
00490       RAN_spar = (float *) malloc(sizeof(float)*nrand*p) ;
00491       RAN_sts  = (float *) malloc(sizeof(float)*nrand*ts_length) ;
00492 
00493       /* compute new random signal vectors */
00494 
00495       for( ipt=0 ; ipt < nrand ; ipt++ ){
00496 
00497          par = RAN_spar + (ipt*p) ;         /* ipt-th parameter vector */
00498          ts  = RAN_sts  + (ipt*ts_length) ; /* ipt-th signal model time series */
00499 
00500          for( ip=0 ; ip < p ; ip++ )                 /* parameter vector */
00501             par[ip] = get_random_value(min_sconstr[ip], max_sconstr[ip]) ;
00502 
00503 #if 0
00504          smodel( par , ts_length , x_array , ts ) ;  /* time series vector */
00505 #else
00506          AFNI_CALL_VOID_4ARG(smodel ,
00507                              float *,par , int,ts_length,
00508                              float **,x_array , float *,ts ) ;
00509 #endif
00510       }
00511 
00512 #if 0
00513       fprintf(stderr," - done\n") ;
00514 #endif
00515    } /* end of signal model stowage */
00516 
00517    return ;
00518 }
00519 /*===========================================================================*/
00520 #endif /* SAVE_RAN */
00521 
00522 /*---------------------------------------------------------------------------*/
00523 /*
00524   Calculate the estimated time series using the full model.
00525 */
00526 
00527 void full_model 
00528 (
00529   vfp nmodel,                 /* pointer to noise model */
00530   vfp smodel,                 /* pointer to signal model */
00531   float * gn,                 /* parameters for noise model */
00532   float * gs,                 /* parameters for signal model */
00533   int ts_length,              /* length of time series data */
00534   float ** x_array,           /* independent variable matrix */
00535   float * yhat_array          /* output estimated time series */
00536 )
00537 
00538 {
00539   int it;                     /* time index */
00540   float * y_array = NULL;     /* estimated signal time series */
00541 #ifdef SAVE_RAN
00542   int use_model = ( RAN_sind < 0 || ts_length != OLD_ts_length ) ;
00543 #endif
00544 
00545   /*----- generate time series corresponding to signal model -----*/
00546 
00547 #ifdef SAVE_RAN
00548   if( use_model ){  /* don't use saved time series */
00549 #endif
00550      y_array = (float *) malloc (sizeof(float) * (ts_length));
00551      if (y_array == NULL)
00552        NLfit_error ("Unable to allocate memory for y_array");
00553 #if 0
00554      smodel (gs, ts_length, x_array, y_array);
00555 #else
00556      AFNI_CALL_VOID_4ARG(smodel ,
00557                          float *,gs , int,ts_length,
00558                          float **,x_array , float *,y_array ) ;
00559 #endif
00560 
00561 #ifdef SAVE_RAN
00562   } else            /* recall a saved time series */
00563      y_array = RAN_sts + (ts_length*RAN_sind) ;
00564 #endif
00565 
00566   /*----- generate time series corresponding to the noise model -----*/
00567 #if 0
00568   nmodel (gn, ts_length, x_array, yhat_array);
00569 #else
00570   AFNI_CALL_VOID_4ARG(nmodel ,
00571                       float *,gn , int,ts_length,
00572                       float **,x_array , float *,yhat_array ) ;
00573 #endif
00574 
00575   /*----- add signal and noise model time series -----*/
00576   for (it = 0;  it < ts_length;  it++)
00577     yhat_array[it] += y_array[it];
00578   
00579 
00580   /*----- deallocate memory -----*/
00581 #ifdef SAVE_RAN
00582   if( use_model )
00583 #endif
00584     free (y_array) ;
00585 
00586   y_array = NULL;  /* not really needed since this array is 'auto' */
00587 }
00588 
00589 
00590 /*---------------------------------------------------------------------------*/
00591 /*
00592   This routine calculates the error sum of squares corresponding to 
00593   the given vector of parameters for the full model.
00594 */
00595 
00596 float calc_sse 
00597 (
00598   vfp nmodel,             /* pointer to noise model */
00599   vfp smodel,             /* pointer to signal model */
00600   int r,                  /* number of parameters in the noise model */
00601   int p,                  /* number of parameters in the signal model */
00602   int nabs,               /* use absolute constraints for noise parameters */
00603   float * min_nconstr,    /* minimum parameter constraints for noise model */
00604   float * max_nconstr,    /* maximum parameter constraints for noise model */
00605   float * min_sconstr,    /* minimum parameter constraints for signal model */
00606   float * max_sconstr,    /* maximum parameter constraints for signal model */
00607   float * b_array,        /* estimated parameters for the reduced model */
00608   float * vertex,         /* test parameter vector */
00609   int ts_length,          /* length of time series array */
00610   float ** x_array,       /* independent variable matrix */
00611   float * ts_array        /* observed time series */
00612 )
00613 
00614 {
00615   const float BIG_NUMBER = 1.0e+10;   /* return large number if constraints
00616                                          are violated */
00617   int i;                              /* time point index */
00618   float t;                            /* time */
00619   float diff;                         /* error at individual time point */
00620   float sse;                          /* error sum of squares */
00621   float * y_array = NULL;             /* estimated time series */
00622 
00623 
00624   /*----- allocate memory for estimated time series -----*/
00625   y_array = (float *) malloc (sizeof(float) * ts_length);
00626 
00627   /*----- apply constraints? -----*/
00628   if (calc_constraints (r, p, nabs, b_array, min_nconstr, max_nconstr, 
00629                         min_sconstr, max_sconstr, vertex)) 
00630     { 
00631       free (y_array);   y_array = NULL;
00632       return (BIG_NUMBER);
00633     }    
00634   
00635   /*----- create estimated time series using the full model parameters -----*/
00636   full_model (nmodel, smodel, vertex, vertex + r, ts_length, x_array, y_array);
00637 
00638   /*----- calculate error sum of squares -----*/
00639   sse = 0.0;
00640   for (i = 0;  i < ts_length;  i++)
00641     {
00642       diff = ts_array[i] - y_array[i];
00643       sse += diff * diff;
00644     }
00645 
00646   /*----- release memory -----*/
00647   free (y_array);   y_array = NULL;
00648 
00649   /*----- return error sum of squares -----*/
00650   return (sse);
00651 }
00652 
00653 
00654 /*---------------------------------------------------------------------------*/
00655 /*
00656   Select and evaluate randomly chosen vectors in the parameter space.
00657 */
00658 
00659 void random_search 
00660 (
00661   vfp nmodel,             /* pointer to noise model */
00662   vfp smodel,             /* pointer to signal model */
00663   int r,                  /* number of parameters in the noise model */
00664   int p,                  /* number of parameters in the signal model */
00665   int nabs,               /* use absolute constraints for noise parameters */
00666   float * min_nconstr,    /* minimum parameter constraints for noise model */
00667   float * max_nconstr,    /* maximum parameter constraints for noise model */
00668   float * min_sconstr,    /* minimum parameter constraints for signal model */
00669   float * max_sconstr,    /* maximum parameter constraints for signal model */
00670   int ts_length,          /* length of time series array */
00671   float ** x_array,       /* independent variable matrix */
00672   float * ts_array,       /* observed time series */
00673   float * par_rdcd,       /* estimated parameters for the reduced model */
00674   int nrand,              /* number of random vectors to generate */
00675   int nbest,              /* number of random vectors to keep */
00676   float ** parameters,    /* array of best random vectors */
00677   float * response        /* array of best sse's */
00678 )
00679 
00680 {
00681   int ip;                 /* parameter index */
00682   int iv, ipt;            /* vector indices */
00683   float sse;              /* error sum of squares */
00684   float * par = NULL;     /* test parameter vector */
00685 
00686 
00687 #ifdef SAVE_RAN
00688   /*** 28 July 1998: set up to choose the best noise model for each case,
00689                      rather than a completely random one -- RWCox.       ***/
00690 #undef  BEST_NOISE
00691 #ifdef  BEST_NOISE
00692    float * qts = (float *) malloc(sizeof(float)*ts_length) ;
00693    float junk ;
00694 #endif
00695 #endif /* SAVE_RAN */
00696 
00697   /*** 27 July 1998: generate some random signal parameters,
00698                      but only if the signal model is new this time. ***/
00699 #ifdef SAVE_RAN
00700    RAN_setup( nmodel , smodel , r , p , nabs ,
00701               min_nconstr, max_nconstr,
00702               min_sconstr, max_sconstr, 
00703               ts_length, x_array, nrand ) ;
00704 #endif
00705 
00706   /*----- allocate memory for test parameter vector -----*/
00707   par = (float *) malloc (sizeof(float) * (r+p));
00708 
00709   /*----- initialize response values -----*/
00710   for (iv = 0;  iv < nbest;  iv++)
00711     response[iv] = 1.0e+30;
00712 
00713 
00714   /*----- loop over random vectors -----*/
00715   for (ipt = 0;  ipt < nrand;  ipt++)
00716     {
00717 
00718       /*** 28 July 1998: get the best noise parameter model ***/
00719 #ifdef BEST_NOISE
00720       for( ip=0 ; ip < ts_length ; ip++ )  /* subtract signal model from data */
00721          qts[ip] = ts_array[ip] - RAN_sts[ts_length*ipt+ip] ;
00722       calc_reduced_model( ts_length , r , x_array , qts , par , &junk ) ;
00723 #else
00724       /*----- select random parameters -----*/
00725       if (nabs)   /*--- absolute noise parameter constraints ---*/
00726         for (ip = 0;  ip < r;  ip++)
00727           par[ip] = get_random_value (min_nconstr[ip], max_nconstr[ip]);
00728       else        /*--- relative noise parameter constraints ---*/
00729         for (ip = 0;  ip < r;  ip++)
00730           par[ip] = get_random_value (min_nconstr[ip] + par_rdcd[ip], 
00731                                       max_nconstr[ip] + par_rdcd[ip]);
00732 #endif /* BEST_NOISE */
00733 
00734       /*** 27 July 1998: get the signal parameters from the saved set ***/
00735 #ifdef SAVE_RAN
00736       for( ip=0 ; ip < p ; ip++ )
00737          par[ip+r] = RAN_spar[ipt*p+ip] ;
00738 #else
00739       for (ip = 0;  ip < p;  ip++)
00740           par[ip+r] = get_random_value (min_sconstr[ip], max_sconstr[ip]);
00741 #endif
00742 
00743 
00744       /*----- evaluate this random vector -----*/
00745 #ifdef SAVE_RAN
00746       RAN_sind = ipt ;
00747 #endif
00748 
00749       sse = calc_sse (nmodel, smodel, r, p, nabs, min_nconstr, max_nconstr, 
00750                       min_sconstr, max_sconstr, par_rdcd, par, 
00751                       ts_length, x_array, ts_array);
00752 
00753       /*----- save best random vectors -----*/
00754 
00755       /*** The previous code for this is erroneous, since it does not
00756            allow for the case where an in-between best and worst vector
00757            is found.  For example, suppose the response array is now
00758            [ 1 3 5 ] (with nbest==3).  Then we get sse=2.  This would
00759            replace the 3 value, leaving us with [1 2 5], instead of
00760            [ 1 2 3 ], which is a better set.  The solution is to push
00761            the displaced vectors up in the array.  In this way, the
00762            response array is always sorted, and keeps the truly best. ***/
00763 
00764 #define PUSH_BEST  /* 27 July 1998 -- RWCox */
00765 #ifdef  PUSH_BEST
00766       for (iv = 0;  iv < nbest;  iv++)
00767         if (sse < response[iv]){
00768            int jv ;
00769            for( jv=nbest-1 ; jv > iv ; jv-- ){  /* the push-up code */
00770               response[jv] = response[jv-1] ;
00771               for( ip=0 ; ip < r+p ; ip++ )
00772                  parameters[jv][ip] = parameters[jv-1][ip] ;
00773            }
00774 
00775            response[iv] = sse ;                 /* now can copy new */
00776            for( ip=0 ; ip < r+p ;  ip++ )       /* values into place */
00777               parameters[iv][ip] = par[ip] ;
00778            break ;
00779         }
00780 #else
00781       for (iv = 0;  iv < nbest;  iv++)          /* this is the old code */
00782         if (sse < response[iv])
00783           {
00784             for (ip = 0;  ip < r+p;  ip++)
00785               parameters[iv][ip] = par[ip];
00786             response[iv] = sse;
00787             break;
00788           }
00789 #endif
00790   } /* end of loop over random vectors */
00791 
00792   /*----- release memory space -----*/
00793   free (par);       par = NULL;
00794 
00795 #ifdef SAVE_RAN
00796       RAN_sind = -1 ;
00797 #endif
00798 
00799 #if 0
00800    /*** Some debugging stuff -- 28 July 1998 ***/
00801    { static count=-1 ;
00802      count++ ;
00803      if( count%10 == 0 ){
00804         printf("Random search sse:\n") ;
00805         for( iv=0 ; iv < nbest ; iv++ ) printf(" %g",response[iv]) ;
00806         printf("\n") ;
00807      }
00808    }
00809 #endif
00810  
00811 }
00812 
00813 
00814 /*---------------------------------------------------------------------------*/
00815 /*
00816   Estimate the parameters for the full model.
00817 */
00818 
00819 void calc_full_model 
00820 (
00821   vfp nmodel,             /* pointer to noise model */
00822   vfp smodel,             /* pointer to signal model */
00823   int r,                  /* number of parameters in the noise model */
00824   int p,                  /* number of parameters in the signal model */
00825   float * min_nconstr,    /* minimum parameter constraints for noise model */
00826   float * max_nconstr,    /* maximum parameter constraints for noise model */
00827   float * min_sconstr,    /* minimum parameter constraints for signal model */
00828   float * max_sconstr,    /* maximum parameter constraints for signal model */
00829   int ts_length,          /* length of time series array */
00830   float ** x_array,       /* independent variable matrix */
00831   float * ts_array,       /* observed time series */
00832   float * par_rdcd,       /* estimated parameters for the reduced model */
00833   float sse_rdcd,         /* error sum of squares for the reduced model */
00834   int nabs,               /* use absolute constraints for noise parameters */
00835   int nrand,              /* number of random vectors to generate */
00836   int nbest,              /* number of random vectors to keep */
00837   float rms_min,          /* minimum rms error required to reject rdcd model */
00838   float * par_full,       /* estimated parameters for the full model */
00839   float * sse_full,       /* error sum of squares for the full model */
00840   int * novar             /* flag for insufficient variation in data */
00841 )
00842 
00843 {
00844   int iv;                        /* vector index */
00845   int ip;                        /* parameter index */ 
00846   float ** parameters = NULL;    /* array of parameter vectors */
00847   float * sse = NULL;            /* array of sse's */
00848 
00849 
00850   /*----- if this is the null signal model, 
00851           or if rms error for reduced model is very small, 
00852           just use the reduced model -----*/
00853   if ( (p < 1) || (sqrt(sse_rdcd/(ts_length - r)) < rms_min) ) 
00854     {
00855       *novar = 1;
00856       for (ip = 0;  ip < r;  ip++)
00857         par_full[ip] = par_rdcd[ip];
00858       for (ip = r;  ip < r+p;  ip++)
00859         par_full[ip] = 0.0;
00860       *sse_full = sse_rdcd;
00861       return;
00862     }
00863   else
00864     *novar = 0;
00865 
00866 
00867   /*----- initialize random number generator -----*/
00868   srand48 (1234567);
00869 
00870   /*----- allocate memory for calculation of full model-----*/
00871   initialize_full_model (r+p, nbest, &parameters, &sse); 
00872 
00873   /*----- evaluate randomly chosen vectors in the parameter space -----*/
00874   random_search (nmodel, smodel, r, p, nabs,
00875                  min_nconstr, max_nconstr, min_sconstr, max_sconstr,
00876                  ts_length, x_array, ts_array, par_rdcd, nrand, nbest, 
00877                  parameters, sse);
00878   
00879   /*----- use the best random vectors as the starting points for 
00880     nonlinear optimization -----*/
00881   for (iv = 0;  iv < nbest;  iv++)
00882     {
00883       simplex_optimization (nmodel, smodel, r, p, 
00884                             min_nconstr, max_nconstr, min_sconstr, max_sconstr,
00885                             nabs, ts_length, x_array, ts_array, par_rdcd,
00886                             parameters[iv], &sse[iv]);
00887     }
00888 
00889   /*----- save the best result from the nonlinear optimization -----*/      
00890   *sse_full = 1.0e+30;
00891   for (iv = 0;  iv < nbest;  iv++)
00892     {
00893       if (sse[iv] < *sse_full)
00894         {
00895           *sse_full = sse[iv];
00896           for (ip = 0;  ip < r+p;  ip++)
00897             par_full[ip] = parameters[iv][ip];
00898         }
00899     }
00900 
00901 
00902   /*----- release memory space -----*/
00903   for (iv = 0;  iv < nbest;  iv++)
00904     {
00905       free (parameters[iv]);
00906       parameters[iv] = NULL;
00907     }
00908   free (parameters);       parameters = NULL;     
00909   free (sse);              sse = NULL;
00910 
00911 }
00912 
00913 
00914 /*---------------------------------------------------------------------------*/
00915 /*
00916   Routine to calculate the partial derivative matrix, evaluated at the
00917   estimated parameter location.
00918 */
00919 
00920 void calc_partial_derivatives 
00921 (
00922   vfp nmodel,             /* pointer to noise model */
00923   vfp smodel,             /* pointer to signal model */
00924   int r,                  /* number of parameters in the noise model */
00925   int p,                  /* number of parameters in the signal model */
00926   float * min_nconstr,    /* minimum parameter constraints for noise model */
00927   float * max_nconstr,    /* maximum parameter constraints for noise model */
00928   float * min_sconstr,    /* minimum parameter constraints for signal model */
00929   float * max_sconstr,    /* maximum parameter constraints for signal model */
00930   int ts_length,          /* length of time series data */
00931   float ** x_array,       /* independent variable matrix */
00932   float * par_full,       /* estimated parameters for the full model */
00933   matrix d                /* matrix of estimated partial derivatives */
00934 )
00935 
00936 {
00937   const float EPSILON = 1.0e-10;
00938   int dimension;          /* dimension of full model */
00939   int ip, jp;             /* parameter indices */ 
00940   int it;                 /* time index */
00941   float delp;             /* delta parameter */
00942   float * par = NULL;            /* perturbed parameter array */
00943   float * ts1_array = NULL;      /* estimated time series */
00944   float * ts2_array = NULL;      /* perturbed estimated time series */
00945 
00946 
00947   /*----- dimension of full model -----*/
00948   dimension = r + p;
00949 
00950   /*----- allocate memory -----*/
00951   ts1_array = (float *) malloc (sizeof(float) * ts_length);
00952   ts2_array = (float *) malloc (sizeof(float) * ts_length);
00953   par = (float *) malloc (sizeof(float) * dimension);
00954 
00955   /*----- fitted time series at estimated parameter vector -----*/
00956   full_model (nmodel, smodel, par_full, par_full + r, 
00957               ts_length, x_array, ts1_array);
00958 
00959   /*----- loop over parameters in model -----*/
00960   for (jp = 0;  jp < dimension;  jp++)
00961     {
00962       /*----- initialize parameters -----*/
00963       for (ip = 0;  ip < dimension;  ip++)
00964         par[ip] = par_full[ip];
00965 
00966       /*----- add small increment to the jpth parameter -----*/
00967       if (jp < r)
00968         delp = (max_nconstr[jp] - min_nconstr[jp]) / 1000.0;
00969       else
00970         delp = (max_sconstr[jp-r] - min_sconstr[jp-r]) / 1000.0;
00971       par[jp] += delp;
00972 
00973       /*----- fit time series for perturbed parameter -----*/
00974       full_model (nmodel, smodel, par, par + r, 
00975                   ts_length, x_array, ts2_array);
00976 
00977       /*----- estimate partial derivative -----*/
00978       if (delp > EPSILON)
00979         for (it = 0;  it < ts_length;  it++)
00980           d.elts[it][jp] = (ts2_array[it] - ts1_array[it]) / delp;
00981       else
00982         for (it = 0;  it < ts_length;  it++)
00983           d.elts[it][jp] = 0.0;
00984       
00985     }
00986 
00987   /*----- free memory -----*/
00988   free (par);        par = NULL;
00989   free (ts2_array);  ts2_array = NULL;
00990   free (ts1_array);  ts1_array = NULL;
00991 
00992 }  
00993 
00994 
00995 /*---------------------------------------------------------------------------*/
00996 /*
00997   Calculate the coefficient of multiple determination R^2.
00998 */
00999 
01000 float calc_rsqr 
01001 (
01002   float ssef,                 /* error sum of squares from full model */
01003   float sser                  /* error sum of squares from reduced model */
01004 )
01005 
01006 {
01007   const float EPSILON = 1.0e-5;     /* protection against divide by zero */
01008   float rsqr;                       /* coeff. of multiple determination R^2  */
01009 
01010 
01011   /*----- coefficient of multiple determination R^2 -----*/
01012   if (sser < EPSILON)
01013     rsqr = 0.0;
01014   else
01015     rsqr = (sser - ssef) / sser;
01016 
01017 
01018   /*----- Limit range of values for R^2 -----*/
01019   if (rsqr < 0.0)   rsqr = 0.0;
01020   if (rsqr > 1.0)   rsqr = 1.0;
01021 
01022 
01023   /*----- Return coefficient of multiple determination R^2 -----*/
01024   return (rsqr);
01025 }
01026 
01027 
01028 /*---------------------------------------------------------------------------*/
01029 /*
01030   Perform additional calculations following the least squares parameter
01031   estimates.
01032 */
01033 
01034 void analyze_results 
01035 (
01036   vfp nmodel,             /* pointer to noise model */
01037   vfp smodel,             /* pointer to signal model */
01038   int r,                  /* number of parameters in the noise model */
01039   int p,                  /* number of parameters in the signal model */
01040   int novar,              /* flag for insufficient variation in the data */
01041   float * min_nconstr,    /* minimum parameter constraints for noise model */
01042   float * max_nconstr,    /* maximum parameter constraints for noise model */
01043   float * min_sconstr,    /* minimum parameter constraints for signal model */
01044   float * max_sconstr,    /* maximum parameter constraints for signal model */
01045   int ts_length,          /* length of time series data */
01046   float ** x_array,       /* independent variable matrix */
01047   float * par_rdcd,       /* estimated parameters for the reduced model */
01048   float sse_rdcd,         /* error sum of squares for the reduced model */ 
01049   float * par_full,       /* estimated parameters for the full model */
01050   float sse_full,         /* error sum of squares for the full model */
01051   float * rmsreg,         /* root-mean-square for the full regression model */
01052   float * freg,           /* f-statistic for the full regression model */
01053   float * rsqr,           /* R^2 (coef. of multiple determination) */
01054   float * smax,           /* signed maximum of signal */
01055   float * tmax,           /* epoch of signed maximum of signal */
01056   float * pmax,           /* percentage change due to signal */
01057   float * area,           /* area between signal and baseline */
01058   float * parea,          /* percentage area between signal and baseline */
01059   float * tpar_full       /* t-statistic of the parameters in the full model */
01060 )
01061 
01062 {
01063   const float EPSILON = 1.0e-10;
01064   int dimension;                 /* dimension of full model */
01065   int ip;                        /* parameter index */
01066   int df_rdcd;                   /* degrees of freedom for reduced model */
01067   int df_full;                   /* degrees of freedom for full model */
01068   float mse_full;                /* MSE for full model */
01069   float mse_reg;                 /* MSE for the regression */
01070   int ok;                        /* boolean for successful matrix inversion */
01071   float * y_array = NULL;        /* estimated signal time series */
01072   float * base_array = NULL;     /* baseline time series */
01073   float barea;                   /* area under baseline */
01074   int it;                        /* time series index */
01075   float y1, y2, y3;              /* temporary values for calculating area */
01076 
01077 
01078   /*----- initialization -----*/
01079   dimension = r + p;
01080   *rmsreg = *freg = *rsqr = *smax = *tmax = *pmax = *area = *parea = 0.0;
01081   for (ip = 0;  ip < dimension;  ip++)
01082     tpar_full[ip] = 0.0;
01083 
01084 
01085   /*----- check for insufficient variation in the data -----*/
01086   if (novar)  return;
01087 
01088 
01089   /*----- Adjust dof if constraints force a parameter to be a constant -----*/
01090   df_rdcd = ts_length - r;
01091   df_full = ts_length - dimension;
01092 
01093   for (ip = 0;  ip < r;  ip++)
01094     if (min_nconstr[ip] == max_nconstr[ip])
01095       {
01096         df_rdcd++;
01097         df_full++;
01098       }
01099 
01100   for (ip = 0;  ip < p;  ip++)
01101     if (min_sconstr[ip] == max_sconstr[ip])
01102       df_full++;
01103 
01104 
01105   /*----- MSE for full model -----*/
01106   mse_full = sse_full / df_full;
01107 
01108   /*----- MSE due to regression -----*/
01109   if (df_rdcd == df_full)
01110     mse_reg = 0.0;
01111   else
01112     mse_reg = (sse_rdcd - sse_full) / (df_rdcd - df_full);
01113   if (mse_reg < 0.0) mse_reg = 0.0;
01114 
01115   /*----- f-statistic for significance of the regression -----*/
01116   if (mse_full > EPSILON)
01117     *freg = mse_reg / mse_full;
01118   else
01119     *freg = 0.0;
01120 
01121   /*----- root-mean-square error for the regression -----*/
01122   *rmsreg = sqrt(mse_full);
01123 
01124   /*----- R^2 (coefficient of multiple determination) -----*/
01125   *rsqr = calc_rsqr (sse_full, sse_rdcd);
01126 
01127   /*----- generate time series corresponding to the signal model -----*/
01128   y_array = (float *) malloc (sizeof(float) * (ts_length));
01129   if (y_array == NULL)
01130     NLfit_error ("Unable to allocate memory for y_array");
01131 #if 0
01132   smodel (par_full+r, ts_length, x_array, y_array);
01133 #else
01134   AFNI_CALL_VOID_4ARG(smodel ,
01135                       float *,(par_full+r) , int,ts_length,
01136                       float **,x_array , float *,y_array ) ;
01137 #endif
01138 
01139   /*----- generate time series corresponding to the noise model -----*/
01140   base_array = (float *) malloc (sizeof(float) * (ts_length));
01141   if (base_array == NULL)
01142     NLfit_error ("Unable to allocate memory for base_array");
01143 #if 0
01144   nmodel (par_full, ts_length, x_array, base_array);
01145 #else
01146   AFNI_CALL_VOID_4ARG(nmodel ,
01147                       float *,par_full , int,ts_length,
01148                       float **,x_array , float *,base_array ) ;
01149 #endif
01150 
01151   /*----- initialize signal parameters -----*/
01152   *tmax = x_array[0][1];
01153   *smax = y_array[0];
01154   if (fabs(base_array[0]) > EPSILON)
01155     *pmax = 100.0 * y_array[0] / fabs(base_array[0]);
01156   else
01157     *pmax = 0.0;
01158   *area = 0.0;
01159   barea = 0.0;
01160   *parea = 0.0;
01161 
01162   /*----- calculate signed maximum of the signal, percent change, area -----*/
01163   for (it = 1;  it < ts_length;  it++)
01164     {
01165       /*----- calculate signed maximum of signal, and
01166         calculate max. percentage change of signal wrt baseline -----*/
01167       if (fabs(y_array[it]) > fabs(*smax))
01168         {
01169           *tmax = x_array[it][1];
01170           *smax = y_array[it];
01171           if (fabs(base_array[it]) > EPSILON)
01172             *pmax = 100.0 * y_array[it] / fabs(base_array[it]);
01173         }
01174       
01175       /*----- calculate area and percentage area under the signal -----*/
01176       y1 = y_array[it-1];   y2 = y_array[it];
01177       if ((y1 > 0.0) && (y2 > 0.0))
01178         {
01179           *area += (y1 + y2) / 2.0;
01180           *parea += (y1 + y2) / 2.0;
01181         }
01182       else if ((y1 < 0.0) && (y2 < 0.0))
01183         {
01184           *area -= (y1 + y2) / 2.0;
01185           *parea += (y1 + y2) / 2.0;
01186         }
01187       else
01188         {
01189           y3 = fabs(y1) + fabs(y2);
01190           if (y3 > EPSILON)
01191             {
01192               *area += (y1*y1 + y2*y2) / (2.0 * y3);
01193               if (y1 > y2)
01194                 *parea += (y1*y1 - y2*y2) / (2.0 * y3);
01195               else
01196                 *parea -= (y1*y1 - y2*y2) / (2.0 * y3);
01197             }
01198         }
01199       
01200       y1 = base_array[it-1];   y2 = base_array[it];
01201       if ((y1 > 0.0) && (y2 > 0.0))
01202         {
01203           barea += (y1 + y2) / 2.0;
01204         }
01205       else if ((y1 < 0.0) && (y2 < 0.0))
01206         {
01207           barea -= (y1 + y2) / 2.0;
01208         }
01209       else
01210         {
01211           y3 = fabs(y1) + fabs(y2);
01212           if (y3 > EPSILON)
01213             {
01214               barea += (y1*y1 + y2*y2) / (2.0 * y3);
01215             }
01216         }
01217     }  /* it */
01218 
01219   if (barea > EPSILON)
01220     *parea *= 100.0/barea;
01221   else
01222     *parea = 0.0;
01223 
01224   free (base_array);   base_array = NULL;
01225   free (y_array);      y_array = NULL;
01226 
01227 
01228   /*----- Calculate t-statistics? -----*/
01229   if (calc_tstats)
01230     {
01231       float stddev;                 /* est. std. dev. for a parameter */
01232       matrix d, dt, dtd, dtdinv;    /* matrices used for calc. of covariance */
01233 
01234       /*----- initialize matrices -----*/
01235       matrix_initialize (&d);
01236       matrix_initialize (&dt);
01237       matrix_initialize (&dtd);
01238       matrix_initialize (&dtdinv);
01239       
01240       /*----- calculate the matrix of partial derivatives D -----*/
01241       matrix_create (ts_length, dimension, &d);
01242       calc_partial_derivatives (nmodel, smodel, r, p, 
01243                                 min_nconstr, max_nconstr, 
01244                                 min_sconstr, max_sconstr, 
01245                                 ts_length, x_array, par_full, d);
01246       
01247       /*----- calculate variance-covariance matrix -----*/
01248       matrix_transpose (d, &dt);
01249       matrix_multiply (dt, d, &dtd);
01250       ok = matrix_inverse (dtd, &dtdinv);
01251       if (ok)
01252         for (ip = 0;  ip < dimension;  ip++)
01253           {
01254             stddev = sqrt((sse_full/(df_full)) * dtdinv.elts[ip][ip]);
01255             if (stddev > EPSILON)
01256               tpar_full[ip] = par_full[ip] / stddev;
01257             else
01258               tpar_full[ip] = 0.0;
01259           }
01260       else
01261         for (ip = 0;  ip < dimension;  ip++)
01262           {
01263             tpar_full[ip] = 0.0;
01264           }
01265       
01266       
01267       /*----- dispose of matrices -----*/
01268       matrix_destroy (&dtdinv);
01269       matrix_destroy (&dtd);
01270       matrix_destroy (&dt); 
01271       matrix_destroy (&d);
01272 
01273     }  /* if (calc_tstats) */    
01274 
01275 }
01276 
01277 
01278 /*---------------------------------------------------------------------------*/
01279 /*
01280   Convert F-value to p-value.  
01281   This routine was copied from: mri_stats.c
01282 */
01283 
01284 
01285 double fstat_t2p( double ff , double dofnum , double dofden )
01286 {
01287    int which , status ;
01288    double p , q , f , dfn , dfd , bound ;
01289 
01290    which  = 1 ;
01291    p      = 0.0 ;
01292    q      = 0.0 ;
01293    f      = ff ;
01294    dfn    = dofnum ;
01295    dfd    = dofden ;
01296 
01297    cdff( &which , &p , &q , &f , &dfn , &dfd , &status , &bound ) ;
01298 
01299    if( status == 0 ) return q ;
01300    else              return 1.0 ;
01301 }
01302 
01303 /*---------------------------------------------------------------------------*/
01304 /*
01305   Report results for this voxel.
01306 */
01307 
01308 void report_results 
01309 (
01310   char * nname,            /* name of noise model */
01311   char * sname,            /* name of signal model */
01312   int r,                   /* number of parameters in the noise model */
01313   int p,                   /* number of parameters in the signal model */
01314   char ** npname,          /* noise parameter names */
01315   char ** spname,          /* signal parameter names */
01316   int ts_length,           /* length of time series array */
01317   float * par_rdcd,        /* estimated parameters for the reduced model */
01318   float sse_rdcd,          /* error sum of squares for the reduced model */
01319   float * par_full,        /* estimated parameters for the full model */
01320   float sse_full,          /* error sum of squares for the full model */
01321   float * tpar_full,       /* t-statistic of the parameters in full model */
01322   float rmsreg,            /* root-mean-square error for the full model */
01323   float freg,              /* f-statistic for the full regression model */
01324   float rsqr,              /* R^2 (coef. of multiple determination) */
01325   float smax,              /* signed maximum of signal */
01326   float tmax,              /* epoch of signed maximum of signal */
01327   float pmax,              /* percentage change due to signal */
01328   float area,              /* area between signal and baseline */
01329   float parea,             /* percentage area between signal and baseline */
01330   char ** label            /* label output for this voxel */
01331 )
01332 
01333 {
01334   int ip;                  /* parameter index */
01335   double pvalue;
01336 
01337 
01338   /*----- create label if desired by calling program -----*/
01339   if (label == NULL)  return;
01340 
01341   /*----- make this a 0 length string to start -----*/
01342   lbuf[0] = '\0';
01343       
01344   /*----- write reduced model parameters -----*/
01345   sprintf (sbuf, "Reduced (%s) Model: \n", nname);
01346   strcat (lbuf, sbuf);
01347   for (ip = 0;  ip < r;  ip++)
01348     {
01349       sprintf (sbuf, "b[%d]= %12.6f  %s \n", ip, par_rdcd[ip], npname[ip]); 
01350       strcat (lbuf, sbuf);
01351     }
01352     
01353   /*----- write full model parameters -----*/  
01354   sprintf (sbuf, "\nFull (%s + %s) Model: \n", nname, sname);
01355   strcat (lbuf, sbuf);
01356   for (ip = 0;  ip < r;  ip++)
01357     {
01358       sprintf (sbuf, "gn[%d]=%12.6f  %s \n", ip, par_full[ip], npname[ip]);
01359       strcat (lbuf, sbuf);
01360                     
01361     }
01362   for (ip = 0;  ip < p;  ip++)
01363     {
01364       sprintf (sbuf, "gs[%d]=%12.6f  %s \n", ip, par_full[ip+r], spname[ip]);
01365       strcat (lbuf, sbuf);            
01366     }
01367   
01368   sprintf (sbuf, "\nSignal Tmax  = %12.3f \n", tmax);
01369   strcat (lbuf, sbuf);
01370   sprintf (sbuf,   "Signal Smax  = %12.3f \n", smax);
01371   strcat (lbuf, sbuf);
01372   sprintf (sbuf,   "Signal PSmax = %12.3f \n", pmax);
01373   strcat (lbuf, sbuf);
01374   sprintf (sbuf,   "Signal Area  = %12.3f \n", area);
01375   strcat (lbuf, sbuf);
01376   sprintf (sbuf,   "Signal PArea = %12.3f \n", parea);
01377   strcat (lbuf, sbuf);
01378 
01379   sprintf (sbuf, "\nRMSE Rdcd = %8.3f \n", sqrt(sse_rdcd/(ts_length-r)));
01380   strcat (lbuf, sbuf);
01381   sprintf (sbuf, "RMSE Full = %8.3f \n", sqrt(sse_full/(ts_length-r-p)));
01382   strcat (lbuf, sbuf);
01383             
01384   sprintf (sbuf, "\nR^2       = %7.3f \n", rsqr);
01385   strcat (lbuf, sbuf);
01386   sprintf (sbuf, "F[%2d,%3d] = %7.3f \n", p, ts_length-r-p, freg);
01387   strcat (lbuf, sbuf);
01388   pvalue = fstat_t2p ( (double) freg, (double) p, (double) ts_length-r-p);
01389   sprintf (sbuf, "p-value   = %e  \n", pvalue);
01390   strcat (lbuf, sbuf);
01391  
01392 
01393   /*----- send address of lbuf back in what label points to -----*/
01394   *label = lbuf;
01395   
01396 }
01397 
01398 
01399 
01400 
01401 
01402 
01403 
01404 
01405 
01406 
01407 
01408 
01409