#include "pomplink.h"
#include "internal.h"

Include dependency graph for sirs_pomp.c:

Macros
#define	Beta (__p[__parindex[0]])

#define	gamma (__p[__parindex[1]])

#define	psi (__p[__parindex[2]])

#define	chi (__p[__parindex[3]])

#define	omega (__p[__parindex[4]])

#define	S0 (__p[__parindex[5]])

#define	I0 (__p[__parindex[6]])

#define	R0 (__p[__parindex[7]])

#define	POP (__p[__parindex[8]])

#define	S (__x[__stateindex[0]])

#define	I (__x[__stateindex[1]])

#define	R (__x[__stateindex[2]])

#define	ll (__x[__stateindex[3]])

#define	ellI (__x[__stateindex[4]])

#define	node (__x[__stateindex[5]])

#define	EVENT_RATES

#define	lik (__lik[0])

Functions
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 penalty)

void	sirs_rinit (double __x, const double __p, double t, const int __stateindex, const int __parindex, const int __covindex, const double __covars)
	Latent-state initializer (rinit).

void	sirs_gill (double __x, const double __p, const int __stateindex, const int __parindex, const int __covindex, const double __covars, double t, double dt)

void	sirs_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)
	Measurement model likelihood (dmeasure).

Variables
static const int	nrate = 3

Macro Definition Documentation

◆ Beta

#define Beta (__p[__parindex[0]])

Definition at line 6 of file sirs_pomp.c.

◆ chi

#define chi (__p[__parindex[3]])

Definition at line 9 of file sirs_pomp.c.

◆ ellI

#define ellI (__x[__stateindex[4]])

Definition at line 19 of file sirs_pomp.c.

◆ EVENT_RATES

#define EVENT_RATES

Value:

  event_rates(__x,__p,t,                                \
              __stateindex,__parindex,__covindex,       \
              __covars,rate,&penalty)                   \

Definition at line 22 of file sirs_pomp.c.

#define EVENT_RATES                                     \
  event_rates(__x,__p,t,                                \
              __stateindex,__parindex,__covindex,       \
              __covars,rate,&penalty)                   \
 

◆ gamma

#define gamma (__p[__parindex[1]])

Definition at line 7 of file sirs_pomp.c.

◆ I

#define I (__x[__stateindex[1]])

Definition at line 16 of file sirs_pomp.c.

◆ I0

#define I0 (__p[__parindex[6]])

Definition at line 12 of file sirs_pomp.c.

◆ lik

#define lik (__lik[0])

Definition at line 195 of file sirs_pomp.c.

◆ ll

#define ll (__x[__stateindex[3]])

Definition at line 18 of file sirs_pomp.c.

◆ node

#define node (__x[__stateindex[5]])

Definition at line 20 of file sirs_pomp.c.

◆ omega

#define omega (__p[__parindex[4]])

Definition at line 10 of file sirs_pomp.c.

◆ POP

#define POP (__p[__parindex[8]])

Definition at line 14 of file sirs_pomp.c.

◆ psi

#define psi (__p[__parindex[2]])

Definition at line 8 of file sirs_pomp.c.

◆ R

#define R (__x[__stateindex[2]])

Definition at line 17 of file sirs_pomp.c.

◆ R0

#define R0 (__p[__parindex[7]])

Definition at line 13 of file sirs_pomp.c.

◆ S

#define S (__x[__stateindex[0]])

Definition at line 15 of file sirs_pomp.c.

◆ S0

#define S0 (__p[__parindex[5]])

Definition at line 11 of file sirs_pomp.c.

Function Documentation

◆ event_rates()

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 *	penalty )

static

Definition at line 27 of file sirs_pomp.c.

   {
  double event_rate = 0;
  double alpha, disc;
  *penalty = 0;
  assert(I >= ellI);
  assert(ellI >= 0);
  assert(S >= 0);
  // 0: transmission with saturation 0 or 1
  alpha = Beta*S*I/POP;
  disc = (I > 0) ? ellI*(ellI-1)/I/(I+1) : 1;
  event_rate += (*rate = alpha*(1-disc)); rate++;
  *penalty += alpha*disc;
  // 1: recovery
  alpha = gamma*I;
  if (I > ellI) {
    event_rate += (*rate = alpha); rate++;
  } else {
    *rate = 0; rate++;
    *penalty += alpha;
  }
  // 2: loss of immunity
  alpha = omega*R;
  event_rate += (*rate = alpha); rate++;
  // sampling
  alpha = (psi+chi)*I;
  *penalty += alpha;
  assert(R_FINITE(event_rate));
  return event_rate;
}

◆ sirs_dmeas()

void sirs_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 )

Measurement model likelihood (dmeasure).

Definition at line 198 of file sirs_pomp.c.

  {
  assert(!ISNAN(ll));
  lik = (give_log) ? ll : exp(ll);
}

◆ sirs_gill()

void sirs_gill	(	double *	__x,
		const double *	__p,
		const int *	__stateindex,
		const int *	__parindex,
		const int *	__covindex,
		const double *	__covars,
		double	t,
		double	dt )

Latent-state process simulator (rprocess).

This integrates the filter equation.

Definition at line 91 of file sirs_pomp.c.

  {
  double tstep = 0.0, tmax = t + dt;
  const int *nodetype = get_userdata_int("nodetype");
  const int *sat = get_userdata_int("saturation");
 
  int parent = (int) nearbyint(node);
 
#ifndef NDEBUG
  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
    ellI += 1;
    break;
  case 1:                       // sample
    assert(I >= ellI);
    assert(ellI >= 0);
    if (sat[parent] == 1) {
      ll += log(psi);
    } else if (sat[parent] == 0) {
      ellI -= 1;
      ll += log(psi+chi);
      if (psi+chi <= 0)
    ll += R_NegInf;
      else if (unif_rand() < psi/(psi+chi)) { // non-destructive sample
    ll += log(I-ellI);
      } else {          // destructive sample
    ll += log(I);
    I -= 1;
      }
    } else {
      assert(0);                // #nocov
      ll += R_NegInf;           // #nocov
    }
    break;
  case 2:                       // branch point s=(1,1)
    assert(S >= 0);
    assert(I >= 0);
    assert(ellI > 0);
    assert(sat[parent]==2);
    ll += (I > 0 && I >= ellI) ? log(Beta*S*I/POP) : R_NegInf;
    S -= 1; I += 1;
    ellI += 1;
    ll -= log(I*(I-1)/2);
    S = (S > 0) ? S : 0;
    break;
  }
 
  if (tmax > t) {
 
    // take Gillespie steps to the end of the interval:
    int event;
    double penalty = 0;
    double rate[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;
      switch (event) {
      case 0:                     // transmission
        S -= 1; I += 1;
        break;
      case 1:                     // recovery
        I -= 1; R += 1;
        break;
      case 2:                     // loss of immunity
        R -= 1; S += 1;
        break;
      default:                    // #nocov
        assert(0);                // #nocov
        break;                    // #nocov
      }
      t += tstep;
      event_rate = EVENT_RATES;
      tstep = exp_rand()/event_rate;
    }
    tstep = tmax - t;
    ll -= penalty*tstep;
  }
  node += 1;
}

Here is the call graph for this function:

◆ sirs_rinit()

void sirs_rinit	(	double *	__x,
		const double *	__p,
		double	t,
		const int *	__stateindex,
		const int *	__parindex,
		const int *	__covindex,
		const double *	__covars )

Latent-state initializer (rinit).

Definition at line 69 of file sirs_pomp.c.

  {
  double m = POP/(S0+I0+R0);
  S = nearbyint(S0*m);
  I = nearbyint(I0*m);
  R = nearbyint(R0*m);
  ll = 0;
  ellI = 0;
  node = 0;
}

Variable Documentation

◆ nrate

const int nrate = 3