/***************************************************************************** Major portions of this software are copyrighted by the Medical College of Wisconsin, 1994-2000, and are released under the Gnu General Public License, Version 2. See the file README.Copyright for details. ******************************************************************************/ #include "NLfit_model.h" static int refnum = 0 ; /* # pts in refts */ static int refnz = 0 ; /* # of nonzero pts */ static float * refts = NULL ; /* reference time series */ static int * refin = NULL ; /* indexes of nonzero pts */ void gamma_model( float * , int , float ** , float * ) ; #define ERREX(str) ( fprintf(stderr,"\n*** %s\a\n",str) , exit(1) ) /*---------------------------------------------------------------------- Function to set the reference time series, with which the model function is convolved to produce the simulated data. ------------------------------------------------------------------------*/ void conv_set_ref( int num , float * ref ) { if( num > 0 && ref != NULL ){ /*** if have inputs, make space & copy in ***/ int ii ; /* get rid of old data */ if( refts != NULL ){ free(refts); refts = NULL; free(refin); refin = NULL; } refnum = num ; refts = (float *) malloc( sizeof(float) * num ) ; refin = (int *) malloc( sizeof(int) * num ) ; memcpy( refts , ref , sizeof(float) * num ) ; for( ii=0,refnz=0 ; ii < num ; ii++ ) /* build list of nonzero */ if( refts[ii] != 0 ) refin[refnz++] = ii ; /* points in refts */ if( refnz == 0 ) ERREX("model_convgamma: All zero reference timeseries!") ; #if 0 fprintf(stderr,"conv_set_ref: num=%d nonzero=%d\n",num,refnz) ; #endif return ; } else { /*** if no inputs, do something special ***/ char * cp ; MRI_IMAGE * flim ; float one = 1.0 ; cp = my_getenv("AFNI_CONVMODEL_REF") ; /* get name of reference file */ if( cp == NULL ) ERREX("model_convgamma: Can't read AFNI_CONVMODEL_REF from environment") ; flim = mri_read_1D(cp) ; /* 16 Nov 1999: replaces mri_read_ascii */ if( flim == NULL ){ char buf[256] ; sprintf(buf,"model_convgamma: Can't read timeseries file %s",cp) ; ERREX(buf) ; } #if 0 fprintf(stderr,"conv_set_ref: refts=%s nx=%d\n",cp,flim->ny) ; #endif conv_set_ref( flim->nx , MRI_FLOAT_PTR(flim) ) ; /* recursion! */ mri_free(flim) ; } return ; } /*----------------------------------------------------------------------- Function to compute the simulated time series. -------------------------------------------------------------------------*/ void conv_model( float * gs , int ts_length , float ** x_array , float * ts_array ) { int ii, jj,jbot,jtop , kk , nid_top,nid_bot ; float top , val ; static int nid = 0 ; /* number of pts in impulse */ static float * fid = NULL ; /* impulse response function */ /*** make sure there is a reference function to convolve with ***/ if( refnum <= 0 ) conv_set_ref( 0 , NULL ) ; /*** initialize the output ***/ for( ii=0 ; ii < ts_length ; ii++ ) ts_array[ii] = 0.0 ; /*** initialize the impulse response ***/ if( nid < ts_length ){ /* make some space for it */ if( fid != NULL ) free(fid) ; nid = ts_length ; fid = (float *) malloc( sizeof(float) * nid ) ; } gamma_model( gs , ts_length , x_array , fid ) ; /* compute impulse */ top = 0.0 ; /* find max value */ for( jj=0 ; jj < ts_length ; jj++ ){ val = fabs(fid[jj]) ; if( val > top ) top = val ; } if( top == 0.0 ) fid[0] = 1.0 ; /* very unlikely case */ top *= 0.001 ; for( jj=0 ; jj < ts_length ; jj++ ){ /* discard small values */ if( fabs(fid[jj]) < top ) fid[jj] = 0.0 ; } for( nid_bot=0 ; nid_bot < ts_length ; nid_bot++ ) /* find first nonzero */ if( fid[nid_bot] != 0.0 ) break ; for( nid_top=ts_length-1 ; nid_top > nid_bot ; nid_top-- ) /* and last nonzero */ if( fid[nid_top] != 0.0 ) break ; /*** loop over each nonzero point in the reference ***/ for( ii=0 ; ii < refnz ; ii++ ){ kk = refin[ii] ; if( kk >= ts_length ) break ; val = refts[kk] ; /*** for each point in the impulse ***/ jtop = ts_length - kk ; if( jtop > nid_top ) jtop = nid_top+1 ; for( jj=nid_bot ; jj < jtop ; jj++ ) ts_array[kk+jj] += val * fid[jj] ; } return ; } /*-----------------------------------------------------------------------*/ DEFINE_MODEL_PROTOTYPE MODEL_interface * initialize_model () { MODEL_interface * mi = NULL; /*----- allocate memory space for model interface -----*/ mi = (MODEL_interface *) XtMalloc (sizeof(MODEL_interface)); /*----- name of this model -----*/ strcpy (mi->label, "ConvGamma"); /*----- this is a signal model -----*/ mi->model_type = MODEL_SIGNAL_TYPE; /*----- number of parameters in the model -----*/ mi->params = 4; /*----- parameter labels -----*/ strcpy (mi->plabel[0], "t0"); strcpy (mi->plabel[1], "amp"); strcpy (mi->plabel[2], "r"); strcpy (mi->plabel[3], "b"); /*----- minimum and maximum parameter constraints -----*/ mi->min_constr[0] = 0.0; mi->max_constr[0] = 10.0; mi->min_constr[1] = 0.0; mi->max_constr[1] = 200.0; mi->min_constr[2] = 1.0; mi->max_constr[2] = 19.0; mi->min_constr[3] = 0.1; mi->max_constr[3] = 5.0; /*----- function which implements the model -----*/ mi->call_func = &conv_model; return (mi); } /*----------------------------------------------------------------------*/ /* Routine to calculate the time series which results from using the gamma variate drug response signal model with the specified model parameters. Definition of model parameters: gs[0] = time delay of response (t0) gs[1] = multiplicative constant (k) gs[2] = rise rate exponent (r) gs[3] = decay rate constant (b) */ void gamma_model ( float * gs, /* parameters for signal model */ int ts_length, /* length of time series data */ float ** x_array, /* independent variable matrix */ float * ts_array /* estimated signal model time series */ ) { int it; /* time index */ float t; /* time */ float gsi , fac ; if( gs[3] <= 0.0 || gs[2] <= 0.0 ){ for( it=0 ; it < ts_length ; it++ ) ts_array[it] = 0.0 ; return ; } /* fac is chosen to make the peak value equal to gs[1] */ gsi = 1.0 / gs[3] ; fac = gs[1] * exp( gs[2] * ( 1.0 - log(gs[2]*gs[3]) ) ) ; for( it=0; it < ts_length; it++){ t = x_array[it][1] - gs[0] ; ts_array[it] = (t <= 0.0) ? 0.0 : fac * exp( log(t) * gs[2] - t * gsi ) ; } return ; }