usr/kevin/code/develop/gillespie.c
int gillespie_do_step(state_t * p) {

  double tstep;
  double rs;
  double cumu_a;
  int i, j;
  int ifail = 0;

  /* determine propensity functions */

  dyn.sum_a = 0.0;

  for (i = 0; i < dyn.Nreact; i++) {

    if (dyn.R[i].Nreact == 0) { 
      dyn.a[i] = dyn.R[i].k;
    }
    else if (dyn.R[i].Nreact == 1) { 
      dyn.a[i] = dyn.R[i].k * p->nx[dyn.R[i].react[0].index];
    }
    else if (dyn.R[i].react[0].index == dyn.R[i].react[1].index) {
      dyn.a[i] = dyn.R[i].k * p->nx[dyn.R[i].react[0].index]
	* (p->nx[dyn.R[i].react[1].index] - 1);
    }
    else {
      dyn.a[i] = dyn.R[i].k * p->nx[dyn.R[i].react[0].index]
	* p->nx[dyn.R[i].react[1].index];
    }

    dyn.sum_a += dyn.a[i];
  }

  /* propagate time  */

  if (dyn.sum_a < FLT_EPSILON) { 
    /* No reactions are possible. */
    ifail = 1;
  }
  else {

    /* We rely here on the fact that ranlcg does not produce zero. */
    rs = ranlcg_reep(random_generator_);
    tstep = log(1./rs)/dyn.sum_a;

    p->t += tstep;

    /* select reaction */

    rs = ranlcg_reep(random_generator_); 
    rs *= dyn.sum_a;

    j = 0;
    cumu_a = dyn.a[j];
    while (cumu_a < rs) cumu_a += dyn.a[++j];

    /* update concentrations */

    for (i = 0; i < dyn.R[j].Nreact; i++) {
      p->nx[dyn.R[j].react[i].index] --;
    }

    for (i = 0; i < dyn.R[j].Nprod; i++) {
      p->nx[dyn.R[j].prod[i].index] += dyn.R[j].prod[i].change;
    }
  }

  return ifail;
}
/*****************************************************************************
 *
 *  gillespie_do_step
 *
 *  Advance state p by one step.
 *
 *****************************************************************************/
gillespie_do_step
Advance state p by one step.
gillespie_do_step Advance state p by one step .
int gillespie_do_step(state_t * p) {
  double tstep;
  double rs;
  double cumu_a;
  int i, j;
  int ifail = 0;
  dyn.sum_a = 0.0;
  for (i = 0; i < dyn.Nreact; i++) {
    if (dyn.R[i].Nreact == 0) { 
      dyn.a[i] = dyn.R[i].k;
    }
    else if (dyn.R[i].Nreact == 1) { 
      dyn.a[i] = dyn.R[i].k * p->nx[dyn.R[i].react[0].index];
    }
    else if (dyn.R[i].react[0].index == dyn.R[i].react[1].index) {
      dyn.a[i] = dyn.R[i].k * p->nx[dyn.R[i].react[0].index]
	* (p->nx[dyn.R[i].react[1].index] - 1);
    }
    else {
      dyn.a[i] = dyn.R[i].k * p->nx[dyn.R[i].react[0].index]
	* p->nx[dyn.R[i].react[1].index];
    }
    dyn.sum_a += dyn.a[i];
  }
  if (dyn.sum_a < FLT_EPSILON) { 
    ifail = 1;
  }
  else {
    rs = ranlcg_reep(random_generator_);
    tstep = log(1./rs)/dyn.sum_a;
    p->t += tstep;
    rs = ranlcg_reep(random_generator_); 
    rs *= dyn.sum_a;
    j = 0;
    cumu_a = dyn.a[j];
    while (cumu_a < rs) cumu_a += dyn.a[++j];
    for (i = 0; i < dyn.R[j].Nreact; i++) {
      p->nx[dyn.R[j].react[i].index] --;
    }
    for (i = 0; i < dyn.R[j].Nprod; i++) {
      p->nx[dyn.R[j].prod[i].index] += dyn.R[j].prod[i].change;
    }
  }
  return ifail;
}
int gillespie_do_step ( state_t * p ) { double tstep ; double rs ; double cumu_a ; int i , j ; int ifail = 0 ; dyn . sum_a = 0.0 ; for ( i = 0 ; i < dyn . Nreact ; i ++ ) { if ( dyn . R [ i ] . Nreact == 0 ) { dyn . a [ i ] = dyn . R [ i ] . k ; } else if ( dyn . R [ i ] . Nreact == 1 ) { dyn . a [ i ] = dyn . R [ i ] . k * p -> nx [ dyn . R [ i ] . react [ 0 ] . index ] ; } else if ( dyn . R [ i ] . react [ 0 ] . index == dyn . R [ i ] . react [ 1 ] . index ) { dyn . a [ i ] = dyn . R [ i ] . k * p -> nx [ dyn . R [ i ] . react [ 0 ] . index ] * ( p -> nx [ dyn . R [ i ] . react [ 1 ] . index ] - 1 ) ; } else { dyn . a [ i ] = dyn . R [ i ] . k * p -> nx [ dyn . R [ i ] . react [ 0 ] . index ] * p -> nx [ dyn . R [ i ] . react [ 1 ] . index ] ; } dyn . sum_a += dyn . a [ i ] ; } if ( dyn . sum_a < FLT_EPSILON ) { ifail = 1 ; } else { rs = ranlcg_reep ( random_generator_ ) ; tstep = log ( 1. / rs ) / dyn . sum_a ; p -> t += tstep ; rs = ranlcg_reep ( random_generator_ ) ; rs *= dyn . sum_a ; j = 0 ; cumu_a = dyn . a [ j ] ; while ( cumu_a < rs ) cumu_a += dyn . a [ ++ j ] ; for ( i = 0 ; i < dyn . R [ j ] . Nreact ; i ++ ) { p -> nx [ dyn . R [ j ] . react [ i ] . index ] -- ; } for ( i = 0 ; i < dyn . R [ j ] . Nprod ; i ++ ) { p -> nx [ dyn . R [ j ] . prod [ i ] . index ] += dyn . R [ j ] . prod [ i ] . change ; } } return ifail ; }
gillespie_do_step
  Advance state p by one step.
gillespie_do_step Advance state p by one step .
/* determine propensity functions */ /* propagate time  */ /* No reactions are possible. */ /* We rely here on the fact that ranlcg does not produce zero. */ /* select reaction */ /* update concentrations */