strains_pomp.c

Go to the documentation of this file.

#include "pomplink.h"
#include "internal.h"
 
static const int nrate = 6;
 
#define STRAIN1 1
#define STRAIN2 2
#define STRAIN3 3
 
#define Beta1     (__p[__parindex[0]])
#define Beta2     (__p[__parindex[1]])
#define Beta3     (__p[__parindex[2]])
#define gamma     (__p[__parindex[3]])
#define chi       (__p[__parindex[4]])
#define S_0       (__p[__parindex[5]])
#define I1_0      (__p[__parindex[6]])
#define I2_0      (__p[__parindex[7]])
#define I3_0      (__p[__parindex[8]])
#define R_0       (__p[__parindex[9]])
#define POP       (__p[__parindex[10]])
#define S         (__x[__stateindex[0]])
#define I1        (__x[__stateindex[1]])
#define I2        (__x[__stateindex[2]])
#define I3        (__x[__stateindex[3]])
#define R         (__x[__stateindex[4]])
#define ll        (__x[__stateindex[5]])
#define ellI1     (__x[__stateindex[6]])
#define ellI2     (__x[__stateindex[7]])
#define ellI3     (__x[__stateindex[8]])
#define node      (__x[__stateindex[9]])
 
#define EVENT_RATES                                     \
  event_rates(__x,__p,t,                                \
              __stateindex,__parindex,__covindex,       \
              __covars,rate,logpi,&penalty)             \
 
static double event_rates
(
 double *__x,
 const double *__p,
 double t,
 const int *__stateindex,
 const int *__parindex,
 const int *__covindex,
 const double *__covars,
 double *rate,
 double *logpi,
 double *penalty
 ) {
  double event_rate = 0;
  double alpha, pi;
  *penalty = 0;
  assert(S >= 0);
  assert(I1 >= ellI1);
  assert(I2 >= ellI2);
  assert(I3 >= ellI3);
  assert(ellI1 >= 0);
  assert(ellI2 >= 0);
  assert(ellI3 >= 0);
  // 0: strain-1 transmission with saturation 0 or 1
  alpha = Beta1*S*I1/POP;
  pi = (I1 > 0) ? 1-ellI1*(ellI1-1)/I1/(I1+1) : 0;
  event_rate += (*rate = alpha*pi); rate++;
  *logpi = 0; logpi++;      // this cancels the boost exactly
  *penalty += alpha*(1-pi);
  // 1: strain 1 recovery
  alpha = gamma*I1;
  pi = (I1 > ellI1) ? 1 : 0;
  event_rate += (*rate = alpha*pi); rate++;
  *logpi = log(pi); logpi++;
  *penalty += alpha*(1-pi);
  // strain 1 sampling
  alpha = chi*I1;
  *penalty += alpha;
  // 2: strain-2 transmission with saturation 0 or 1
  alpha = Beta2*S*I2/POP;
  pi = (I2 > 0) ? 1-ellI2*(ellI2-1)/I2/(I2+1) : 0;
  event_rate += (*rate = alpha*pi); rate++;
  *logpi = 0; logpi++;      // this cancels the boost exactly
  *penalty += alpha*(1-pi);
  // 3: strain 2 recovery
  alpha = gamma*I2;
  pi = (I2 > ellI2) ? 1 : 0;
  event_rate += (*rate = alpha*pi); rate++;
  *logpi = log(pi); logpi++;
  *penalty += alpha*(1-pi);
  // strain 2 sampling
  alpha = chi*I2;
  *penalty += alpha;
  // 4: strain-3 transmission with saturation 0 or 1
  alpha = Beta3*S*I3/POP;
  pi = (I3 > 0) ? 1-ellI3*(ellI3-1)/I3/(I3+1) : 0;
  event_rate += (*rate = alpha*pi); rate++;
  *logpi = 0; logpi++;      // this cancels the boost exactly
  *penalty += alpha*(1-pi);
  // 5: strain 3 recovery
  alpha = gamma*I3;
  pi = (I3 > ellI3) ? 1 : 0;
  event_rate += (*rate = alpha*pi); rate++;
  *logpi = log(pi); logpi++;
  *penalty += alpha*(1-pi);
  // strain 3 sampling
  alpha = chi*I3;
  *penalty += alpha;
  assert(R_FINITE(event_rate));
  return event_rate;
}

void strains_rinit
(
 double *__x,
 const double *__p,
 double t,
 const int *__stateindex,
 const int *__parindex,
 const int *__covindex,
 const double *__covars
 ){
  double m = POP/(S_0+I1_0+I2_0+I3_0+R_0);
  S = nearbyint(S_0*m);
  I1 = nearbyint(I1_0*m);
  I2 = nearbyint(I2_0*m);
  I3 = nearbyint(I3_0*m);
  R = nearbyint(R_0*m);
  ll = 0;
  ellI1 = 0;
  ellI2 = 0;
  ellI3 = 0;
  node = 0;
}

void strains_gill
(
 double *__x,
 const double *__p,
 const int *__stateindex,
 const int *__parindex,
 const int *__covindex,
 const double *__covars,
 double t,
 double dt
 ){
  double tstep = 0.0, tmax = t + dt;
  const int *nodetype = get_userdata_int("nodetype");
  const int *sat = get_userdata_int("saturation");
  const int *deme = get_userdata_int("deme");
  const int *index = get_userdata_int("index");
  const int *child = get_userdata_int("child");
 
  int parent = (int) nearbyint(node);
  int c = child[index[parent]];
 
#ifndef NDEBUG
  const int *lineage = get_userdata_int("lineage");
  int nnode = *get_userdata_int("nnode");
  assert(parent>=0);
  assert(parent<=nnode);
#endif
 
  ll = 0;
 
  // singular portion of filter equation
  switch (nodetype[parent]) {
  default:                      // non-genealogical event #nocov
    break;                      // #nocov
  case 0:                       // root
    assert(sat[parent]==1);
    assert(lineage[parent]==lineage[c]);
    switch (deme[c]) {
    case STRAIN1:
      ellI1 += 1; break;
    case STRAIN2:
      ellI2 += 1; break;
    case STRAIN3:
      ellI3 += 1; break;
    default:                    // #nocov
      assert(0); break;         // #nocov
    }
    break;
  case 1:                       // sample
    assert(sat[parent]==0);
    switch (deme[parent]) {
    case STRAIN1:
      assert(I1 >= ellI1);
      assert(ellI1 >= 0);
      ll += log(chi*I1);
      ellI1 -= 1; I1 -= 1;
      break;
    case STRAIN2:
      assert(I2 >= ellI2);
      assert(ellI2 >= 0);
      ll += log(chi*I2);
      ellI2 -= 1; I2 -= 1;
      break;
    case STRAIN3:
      assert(I3 >= ellI3);
      assert(ellI3 >= 0);
      ll += log(chi*I3);
      ellI3 -= 1; I3 -= 1;
      break;
    default:                    // #nocov
      assert(0); break;         // #nocov
    }
    break;
  case 2:                       // branch point s=(1,1)
    assert(S >= 0);
    if (sat[parent]!=2) break;
    switch (deme[parent]) {
    case STRAIN1:
      assert(I1 >= 0);
      assert(ellI1 > 0);
      ll += (I1 > 0 && I1 >= ellI1) ? log(Beta1*S*I1/POP) : R_NegInf;
      S -= 1; I1 += 1; ellI1 += 1;
      ll -= log(I1*(I1-1)/2);
      break;
    case STRAIN2:
      assert(I2 >= 0);
      assert(ellI2 > 0);
      ll += (I2 > 0 && I2 >= ellI2) ? log(Beta2*S*I2/POP) : R_NegInf;
      S -= 1; I2 += 1; ellI2 += 1;
      ll -= log(I2*(I2-1)/2);
      break;
    case STRAIN3:
      assert(I3 >= 0);
      assert(ellI3 > 0);
      ll += (I3 > 0 && I3 >= ellI3) ? log(Beta3*S*I3/POP) : R_NegInf;
      S -= 1; I3 += 1; ellI3 += 1;
      ll -= log(I3*(I3-1)/2);
      break;
    default:                    // #nocov
      assert(0); break;         // #nocov
    }
    S = (S > 0) ? S : 0;
    break;
  }
 
  if (tmax > t && R_FINITE(ll)) {
 
    // take Gillespie steps to the end of the interval:
    int event;
    double penalty = 0;
    double rate[nrate], logpi[nrate];
 
    double event_rate = EVENT_RATES;
    tstep = exp_rand()/event_rate;
 
    while (t + tstep < tmax) {
      event = rcateg(event_rate,rate,nrate);
      assert(event>=0 && event<nrate);
      ll -= penalty*tstep + logpi[event];
      switch (event) {
      case 0:                     // strain-1 transmission
        S -= 1; I1 += 1;
        break;
      case 1:                     // strain-1 recovery
        I1 -= 1; R += 1;
        break;
      case 2:                     // strain-2 transmission
        S -= 1; I2 += 1;
        break;
      case 3:                     // strain-2 recovery
        I2 -= 1; R += 1;
        break;
      case 4:                     // strain-3 transmission
        S -= 1; I3 += 1;
        break;
      case 5:                     // strain-3 recovery
        I3 -= 1; R += 1;
        break;
      default:                  // #nocov
        assert(0); break;       // #nocov
      }
      t += tstep;
      event_rate = EVENT_RATES;
      tstep = exp_rand()/event_rate;
    }
    tstep = tmax - t;
    ll -= penalty*tstep;
  }
  node += 1;
}
 
# define lik  (__lik[0])

void strains_dmeas
(
 double *__lik,
 const double *__y,
 const double *__x,
 const double *__p,
 int give_log,
 const int *__obsindex,
 const int *__stateindex,
 const int *__parindex,
 const int *__covindex,
 const double *__covars,
 double t
 ){
  assert(!ISNAN(ll));
  lik = (give_log) ? ll : exp(ll);
}