/*****************************************************************************
   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"

/****************************************************************************
     Model function for 3dNLfim:
     gamma variate impulse response, convolved with NREF input time series,
     each with an independent scaling amplitude, but all other parameters
     of the gamma variate fixed.  There are NREF+3 parameters:
       t0    = time delay
       r     = power
       b     = time scale
       amp_j = amplitude of j-th reference, for j=1..NREF

     28 July 1998 -- RWCox
****************************************************************************/

#define NREF 2

static int     refnum[NREF] ;  /* # pts in refts */
static int     refnz[NREF]  ;  /* # of nonzero pts */
static float * refts[NREF]  ;  /* reference time series */
static int   * refin[NREF]  ;  /* 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.
     num         = length of input time series
     ref[jv][ii] = ii-th point in jv-th time series,
                   for ii=0..num-1, jv=0..NREF-1
   If num==0, will read the reference time series from a file.
------------------------------------------------------------------------*/

void conv_set_ref( int num , float ** ref )
{
   int jv ;

   if( num > 0 && ref != NULL ){ /*** if have inputs, make space & copy in ***/
      int ii ;

      for( jv=0 ; jv < NREF ; jv++ ){

         /* get rid of old data? */

         if( refts[jv] != NULL ){
            free(refts[jv]); refts[jv] = NULL;
            free(refin[jv]); refin[jv] = NULL;
         }

         /* copy new data */

         refnum[jv] = num ;
         refts[jv]  = (float *) malloc( sizeof(float) * num ) ;
         refin[jv]  = (int *)   malloc( sizeof(int)   * num ) ;
         memcpy( refts[jv] , ref[jv] , sizeof(float) * num ) ;

         /* build a list of nonzero entries in this column */

         for( ii=0,refnz[jv]=0 ; ii < num ; ii++ )
            if( refts[jv][ii] != 0.0 ) refin[jv][refnz[jv]++] = ii ;

         if( refnz[jv] == 0 )
            ERREX(__FILE__ ": All zero reference timeseries column!") ;
      }
      return ;

   } else { /*** if no inputs, do something special ***/

     char * cp ;
     MRI_IMAGE * flim ;
     int jv , nx ;
     float * ref[NREF] ;

     cp = my_getenv("AFNI_CONVMODEL_REF") ;  /* get name of reference file */
     if( cp == NULL )
        ERREX(__FILE__ ": 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,__FILE__ ": Can't read timeseries file %s",cp) ;
        ERREX(buf) ;
     } else {
        fprintf(stderr,__FILE__ ": Read reference file %s\n",cp) ;
     }

     if( flim->ny < NREF )
        ERREX(__FILE__ ": reference file has too few columns!") ;
     else if( flim->nv > NREF )
        fprintf(stderr,__FILE__ " WARNING: reference file has too many columns!\n") ;

     nx = flim->nx ;
     for( jv=0 ; jv < NREF ; jv++ )
        ref[jv] = MRI_FLOAT_PTR(flim) + jv*nx ;

     conv_set_ref( nx , ref ) ;  /* recursion! */
     mri_free(flim) ;
   }
   return ;
}

/*-----------------------------------------------------------------------
  Function to compute the simulated time series:
    gs[0..2]      = gamma variate parameters t0, r, b
    gs[3..NREF+3] = amplitudes for each reference 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 , jv ;
   float top , val , amp ;

   static int     nid = 0 ;     /* number of pts in fid array */
   static float * fid = NULL ;  /* impulse response function */

   /*** make sure there is a reference function to convolve with ***/

   if( refnum[0] <= 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 */

   /*** discard small values from the fid ***/

#define EPS 0.001  /* definition of small */

#undef FIND_TOP  /* don't need this, since our fid amplitude is set to 1.0 */
#ifdef FIND_TOP
   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 *= EPS ;
#else
   top = EPS ;
#endif
   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 reference ***/

   for( jv=0 ; jv < NREF ; jv++ ){

      amp = gs[jv+3] ; if( amp == 0.0 ) continue ;

      /*** loop over each nonzero point in the reference ***/

      for( ii=0 ; ii < refnz[jv] ; ii++ ){
         kk  = refin[jv][ii] ; if( kk >= ts_length ) break ;
         val = amp * refts[jv][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;
  int jv ;
  char buf[32] ;

  /*----- allocate memory space for model interface -----*/

  mi = (MODEL_interface *) RwcMalloc (sizeof(MODEL_interface));

  /*----- name of this model -----*/

  sprintf( mi->label , "ConvGamma%da" , NREF ) ;

  /*----- this is a signal model -----*/

  mi->model_type = MODEL_SIGNAL_TYPE;

  /*----- number of parameters in the model -----*/

  mi->params = 3 + NREF ;

  /*----- parameter labels -----*/

  strcpy (mi->plabel[0], "t0");  /* impulse response is proportional to */
  strcpy (mi->plabel[1], "r");   /*       r  -(t-t0)/b  */
  strcpy (mi->plabel[2], "b");   /* (t-t0)  e           */

  for( jv=0 ; jv < NREF ; jv++ )
     sprintf( mi->plabel[jv+3] , "amp_%d" , jv+1 ) ;

  /*----- minimum and maximum parameter constraints -----*/

  mi->min_constr[0] = 0.0;  mi->max_constr[0] = 10.0;  /* delay */
  mi->min_constr[1] = 1.0;  mi->max_constr[1] = 19.0;  /* power */
  mi->min_constr[2] = 0.1;  mi->max_constr[2] =  5.0;  /* time scale */

  for( jv=0 ; jv < NREF ; jv++ ){
     mi->min_constr[jv+3] = 0.0;  mi->max_constr[jv+3] = 200.0;
  }

  /*----- function which implements the model -----*/
  mi->call_func = (void_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] = rise rate exponent (r)
	 gs[2] = decay rate constant (b)

  Time to peak is r * b ;
  FWHM of the peak is about 2.3 * sqrt(r) * b, for r > 1.

*/

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[2] <= 0.0 || gs[1] <= 0.0 ){
     ts_array[0] = 1.0 ;
     for( it=1 ; it < ts_length ; it++ ) ts_array[it] = 0.0 ;
     return ;
  }

  /* fac is chosen to make the peak value equal to 1 (at t = gs[1]*gs[2]) */

  gsi = 1.0 / gs[2] ;
  fac = exp( gs[1] * ( 1.0 - log(gs[1]*gs[2]) ) ) ;
  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[1] - t * gsi ) ;
  }
  return ;
}
