src/sim/sim_dmc.c
int dmc_do_step(dynam_t * dyn) {

  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->nreactions; i++) {

    if (dyn->R[i].nreactant == 0) { 
      dyn->a[i] = dyn->R[i].k;
    }
    else if (dyn->R[i].nreactant == 1) { 
      dyn->a[i] = dyn->R[i].k * dyn->state.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 * dyn->state.nx[dyn->R[i].react[0].index]
	* (dyn->state.nx[dyn->R[i].react[1].index] - 1);
    }
    else {
      dyn->a[i] = dyn->R[i].k * dyn->state.nx[dyn->R[i].react[0].index]
	* dyn->state.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. */

    ranlcg_reep(dyn->rng, &rs);
    tstep = log(1./rs)/dyn->sum_a;

    dyn->state.t += tstep;

    /* select reaction */

    ranlcg_reep(dyn->rng, &rs); 
    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].nreactant; i++) {
      dyn->state.nx[dyn->R[j].react[i].index] --;
    }

    for (i = 0; i < dyn->R[j].nproduct; i++) {
      dyn->state.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 dmc_do_step(dynam_t * dyn) {
  double tstep;
  double rs;
  double cumu_a;
  int i, j;
  int ifail = 0;
  dyn->sum_a = 0.0;
  for (i = 0; i < dyn->nreactions; i++) {
    if (dyn->R[i].nreactant == 0) { 
      dyn->a[i] = dyn->R[i].k;
    }
    else if (dyn->R[i].nreactant == 1) { 
      dyn->a[i] = dyn->R[i].k * dyn->state.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 * dyn->state.nx[dyn->R[i].react[0].index]
	* (dyn->state.nx[dyn->R[i].react[1].index] - 1);
    }
    else {
      dyn->a[i] = dyn->R[i].k * dyn->state.nx[dyn->R[i].react[0].index]
	* dyn->state.nx[dyn->R[i].react[1].index];
    }
    dyn->sum_a += dyn->a[i];
  }
  if (dyn->sum_a < FLT_EPSILON) { 
    ifail = 1;
  }
  else {
    ranlcg_reep(dyn->rng, &rs);
    tstep = log(1./rs)/dyn->sum_a;
    dyn->state.t += tstep;
    ranlcg_reep(dyn->rng, &rs); 
    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].nreactant; i++) {
      dyn->state.nx[dyn->R[j].react[i].index] --;
    }
    for (i = 0; i < dyn->R[j].nproduct; i++) {
      dyn->state.nx[dyn->R[j].prod[i].index] += dyn->R[j].prod[i].change;
    }
  }
  return ifail;
}
int dmc_do_step ( dynam_t * dyn ) { double tstep ; double rs ; double cumu_a ; int i , j ; int ifail = 0 ; dyn -> sum_a = 0.0 ; for ( i = 0 ; i < dyn -> nreactions ; i ++ ) { if ( dyn -> R [ i ] . nreactant == 0 ) { dyn -> a [ i ] = dyn -> R [ i ] . k ; } else if ( dyn -> R [ i ] . nreactant == 1 ) { dyn -> a [ i ] = dyn -> R [ i ] . k * dyn -> state . 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 * dyn -> state . nx [ dyn -> R [ i ] . react [ 0 ] . index ] * ( dyn -> state . nx [ dyn -> R [ i ] . react [ 1 ] . index ] - 1 ) ; } else { dyn -> a [ i ] = dyn -> R [ i ] . k * dyn -> state . nx [ dyn -> R [ i ] . react [ 0 ] . index ] * dyn -> state . nx [ dyn -> R [ i ] . react [ 1 ] . index ] ; } dyn -> sum_a += dyn -> a [ i ] ; } if ( dyn -> sum_a < FLT_EPSILON ) { ifail = 1 ; } else { ranlcg_reep ( dyn -> rng , & rs ) ; tstep = log ( 1. / rs ) / dyn -> sum_a ; dyn -> state . t += tstep ; ranlcg_reep ( dyn -> rng , & rs ) ; 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 ] . nreactant ; i ++ ) { dyn -> state . nx [ dyn -> R [ j ] . react [ i ] . index ] -- ; } for ( i = 0 ; i < dyn -> R [ j ] . nproduct ; i ++ ) { dyn -> state . 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 */