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 static int     refnum = 0 ;     /* # pts in refts */
00010 static int     refnz  = 0 ;     /* # of nonzero pts */
00011 static float * refts  = NULL ;  /* reference time series */
00012 static int   * refin  = NULL ;  /* indexes of nonzero pts */
00013 
00014 void gamma_model( float * , int , float ** , float * ) ;
00015 
00016 #define ERREX(str) ( fprintf(stderr,"\n*** %s\a\n",str) , exit(1) )
00017 
00018 /*----------------------------------------------------------------------
00019    Function to set the reference time series, with which the
00020    model function is convolved to produce the simulated data.
00021 ------------------------------------------------------------------------*/
00022 
00023 void conv_set_ref( int num , float * ref )
00024 {
00025    if( num > 0 && ref != NULL ){ /*** if have inputs, make space & copy in ***/
00026       int ii ;
00027 
00028       /* get rid of old data */
00029 
00030       if( refts != NULL ){ free(refts); refts = NULL; free(refin); refin = NULL; }
00031 
00032       refnum = num ;
00033       refts  = (float *) malloc( sizeof(float) * num ) ;
00034       refin  = (int *)   malloc( sizeof(int)   * num ) ;
00035       memcpy( refts , ref , sizeof(float) * num ) ;
00036       for( ii=0,refnz=0 ; ii < num ; ii++ )              /* build list of nonzero */
00037          if( refts[ii] != 0 ) refin[refnz++] = ii ;      /* points in refts */
00038       if( refnz == 0 )
00039          ERREX("model_convgamma: All zero reference timeseries!") ;
00040 
00041 #if 0
00042 fprintf(stderr,"conv_set_ref: num=%d nonzero=%d\n",num,refnz) ;
00043 #endif
00044 
00045       return ;
00046 
00047    } else { /*** if no inputs, do something special ***/
00048 
00049      char * cp ;
00050      MRI_IMAGE * flim ;
00051      float one = 1.0 ;
00052 
00053      cp = my_getenv("AFNI_CONVMODEL_REF") ;  /* get name of reference file */
00054      if( cp == NULL )
00055         ERREX("model_convgamma: Can't read AFNI_CONVMODEL_REF from environment") ;
00056 
00057      flim = mri_read_1D(cp) ;            /* 16 Nov 1999: replaces mri_read_ascii */
00058      if( flim == NULL ){
00059         char buf[256] ;
00060         sprintf(buf,"model_convgamma: Can't read timeseries file %s",cp) ;
00061         ERREX(buf) ;
00062      }
00063 
00064 #if 0
00065 fprintf(stderr,"conv_set_ref: refts=%s  nx=%d\n",cp,flim->ny) ;
00066 #endif
00067 
00068      conv_set_ref( flim->nx , MRI_FLOAT_PTR(flim) ) ;  /* recursion! */
00069      mri_free(flim) ;
00070    }
00071    return ;
00072 }
00073 
00074 /*-----------------------------------------------------------------------
00075   Function to compute the simulated time series.
00076 -------------------------------------------------------------------------*/
00077 
00078 void conv_model( float *  gs      , int     ts_length ,
00079                  float ** x_array , float * ts_array   )
00080 {
00081    int ii, jj,jbot,jtop , kk , nid_top,nid_bot ;
00082    float top , val ;
00083 
00084    static int     nid = 0 ;     /* number of pts in impulse */
00085    static float * fid = NULL ;  /* impulse response function */
00086 
00087    /*** make sure there is a reference function to convolve with ***/
00088 
00089    if( refnum <= 0 ) conv_set_ref( 0 , NULL ) ;
00090 
00091    /*** initialize the output ***/
00092 
00093    for( ii=0 ; ii < ts_length ; ii++ ) ts_array[ii] = 0.0 ;
00094 
00095    /*** initialize the impulse response ***/
00096 
00097    if( nid < ts_length ){              /* make some space for it */
00098       if( fid != NULL ) free(fid) ;
00099       nid = ts_length ;
00100       fid = (float *) malloc( sizeof(float) * nid ) ;
00101    }
00102 
00103    gamma_model( gs , ts_length , x_array , fid ) ;  /* compute impulse */
00104 
00105    top = 0.0 ;                                      /* find max value */
00106    for( jj=0 ; jj < ts_length ; jj++ ){
00107       val = fabs(fid[jj]) ; if( val > top ) top = val ;
00108    }
00109    if( top == 0.0 ) fid[0] = 1.0 ;                  /* very unlikely case */
00110    top *= 0.001 ;
00111    for( jj=0 ; jj < ts_length ; jj++ ){             /* discard small values */
00112       if( fabs(fid[jj]) < top ) fid[jj] = 0.0 ;
00113    }
00114    for( nid_bot=0 ; nid_bot < ts_length ; nid_bot++ )         /* find first nonzero */
00115       if( fid[nid_bot] != 0.0 ) break ;
00116    for( nid_top=ts_length-1 ; nid_top > nid_bot ; nid_top-- ) /* and last nonzero */
00117       if( fid[nid_top] != 0.0 ) break ;
00118 
00119    /*** loop over each nonzero point in the reference ***/
00120 
00121    for( ii=0 ; ii < refnz ; ii++ ){
00122       kk  = refin[ii] ; if( kk >= ts_length ) break ;
00123       val = refts[kk] ;
00124 
00125       /*** for each point in the impulse ***/
00126 
00127       jtop = ts_length - kk ; if( jtop > nid_top ) jtop = nid_top+1 ;
00128       for( jj=nid_bot ; jj < jtop ; jj++ )
00129          ts_array[kk+jj] += val * fid[jj] ;
00130    }
00131 
00132    return ;
00133 }
00134 
00135 /*-----------------------------------------------------------------------*/
00136 
00137 DEFINE_MODEL_PROTOTYPE
00138 
00139 MODEL_interface * initialize_model ()
00140 {
00141   MODEL_interface * mi = NULL;
00142 
00143   /*----- allocate memory space for model interface -----*/
00144 
00145   mi = (MODEL_interface *) XtMalloc (sizeof(MODEL_interface));
00146 
00147   /*----- name of this model -----*/
00148 
00149   strcpy (mi->label, "ConvGamma");
00150 
00151   /*----- this is a signal model -----*/
00152 
00153   mi->model_type = MODEL_SIGNAL_TYPE;
00154 
00155   /*----- number of parameters in the model -----*/
00156 
00157   mi->params = 4;
00158 
00159   /*----- parameter labels -----*/
00160 
00161   strcpy (mi->plabel[0], "t0");
00162   strcpy (mi->plabel[1], "amp");
00163   strcpy (mi->plabel[2], "r");
00164   strcpy (mi->plabel[3], "b");
00165 
00166   /*----- minimum and maximum parameter constraints -----*/
00167 
00168   mi->min_constr[0] =     0.0;    mi->max_constr[0] =    10.0;
00169   mi->min_constr[1] =     0.0;    mi->max_constr[1] =   200.0;
00170   mi->min_constr[2] =     1.0;    mi->max_constr[2] =    19.0;
00171   mi->min_constr[3] =     0.1;    mi->max_constr[3] =     5.0;
00172 
00173   /*----- function which implements the model -----*/
00174   mi->call_func = &conv_model;
00175 
00176   return (mi);
00177 }
00178 
00179 /*----------------------------------------------------------------------*/
00180 /*
00181   Routine to calculate the time series which results from using the
00182   gamma variate drug response signal model with the specified
00183   model parameters.
00184 
00185   Definition of model parameters:
00186 
00187          gs[0] = time delay of response (t0)
00188          gs[1] = multiplicative constant (k)
00189          gs[2] = rise rate exponent (r)
00190          gs[3] = decay rate constant (b)
00191 
00192 */
00193 
00194 void gamma_model
00195 (
00196   float * gs,                /* parameters for signal model */
00197   int ts_length,             /* length of time series data */
00198   float ** x_array,          /* independent variable matrix */
00199   float * ts_array           /* estimated signal model time series */
00200 )
00201 
00202 {
00203   int it;                           /* time index */
00204   float t;                          /* time */
00205   float gsi , fac ;
00206 
00207   if( gs[3] <= 0.0 || gs[2] <= 0.0 ){
00208      for( it=0 ; it < ts_length ; it++ ) ts_array[it] = 0.0 ;
00209      return ;
00210   }
00211 
00212   /* fac is chosen to make the peak value equal to gs[1] */
00213 
00214   gsi = 1.0 / gs[3] ;
00215   fac = gs[1] * exp( gs[2] * ( 1.0 - log(gs[2]*gs[3]) ) ) ;
00216   for( it=0;  it < ts_length;  it++){
00217      t = x_array[it][1] - gs[0] ;
00218      ts_array[it] = (t <= 0.0) ? 0.0
00219                                : fac * exp( log(t) * gs[2] - t * gsi ) ;
00220   }
00221   return ;
00222 }