use deme names instead of numbers in yaml model specification

kingaa · Aug 25, 2024 · 9fef413 · 9fef413
1 parent ea24893
commit 9fef413
Show file tree

Hide file tree

Showing 28 changed files with 749 additions and 173 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: phylopomp
 Type: Package
 Title: Phylodynamic Inference for POMP Models
-Version: 0.14.3.0
-Date: 2024-08-14
+Version: 0.14.4.0
+Date: 2024-08-25
 Authors@R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="kingaa@umich.edu",comment=c(ORCID="0000-0001-6159-3207")),
 	  person(given=c("Qianying"),family="Lin",role=c("aut"),comment=c(ORCID="0000-0001-8620-9910"))
 	  )

diff --git a/src/decls.h b/src/decls.h
@@ -12,3 +12,7 @@ extern void seirs_dmeas(double *, const double *, const double *, const double *
 extern void sirs_rinit(double *, const double *, double, const int *, const int *, const int *, const double *);
 extern void sirs_gill(double *, const double *, const int *, const int *, const int *, const double *, double, double);
 extern void sirs_dmeas(double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double);
+/* src/twospecies_pomp.c */
+extern void twospecies_rinit(double *, const double *, double, const int *, const int *, const int *, const double *);
+extern void twospecies_gill(double *, const double *, const int *, const int *, const int *, const double *, double, double);
+extern void twospecies_dmeas(double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double);
diff --git a/src/lbdp.cc b/src/lbdp.cc
@@ -4,6 +4,8 @@
 #include "generics.h"
 #include "internal.h"
 
+static int deme = 0;
+
 //! LBDP process state.
 typedef struct {
   int n;
@@ -63,7 +65,7 @@ double lbdp_proc_t::event_rates (double *rate, int n) const {
 template<>
 void lbdp_genealogy_t::rinit (void) {
   state.n = params.n0;
-  graft(0,params.n0);
+  graft(deme,params.n0);
 }
 
 template<>

diff --git a/src/moran.cc b/src/moran.cc
@@ -4,6 +4,8 @@
 #include "generics.h"
 #include "internal.h"
 
+static int deme = 0;
+
 //! Moran process state.
 typedef struct {
   int m;
@@ -61,7 +63,7 @@ double moran_proc_t::event_rates (double *rate, int n) const {
 template<>
 void moran_genealogy_t::rinit (void) {
   state.m = state.g = 0;
-  graft(0,params.n);
+  graft(deme,params.n);
 }
 
 template<>

diff --git a/src/s2i2r2.cc b/src/s2i2r2.cc
@@ -4,6 +4,10 @@
 #include "generics.h"
 #include "internal.h"
 
+static int host1 = 0;
+static int host2 = 1;
+static int outside = 2;
+
 //! S2I2R2 process state.
 typedef struct {
   int S1;
@@ -151,27 +155,27 @@ void s2i2r2_genealogy_t::rinit (void) {
   state.R2 = params.R2_0;
   state.N1 = double(params.S1_0+params.I1_0+params.R1_0);
   state.N2 = double(params.S2_0+params.I2_0+params.R2_0);
-  graft(0,params.I1_0);
-  graft(1,params.I2_0);
+  graft(host1,params.I1_0);
+  graft(host2,params.I2_0);
 }
 
 template<>
 void s2i2r2_genealogy_t::jump (int event) {
   switch (event) {
   case 0:
-    state.S1 -= 1; state.I1 += 1; birth(0,0);
+    state.S1 -= 1; state.I1 += 1; birth(host1,host1);
     break;
   case 1:
-    state.S2 -= 1; state.I2 += 1; birth(1,1);
+    state.S2 -= 1; state.I2 += 1; birth(host2,host2);
     break;
   case 2:
-    state.S1 -= 1; state.I1 += 1; birth(1,0);
+    state.S1 -= 1; state.I1 += 1; birth(host2,host1);
     break;
   case 3:
-    state.I1 -= 1; state.R1 += 1; death(0);
+    state.I1 -= 1; state.R1 += 1; death(host1);
     break;
   case 4:
-    state.I2 -= 1; state.R2 += 1; death(1);
+    state.I2 -= 1; state.R2 += 1; death(host2);
     break;
   case 5:
     state.R1 -= 1; state.S1 += 1;
@@ -180,16 +184,16 @@ void s2i2r2_genealogy_t::jump (int event) {
     state.R2 -= 1; state.S2 += 1;
     break;
   case 7:
-    sample(0);
+    sample(host1);
     break;
   case 8:
-    sample(1);
+    sample(host2);
     break;
   case 9:
-    state.S1 -= 1; state.I1 += 1; graft(2); migrate(2,0);
+    state.S1 -= 1; state.I1 += 1; graft(outside); migrate(outside,host1);
     break;
   case 10:
-    state.S2 -= 1; state.I2 += 1; graft(2); migrate(2,1);
+    state.S2 -= 1; state.I2 += 1; graft(outside); migrate(outside,host2);
     break;
   case 11:
     state.S1 -= 1; state.N1 -= 1;
@@ -198,10 +202,10 @@ void s2i2r2_genealogy_t::jump (int event) {
     state.S2 -= 1; state.N2 -= 1;
     break;
   case 13:
-    state.I1 -= 1; state.N1 -= 1; death(0);
+    state.I1 -= 1; state.N1 -= 1; death(host1);
     break;
   case 14:
-    state.I2 -= 1; state.N2 -= 1; death(1);
+    state.I2 -= 1; state.N2 -= 1; death(host2);
     break;
   case 15:
     state.R1 -= 1; state.N1 -= 1;

diff --git a/src/seir.cc b/src/seir.cc
@@ -4,6 +4,9 @@
 #include "generics.h"
 #include "internal.h"
 
+static int Exposed = 0;
+static int Infectious = 1;
+
 //! SEIR process state.
 typedef struct {
   int S;
@@ -100,16 +103,16 @@ template<>
 void seir_genealogy_t::jump (int event) {
   switch (event) {
   case 0:
-    state.S -= 1; state.E += 1; birth(1,0);
+    state.S -= 1; state.E += 1; birth(Infectious,Exposed);
     break;
   case 1:
-    state.E -= 1; state.I += 1; migrate(0,1);
+    state.E -= 1; state.I += 1; migrate(Exposed,Infectious);
     break;
   case 2:
-    state.I -= 1; state.R += 1; death(1);
+    state.I -= 1; state.R += 1; death(Infectious);
     break;
   case 3:
-    sample(1);
+    sample(Infectious);
     break;
   case 4:
     state.R -= 1; state.S += 1;

diff --git a/src/si2r.cc b/src/si2r.cc
@@ -4,6 +4,9 @@
 #include "generics.h"
 #include "internal.h"
 
+static int ordinary = 0;
+static int superspreader = 1;
+
 //! SI2R process state.
 typedef struct {
   int S;
@@ -102,42 +105,44 @@ void si2r_genealogy_t::rinit (void) {
   state.I2 = 0;
   state.R = params.R0;
   state.N = double(params.S0+params.I0+params.R0);
-  graft(0,params.I0);
+  graft(ordinary,params.I0);
 }
 
 template<>
 void si2r_genealogy_t::jump (int event) {
   switch (event) {
   case 0:
-    state.S -= 1; state.I1 += 1; birth(0,0);
+    state.S -= 1; state.I1 += 1; birth(ordinary,ordinary);
     break;
   case 1:
     {
       int n = 1+int(rgeom(1.0/params.mu));
       if (state.S >= n) {
-        state.S -= n; state.I1 += n; birth(1,0,n);
+        state.S -= n; state.I1 += n;
+        birth(superspreader,ordinary,n);
       } else {
-        birth(1,0,state.S); state.I1 += state.S; state.S = 0;
+        birth(superspreader,ordinary,state.S);
+        state.I1 += state.S; state.S = 0;
       }
     }
     break;
   case 2:
-    state.I1 -= 1; state.R += 1; death(0);
+    state.I1 -= 1; state.R += 1; death(ordinary);
     break;
   case 3:
-    state.I2 -= 1; state.R += 1; death(1);
+    state.I2 -= 1; state.R += 1; death(superspreader);
     break;
   case 4:
-    sample(0);
+    sample(ordinary);
     break;
   case 5:
-    sample(1);
+    sample(superspreader);
     break;
   case 6:
-    state.I1 -= 1; state.I2 += 1; migrate(0,1);
+    state.I1 -= 1; state.I2 += 1; migrate(ordinary,superspreader);
     break;
   case 7:
-    state.I1 += 1; state.I2 -= 1; migrate(1,0);
+    state.I1 += 1; state.I2 -= 1; migrate(superspreader,ordinary);
     break;
   case 8:
     state.R -= 1; state.S += 1;

diff --git a/src/siir.cc b/src/siir.cc
@@ -4,6 +4,9 @@
 #include "generics.h"
 #include "internal.h"
 
+static int strain1 = 0;
+static int strain2 = 1;
+
 //! SIIR process state.
 typedef struct {
   int S;
@@ -105,36 +108,36 @@ void siir_genealogy_t::rinit (void) {
   state.I2 = params.I2_0;
   state.R = params.R0;
   state.N = double(params.S0+params.I1_0+params.I2_0+params.R0);
-  graft(0,params.I1_0);
-  graft(1,params.I2_0);
+  graft(strain1,params.I1_0);
+  graft(strain2,params.I2_0);
 }
 
 template<>
 void siir_genealogy_t::jump (int event) {
   switch (event) {
   case 0:
-    state.S -= 1; state.I1 += 1; birth(0,0);
+    state.S -= 1; state.I1 += 1; birth(strain1,strain1);
     break;
   case 1:
-    state.S -= 1; state.I2 += 1; birth(1,1);
+    state.S -= 1; state.I2 += 1; birth(strain2,strain2);
     break;
   case 2:
-    state.I1 -= 1; state.R += 1; death(0);
+    state.I1 -= 1; state.R += 1; death(strain1);
     break;
   case 3:
-    state.I2 -= 1; state.R += 1; death(1);
+    state.I2 -= 1; state.R += 1; death(strain2);
     break;
   case 4:
-    sample(0);
+    sample(strain1);
     break;
   case 5:
-    sample(1);
+    sample(strain2);
     break;
   case 6:
-    state.I1 -= 1; state.I2 += 1; migrate(0,1);
+    state.I1 -= 1; state.I2 += 1; migrate(strain1,strain2);
     break;
   case 7:
-    state.I1 += 1; state.I2 -= 1; migrate(1,0);
+    state.I1 += 1; state.I2 -= 1; migrate(strain2,strain1);
     break;
   case 8:
     state.S += 1; state.R -= 1;

diff --git a/src/sir.cc b/src/sir.cc
@@ -4,6 +4,8 @@
 #include "generics.h"
 #include "internal.h"
 
+static int Infected = 0;
+
 //! SIR process state.
 typedef struct {
   int S;
@@ -82,7 +84,7 @@ void sir_genealogy_t::rinit (void) {
   state.I = params.I0;
   state.R = params.R0;
   state.N = double(params.S0+params.I0+params.R0);
-  graft(0,params.I0);
+  graft(Infected,params.I0);
 }
 
 template<>

diff --git a/src/twospecies.cc b/src/twospecies.cc
@@ -4,6 +4,10 @@
 #include "generics.h"
 #include "internal.h"
 
+static int host1 = 0;
+static int host2 = 1;
+static int outside = 2;
+
 //! TwoSpecies process state.
 typedef struct {
   int S1;
@@ -155,30 +159,30 @@ void twospecies_genealogy_t::rinit (void) {
   state.R2 = params.R2_0;
   state.N1 = double(params.S1_0+params.I1_0+params.R1_0);
   state.N2 = double(params.S2_0+params.I2_0+params.R2_0);
-  graft(0,params.I1_0);
-  graft(1,params.I2_0);
+  graft(host1,params.I1_0);
+  graft(host2,params.I2_0);
 }
 
 template<>
 void twospecies_genealogy_t::jump (int event) {
   switch (event) {
   case 0:
-    state.S1 -= 1; state.I1 += 1; birth(0,0);
+    state.S1 -= 1; state.I1 += 1; birth(host1,host1);
     break;
   case 1:
-    state.S2 -= 1; state.I2 += 1; birth(1,1);
+    state.S2 -= 1; state.I2 += 1; birth(host2,host2);
     break;
   case 2:
-    state.S1 -= 1; state.I1 += 1; birth(1,0);
+    state.S1 -= 1; state.I1 += 1; birth(host2,host1);
     break;
   case 3:
-    state.S2 -= 1; state.I2 += 1; birth(0,1);
+    state.S2 -= 1; state.I2 += 1; birth(host1,host2);
     break;
   case 4:
-    state.I1 -= 1; state.R1 += 1; death(0);
+    state.I1 -= 1; state.R1 += 1; death(host1);
     break;
   case 5:
-    state.I2 -= 1; state.R2 += 1; death(1);
+    state.I2 -= 1; state.R2 += 1; death(host2);
     break;
   case 6:
     state.R1 -= 1; state.S1 += 1;
@@ -187,16 +191,16 @@ void twospecies_genealogy_t::jump (int event) {
     state.R2 -= 1; state.S2 += 1;
     break;
   case 8:
-    sample(0);
+    sample(host1);
     break;
   case 9:
-    sample(1);
+    sample(host2);
     break;
   case 10:
-    state.S1 -= 1; state.I1 += 1; graft(2); migrate(2,0);
+    state.S1 -= 1; state.I1 += 1; graft(outside); migrate(outside,host1);
     break;
   case 11:
-    state.S2 -= 1; state.I2 += 1; graft(2); migrate(2,1);
+    state.S2 -= 1; state.I2 += 1; graft(outside); migrate(outside,host2);
     break;
   case 12:
     state.S1 -= 1; state.N1 -= 1;
@@ -205,10 +209,10 @@ void twospecies_genealogy_t::jump (int event) {
     state.S2 -= 1; state.N2 -= 1;
     break;
   case 14:
-    state.I1 -= 1; state.N1 -= 1; death(0);
+    state.I1 -= 1; state.N1 -= 1; death(host1);
     break;
   case 15:
-    state.I2 -= 1; state.N2 -= 1; death(1);
+    state.I2 -= 1; state.N2 -= 1; death(host1);
     break;
   case 16:
     state.R1 -= 1; state.N1 -= 1;