fix bug introduced in previous commit

kingaa · Oct 31, 2024 · ff1c2e8 · ff1c2e8
1 parent 495b854
commit ff1c2e8
Show file tree

Hide file tree

Showing 7 changed files with 18 additions and 24 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: phylopomp
 Type: Package
 Title: Phylodynamic Inference for POMP Models
-Version: 0.14.6.4
-Date: 2024-10-30
+Version: 0.14.6.5
+Date: 2024-10-31
 Authors@R: c(person(given=c("Aaron","A."),family="King",role=c("aut","cre"),email="[email protected]",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/seirs_pomp.c b/src/seirs_pomp.c
@@ -260,9 +260,7 @@ void seirs_gill
 
   // continuous portion of filter equation:
   // take Gillespie steps to the end of the interval
-  // if the state is already incompatible, there is no need for
-  // this, so the state is "frozen".
-  if (tmax > t && R_FINITE(ll)) {
+  if (tmax > t) {
 
     double rate[nrate], logpi[nrate];
     int event;

diff --git a/src/twospecies_pomp.c b/src/twospecies_pomp.c
@@ -417,7 +417,6 @@ void twospecies_gill
       if (unif_rand() < x) {    // s = (2,0)
         color[lineage[c1]] = host1;
         color[lineage[c2]] = host1;
-        ll -= log(x);
         if (S1 > 0) {
           S1 -= 1; I1 += 1; ell1 += 1;
           ll += log(lambda)-log(I1*(I1-1)/2);
@@ -435,7 +434,7 @@ void twospecies_gill
           color[lineage[c1]] = host2;
           color[lineage[c2]] = host1;
         }
-        ll -= log(0.5*(1-x));
+        ll -= log(0.5);
         if (S2 > 0) {
           S2 -= 1; I2 += 1; ell2 += 1;
           ll += log(lambda)-log(I1*I2);
@@ -461,7 +460,6 @@ void twospecies_gill
       if (unif_rand() < x) { // s = (0,2)
         color[lineage[c1]] = host2;
         color[lineage[c2]] = host2;
-        ll -= log(x);
         if (S2 > 0) {
           S2 -= 1; I2 += 1; ell2 += 1;
           ll += log(lambda)-log(I2*(I2-1)/2);
@@ -479,7 +477,7 @@ void twospecies_gill
           color[lineage[c1]] = host2;
           color[lineage[c2]] = host1;
         }
-        ll -= log(0.5*(1-x));
+        ll -= log(0.5);
         if (S1 > 0) {
           S1 -= 1; I1 += 1; ell1 += 1;
           ll += log(lambda)-log(I1*I2);
@@ -502,9 +500,7 @@ void twospecies_gill
 
   // continuous portion of filter equation:
   // take Gillespie steps to the end of the interval.
-  // if the state is already incompatible, there is no need for
-  // this, so the state is "frozen".
-  if (tmax > t && R_FINITE(ll)) {
+  if (tmax > t) {
 
     double rate[nrate], logpi[nrate];
     int event;

diff --git a/tests/seir2-04.png b/tests/seir2-04.png
diff --git a/tests/seir2-05.png b/tests/seir2-05.png
diff --git a/tests/seir2.Rout.save b/tests/seir2.Rout.save
@@ -95,22 +95,22 @@ name    [,1] [,2] [,3] [,4] [,5]
 > 
 > po |> pfilter(Np=1) |> cond_logLik()
  [1]  0.470  0.560  0.693  0.916  0.288  1.631 -0.269 -2.670 -1.538  0.389
-[11] -1.417 -1.069 -0.820 -2.058   -Inf   -Inf -0.905 -0.326 -1.151 -2.791
-[21]  1.535 -3.761   -Inf   -Inf  2.485   -Inf -3.491 -2.972   -Inf -5.262
-[31]   -Inf -0.253 -2.967   -Inf   -Inf   -Inf   -Inf
+[11] -1.417 -1.069 -0.820 -2.058   -Inf   -Inf -1.237 -0.158   -Inf -1.225
+[21]  1.047   -Inf   -Inf -0.841   -Inf   -Inf -5.046   -Inf   -Inf   -Inf
+[31]   -Inf -0.348 -5.306   -Inf   -Inf  1.135   -Inf
 > po |> pfilter(Np=1000) |> replicate(n=20) |> concat() -> pf
 > pf[[1]] |> cond_logLik()
- [1]  0.6846  0.6834  0.6959  0.6396  0.5158  0.1043 -0.7920 -1.1724 -2.1462
-[10] -0.1268 -2.5962 -0.8908 -2.6324 -1.9143 -0.6320 -3.3068 -1.5696 -0.2448
-[19] -1.4413 -1.7519 -0.0458 -1.9268 -1.5491 -0.3602  1.0646 -0.6862 -1.9760
-[28] -2.8387 -0.9644 -5.1106  0.0977 -0.9132 -1.0709  1.1780  1.4694  1.0368
-[37]  1.5009
+ [1]  0.6808  0.6957  0.6683  0.6321  0.5753  0.0884 -0.7797 -1.3926 -2.0394
+[10] -0.0145 -2.1752 -0.7526 -2.4744 -2.0605 -0.6792 -2.7662 -1.7165 -0.2187
+[19] -1.4210 -1.6484  0.1526 -2.2243 -1.7258 -0.2450  1.0954 -0.6398 -2.0069
+[28] -2.5899 -1.2573 -4.1875 -0.1862 -0.8248 -0.6212  1.1347  0.9960  0.8011
+[37]  1.1727
 > pf |> logLik()
- [1] -29.0 -27.6 -27.9 -27.6 -28.3 -28.0 -28.1 -28.0 -29.5 -27.5 -30.7 -28.7
-[13] -30.4 -28.6 -28.7 -28.9 -27.5 -29.4 -28.6 -27.3
+ [1] -28.0 -28.8 -27.4 -26.4 -28.7 -28.7 -28.2 -28.2 -29.2 -27.8 -28.7 -27.3
+[13] -27.3 -28.4 -27.9 -29.8 -29.1 -29.5 -29.3 -28.0
 > pf |> logLik() |> logmeanexp(se=TRUE,ess=TRUE)
     est      se     ess 
--28.199   0.158  13.679 
+-27.969   0.248   9.955 
 > 
 > pf |> plot(type="s")
 > 

diff --git a/tests/twospecies1.Rout.save b/tests/twospecies1.Rout.save
@@ -117,7 +117,7 @@ The following object is masked from 'package:stats':
 +   logLik() |>
 +   logmeanexp(se=TRUE,ess=TRUE)
          est           se          ess 
--119.8921708    0.9266098    1.8998700 
+-203.9170186    0.3959315    3.3136507 
 > 
 > try(
 +   s |>