Skip to content

Commit

Permalink
working with pomp v5.9.1
Browse files Browse the repository at this point in the history
  • Loading branch information
@kingaa
kingaa committed Jun 9, 2024
1 parent 0e75428 commit 3b9d87e
Show file tree
Hide file tree
Showing 7 changed files with 76 additions and 79 deletions.
6 changes: 3 additions & 3 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,15 +1,15 @@
Package: phylopomp
Type: Package
Title: Phylodynamic Inference for POMP Models
Version: 0.14.0.0
Date: 2024-06-06
Version: 0.14.1.0
Date: 2024-06-09
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"))
)
Description: Tools for phylodynamic analysis.
biocViews:
Depends: R(>= 4.2)
Imports: grid, ggtree, yaml, ggplot2, scales, ape, dplyr, tibble, foreach, tidyr, pomp(>= 5.8.4.0), cowplot
Imports: grid, ggtree, yaml, ggplot2, scales, ape, dplyr, tibble, foreach, tidyr, pomp(>= 5.9.1.0), cowplot
LinkingTo: pomp
Suggests: tidyverse, doFuture, broom
License: GPL-3
Expand Down
2 changes: 1 addition & 1 deletion R/seir.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,7 @@ seirs_pomp <- function (
pomp(
data=NULL,
t0=gi$nodetime[1L],
times=gi$nodetime,
times=gi$nodetime[-1L],
params=c(
Beta=Beta,sigma=sigma,gamma=gamma,psi=psi,omega=omega,
ic,N=sum(ic)
Expand Down
57 changes: 21 additions & 36 deletions src/seirs_pomp.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,16 +114,6 @@ void seirs_rinit
const int *__covindex,
const double *__covars
){
double *color = &COLOR;
const int *nodetype = get_userdata_int("nodetype");
const int *lineage = get_userdata_int("lineage");
const int *index = get_userdata_int("index");
const int *child = get_userdata_int("child");
int nnode = *get_userdata_int("nnode");
#ifndef NDEBUG
const int *saturation = get_userdata_int("saturation");
#endif

double adj = N/(S0+E0+I0+R0);
S = nearbyint(S0*adj);
E = nearbyint(E0*adj);
Expand All @@ -132,30 +122,7 @@ void seirs_rinit
ellE = 0;
ellI = 0;
ll = 0;

int parent = 0;
// for all roots:
while (parent < nnode && nodetype[parent] == 0) {
// color lineages by sampling without replacement
assert(saturation[parent]==1);
int c = child[index[parent]];
assert(lineage[parent]==lineage[c]);
double x = (E-ellE)/(E-ellE + I-ellI);
if (unif_rand() < x) {
// lineage is put into E deme
color[lineage[c]] = 0;
ellE += 1;
ll -= log(x);
} else {
color[lineage[c]] = 1; // lineage is put into I deme
ellI += 1;
ll -= log(1-x);
}
parent++;
}
ellE = nearbyint(ellE);
ellI = nearbyint(ellI);
node = parent-1;
node = 0;
}

//! Simulator for the latent-state process (rprocess).
Expand Down Expand Up @@ -192,14 +159,31 @@ void seirs_gill

int parlin = lineage[parent];
assert(parlin >= 0 && parlin < nsample);
assert(!ISNA(color[parlin]));

// singular portion of filter equation
switch (nodetype[parent]) {
case 0: default:
default: // non-genealogical event
break;
case 0: // root
ll = 0;
// color lineages by sampling without replacement
assert(saturation[parent]==1);
int c = child[index[parent]];
assert(lineage[parent]==lineage[c]);
double x = (E-ellE)/(E-ellE + I-ellI);
if (unif_rand() < x) { // lineage is put into E deme
color[lineage[c]] = 0;
ellE += 1;
ll -= log(x);
} else { // lineage is put into I deme
color[lineage[c]] = 1;
ellI += 1;
ll -= log(1-x);
}
break;
case 1: // sample
ll = 0;
assert(!ISNA(color[parlin]));
// If parent is not in deme I, likelihood = 0.
if (nearbyint(color[parlin]) != 1) {
ll += R_NegInf;
Expand All @@ -219,6 +203,7 @@ void seirs_gill
break;
case 2: // branch point s=(1,1)
ll = 0;
assert(!ISNA(color[parlin]));
// If parent is not in deme I, likelihood = 0.
if (nearbyint(color[parlin]) != 1) {
ll += R_NegInf;
Expand Down
Binary file modified tests/seir2-04.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file modified tests/seir2-05.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
2 changes: 2 additions & 0 deletions tests/seir2.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,9 @@ G |>

po |> rinit(nsim=5)

po |> pfilter(Np=1) |> cond_logLik()
po |> pfilter(Np=1000) |> replicate(n=20) |> concat() -> pf
pf[[1]] |> cond_logLik()
pf |> logLik()
pf |> logLik() |> logmeanexp(se=TRUE,ess=TRUE)

Expand Down
88 changes: 49 additions & 39 deletions tests/seir2.Rout.save
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,6 @@ R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Expand Down Expand Up @@ -59,47 +57,59 @@ Error : in 'seirs_pomp': 'S0','E0','I0','R0' must be nonnegative integers.
>
> po |> rinit(nsim=5)
.id
name [,1] [,2] [,3] [,4] [,5]
S 100.00 100.00 100.00 100.00 100.00
E 3.00 3.00 3.00 3.00 3.00
I 5.00 5.00 5.00 5.00 5.00
R 100.00 100.00 100.00 100.00 100.00
ll 4.03 3.33 3.33 3.33 3.33
node 5.00 5.00 5.00 5.00 5.00
ellE 1.00 2.00 2.00 2.00 2.00
ellI 5.00 4.00 4.00 4.00 4.00
color 1.00 1.00 0.00 0.00 0.00
1.00 0.00 1.00 1.00 1.00
1.00 1.00 0.00 1.00 1.00
NA NA NA NA NA
0.00 0.00 1.00 1.00 1.00
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
1.00 1.00 1.00 0.00 0.00
1.00 1.00 1.00 1.00 1.00
name [,1] [,2] [,3] [,4] [,5]
S 100 100 100 100 100
E 3 3 3 3 3
I 5 5 5 5 5
R 100 100 100 100 100
ll 0 0 0 0 0
node 0 0 0 0 0
ellE 0 0 0 0 0
ellI 0 0 0 0 0
color NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
NA NA NA NA NA
>
> po |> pfilter(Np=1) |> cond_logLik()
[1] 0.4700 0.8473 0.4055 0.9163 1.3863 -0.1717 -1.5208 -Inf -1.1688
[10] -Inf -2.7965 0.0362 -2.1279 -2.8974 -0.5341 -2.4754 -0.9432 -0.4969
[19] -0.9961 -2.5295 -Inf -3.4129 -1.7504 -0.8427 -Inf -0.4058 -2.8705
[28] -4.3717 -Inf -Inf -Inf -0.2195 -3.7279 2.0889 -Inf -Inf
[37] -Inf
> po |> pfilter(Np=1000) |> replicate(n=20) |> concat() -> pf
> pf[[1]] |> cond_logLik()
[1] 0.6883 0.6801 0.6848 0.6191 0.5524 0.1313 -0.8176 -1.2669 -2.2199
[10] 0.0114 -2.3487 -0.7804 -2.5436 -2.1143 -0.7066 -3.2491 -1.8492 -0.2606
[19] -1.2850 -1.4912 0.3858 -2.0186 -1.6742 -0.2236 0.9968 -0.7108 -2.0855
[28] -3.1059 -0.9060 -3.9732 -0.1787 -1.0421 -0.5869 0.9857 1.2883 0.9219
[37] 1.2765
> pf |> logLik()
[1] -5.53 -6.20 -7.95 -7.13 -5.71 -6.44 -6.91 -4.81 -5.33 -5.96 -6.78 -6.47
[13] -5.47 -5.11 -7.86 -4.03 -5.16 -4.70 -5.45 -6.15
[1] -28.2 -28.4 -27.2 -27.3 -29.2 -28.0 -29.4 -29.1 -28.4 -28.5 -30.2 -29.8
[13] -27.5 -27.9 -28.4 -28.3 -29.0 -27.8 -27.9 -28.4
> pf |> logLik() |> logmeanexp(se=TRUE,ess=TRUE)
est se ess
-5.510 0.238 10.207
est se ess
-28.168 0.163 13.475
>
> pf |> plot(type="s")
>
Expand Down

0 comments on commit 3b9d87e

Please sign in to comment.