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 #include "NLfit_model.h"
00008 
00009 /****************************************************************************
00010      Model function for 3dNLfim:
00011      gamma variate impulse response, convolved with NREF input time series,
00012      each with an independent scaling amplitude, but all other parameters
00013      of the gamma variate fixed.  There are NREF+3 parameters:
00014        t0    = time delay
00015        r     = power
00016        b     = time scale
00017        amp_j = amplitude of j-th reference, for j=1..NREF
00018 
00019      28 July 1998 -- RWCox
00020 ****************************************************************************/
00021 
00022 #define NREF 2
00023 
00024 static int     refnum[NREF] ;  /* # pts in refts */
00025 static int     refnz[NREF]  ;  /* # of nonzero pts */
00026 static float * refts[NREF]  ;  /* reference time series */
00027 static int   * refin[NREF]  ;  /* indexes of nonzero pts */
00028 
00029 void gamma_model( float * , int , float ** , float * ) ;
00030 
00031 #define ERREX(str) ( fprintf(stderr,"\n*** %s\a\n",str) , exit(1) )
00032 
00033 /*----------------------------------------------------------------------
00034    Function to set the reference time series, with which the
00035    model function is convolved to produce the simulated data.
00036      num         = length of input time series
00037      ref[jv][ii] = ii-th point in jv-th time series,
00038                    for ii=0..num-1, jv=0..NREF-1
00039    If num==0, will read the reference time series from a file.
00040 ------------------------------------------------------------------------*/
00041 
00042 void conv_set_ref( int num , float ** ref )
00043 {
00044    int jv ;
00045 
00046    if( num > 0 && ref != NULL ){ /*** if have inputs, make space & copy in ***/
00047       int ii ;
00048 
00049       for( jv=0 ; jv < NREF ; jv++ ){
00050 
00051          /* get rid of old data? */
00052 
00053          if( refts[jv] != NULL ){
00054             free(refts[jv]); refts[jv] = NULL;
00055             free(refin[jv]); refin[jv] = NULL;
00056          }
00057 
00058          /* copy new data */
00059 
00060          refnum[jv] = num ;
00061          refts[jv]  = (float *) malloc( sizeof(float) * num ) ;
00062          refin[jv]  = (int *)   malloc( sizeof(int)   * num ) ;
00063          memcpy( refts[jv] , ref[jv] , sizeof(float) * num ) ;
00064 
00065          /* build a list of nonzero entries in this column */
00066 
00067          for( ii=0,refnz[jv]=0 ; ii < num ; ii++ )
00068             if( refts[jv][ii] != 0.0 ) refin[jv][refnz[jv]++] = ii ;
00069 
00070          if( refnz[jv] == 0 )
00071             ERREX(__FILE__ ": All zero reference timeseries column!") ;
00072       }
00073       return ;
00074 
00075    } else { /*** if no inputs, do something special ***/
00076 
00077      char * cp ;
00078      MRI_IMAGE * flim ;
00079      int jv , nx ;
00080      float * ref[NREF] ;
00081 
00082      cp = my_getenv("AFNI_CONVMODEL_REF") ;  /* get name of reference file */
00083      if( cp == NULL )
00084         ERREX(__FILE__ ": Can't read AFNI_CONVMODEL_REF from environment") ;
00085 
00086      flim = mri_read_1D(cp) ;            /* 16 Nov 1999: replaces mri_read_ascii */
00087      if( flim == NULL ){
00088         char buf[256] ;
00089         sprintf(buf,__FILE__ ": Can't read timeseries file %s",cp) ;
00090         ERREX(buf) ;
00091      } else {
00092         fprintf(stderr,__FILE__ ": Read reference file %s\n",cp) ;
00093      }
00094 
00095      if( flim->ny < NREF )
00096         ERREX(__FILE__ ": reference file has too few columns!") ;
00097      else if( flim->nv > NREF )
00098         fprintf(stderr,__FILE__ " WARNING: reference file has too many columns!\n") ;
00099 
00100      nx = flim->nx ;
00101      for( jv=0 ; jv < NREF ; jv++ )
00102         ref[jv] = MRI_FLOAT_PTR(flim) + jv*nx ;
00103 
00104      conv_set_ref( nx , ref ) ;  /* recursion! */
00105      mri_free(flim) ;
00106    }
00107    return ;
00108 }
00109 
00110 /*-----------------------------------------------------------------------
00111   Function to compute the simulated time series:
00112     gs[0..2]      = gamma variate paramters t0, r, b
00113     gs[3..NREF+3] = amplitudes for each reference time series
00114 -------------------------------------------------------------------------*/
00115 
00116 void conv_model( float *  gs      , int     ts_length ,
00117                  float ** x_array , float * ts_array   )
00118 {
00119    int ii, jj,jbot,jtop , kk , nid_top,nid_bot , jv ;
00120    float top , val , amp ;
00121 
00122    static int     nid = 0 ;     /* number of pts in fid array */
00123    static float * fid = NULL ;  /* impulse response function */
00124 
00125    /*** make sure there is a reference function to convolve with ***/
00126 
00127    if( refnum[0] <= 0 ) conv_set_ref( 0 , NULL ) ;
00128 
00129    /*** initialize the output ***/
00130 
00131    for( ii=0 ; ii < ts_length ; ii++ ) ts_array[ii] = 0.0 ;
00132 
00133    /*** initialize the impulse response ***/
00134 
00135    if( nid < ts_length ){              /* make some space for it */
00136       if( fid != NULL ) free(fid) ;
00137       nid = ts_length ;
00138       fid = (float *) malloc( sizeof(float) * nid ) ;
00139    }
00140 
00141    gamma_model( gs , ts_length , x_array , fid ) ;  /* compute impulse */
00142 
00143    /*** discard small values from the fid ***/
00144 
00145 #define EPS 0.001  /* definition of small */
00146 
00147 #undef FIND_TOP  /* don't need this, since our fid amplitude is set to 1.0 */
00148 #ifdef FIND_TOP
00149    top = 0.0 ;                                      /* find max value */
00150    for( jj=0 ; jj < ts_length ; jj++ ){
00151       val = fabs(fid[jj]) ; if( val > top ) top = val ;
00152    }
00153    if( top == 0.0 ) fid[0] = 1.0 ;                  /* very unlikely case */
00154    top *= EPS ;
00155 #else
00156    top = EPS ;
00157 #endif
00158    for( jj=0 ; jj < ts_length ; jj++ ){             /* discard small values */
00159       if( fabs(fid[jj]) < top ) fid[jj] = 0.0 ;
00160    }
00161    for( nid_bot=0 ; nid_bot < ts_length ; nid_bot++ )         /* find first nonzero */
00162       if( fid[nid_bot] != 0.0 ) break ;
00163    for( nid_top=ts_length-1 ; nid_top > nid_bot ; nid_top-- ) /* and last nonzero */
00164       if( fid[nid_top] != 0.0 ) break ;
00165 
00166    /*** loop over each reference ***/
00167 
00168    for( jv=0 ; jv < NREF ; jv++ ){
00169 
00170       amp = gs[jv+3] ; if( amp == 0.0 ) continue ;
00171 
00172       /*** loop over each nonzero point in the reference ***/
00173 
00174       for( ii=0 ; ii < refnz[jv] ; ii++ ){
00175          kk  = refin[jv][ii] ; if( kk >= ts_length ) break ;
00176          val = amp * refts[jv][kk] ;
00177 
00178          /*** for each point in the impulse ***/
00179 
00180          jtop = ts_length - kk ; if( jtop > nid_top ) jtop = nid_top+1 ;
00181          for( jj=nid_bot ; jj < jtop ; jj++ )
00182             ts_array[kk+jj] += val * fid[jj] ;
00183       }
00184    }
00185 
00186    return ;
00187 }
00188 
00189 /*-----------------------------------------------------------------------*/
00190 
00191 DEFINE_MODEL_PROTOTYPE
00192 
00193 MODEL_interface * initialize_model ()
00194 {
00195   MODEL_interface * mi = NULL;
00196   int jv ;
00197   char buf[32] ;
00198 
00199   /*----- allocate memory space for model interface -----*/
00200 
00201   mi = (MODEL_interface *) XtMalloc (sizeof(MODEL_interface));
00202 
00203   /*----- name of this model -----*/
00204 
00205   sprintf( mi->label , "ConvGamma%da" , NREF ) ;
00206 
00207   /*----- this is a signal model -----*/
00208 
00209   mi->model_type = MODEL_SIGNAL_TYPE;
00210 
00211   /*----- number of parameters in the model -----*/
00212 
00213   mi->params = 3 + NREF ;
00214 
00215   /*----- parameter labels -----*/
00216 
00217   strcpy (mi->plabel[0], "t0");  /* impulse response is proportional to */
00218   strcpy (mi->plabel[1], "r");   /*       r  -(t-t0)/b  */
00219   strcpy (mi->plabel[2], "b");   /* (t-t0)  e           */
00220 
00221   for( jv=0 ; jv < NREF ; jv++ )
00222      sprintf( mi->plabel[jv+3] , "amp_%d" , jv+1 ) ;
00223 
00224   /*----- minimum and maximum parameter constraints -----*/
00225 
00226   mi->min_constr[0] = 0.0;  mi->max_constr[0] = 10.0;  /* delay */
00227   mi->min_constr[1] = 1.0;  mi->max_constr[1] = 19.0;  /* power */
00228   mi->min_constr[2] = 0.1;  mi->max_constr[2] =  5.0;  /* time scale */
00229 
00230   for( jv=0 ; jv < NREF ; jv++ ){
00231      mi->min_constr[jv+3] = 0.0;  mi->max_constr[jv+3] = 200.0;
00232   }
00233 
00234   /*----- function which implements the model -----*/
00235   mi->call_func = &conv_model;
00236 
00237   return (mi);
00238 }
00239 
00240 /*----------------------------------------------------------------------*/
00241 /*
00242   Routine to calculate the time series which results from using the
00243   gamma variate drug response signal model with the specified
00244   model parameters.
00245 
00246   Definition of model parameters:
00247 
00248          gs[0] = time delay of response (t0)
00249          gs[1] = rise rate exponent (r)
00250          gs[2] = decay rate constant (b)
00251 
00252   Time to peak is r * b ;
00253   FWHM of the peak is about 2.3 * sqrt(r) * b, for r > 1.
00254 
00255 */
00256 
00257 void gamma_model
00258 (
00259   float * gs,                /* parameters for signal model */
00260   int ts_length,             /* length of time series data */
00261   float ** x_array,          /* independent variable matrix */
00262   float * ts_array           /* estimated signal model time series */
00263 )
00264 
00265 {
00266   int it;                           /* time index */
00267   float t;                          /* time */
00268   float gsi , fac ;
00269 
00270   if( gs[2] <= 0.0 || gs[1] <= 0.0 ){
00271      ts_array[0] = 1.0 ;
00272      for( it=1 ; it < ts_length ; it++ ) ts_array[it] = 0.0 ;
00273      return ;
00274   }
00275 
00276   /* fac is chosen to make the peak value equal to 1 (at t = gs[1]*gs[2]) */
00277 
00278   gsi = 1.0 / gs[2] ;
00279   fac = exp( gs[1] * ( 1.0 - log(gs[1]*gs[2]) ) ) ;
00280   for( it=0;  it < ts_length;  it++){
00281      t = x_array[it][1] - gs[0] ;
00282      ts_array[it] = (t <= 0.0) ? 0.0
00283                                : fac * exp( log(t) * gs[1] - t * gsi ) ;
00284   }
00285   return ;
00286 }