curtail now works from left as well as right

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: phylopomp
 Type: Package
 Title: Phylodynamic Inference for POMP Models
-Version: 0.14.8.2
-Date: 2024-11-25
+Version: 0.14.9.0
+Date: 2024-12-02
 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"))
 	  )

R/curtail.R

@@ -5,10 +5,11 @@
 ##' @name curtail
 ##' @include getinfo.R
 ##' @inheritParams getInfo
+##' @param troot new root time for curtailed genealogy
 ##' @return A curtailed genealogy object.
 ##' @example examples/curtail.R
 ##' @rdname curtail
 ##' @export
-curtail <- function (object, time = NA, prune = TRUE, obscure = TRUE) {
-  .Call(P_curtail,geneal(object),time)
+curtail <- function (object, time = NA, troot = NA) {
+  .Call(P_curtail,geneal(object),time,troot)
 }

R/treeplot.R

@@ -5,7 +5,7 @@
 ##' @name treeplot
 ##' @include getinfo.R diagram.R
 ##' @param tree character; tree representation in Newick format.
-##' @param t0 numeric; time of the root.
+##' @param t0 numeric; root time.
 ##' @param time numeric; time of the genealogy.
 ##' @param ladderize Ladderize?
 ##' @param points Show nodes and tips?

TODO.md

@@ -1,5 +1,6 @@
 ## To-do List
 
+- curtail from the left also
 - LBDP and SIRS into new format
 - use integer time?
 

inst/NEWS

@@ -2,6 +2,9 @@ _N_e_w_s _f_o_r _p_a_c_k_a_g_e '_p_h_y_l_o_p_o_m_p'
 
 _C_h_a_n_g_e_s _i_n '_p_h_y_l_o_p_o_m_p' _v_e_r_s_i_o_n _0._1_4:
 
+        • The ‘curtail()’ function now allows one to curtail a
+          genealogy from the left as well as the right.
+
         • The filtering codes for the two-species model are now
           working.
 

inst/NEWS.Rd

@@ -2,6 +2,7 @@
 \title{News for package `phylopomp'}
 \section{Changes in \pkg{phylopomp} version 0.14}{
   \itemize{
+    \item The \code{curtail()} function now allows one to curtail a genealogy from the left as well as the right.
     \item The filtering codes for the two-species model are now working.
     \item A new sample/death event-type has been introduced into the \code{genealogy_t} type.
     \item Corrections of filtering equations for the SEIR and Two-species models.

src/bare.cc

@@ -7,9 +7,14 @@
 extern "C" {
 
   //! curtail the given genealogy
-  SEXP curtail (SEXP State, SEXP Time) {
+  SEXP curtail (SEXP State, SEXP Time, SEXP Troot) {
     genealogy_t A = State;
-    A.curtail(*REAL(AS_NUMERIC(Time)));
+    double t, t0;
+    t = *REAL(AS_NUMERIC(Time));
+    t0 = *REAL(AS_NUMERIC(Troot));
+    if (ISNA(t)) t = A.time();
+    if (ISNA(t0)) t0 = A.timezero();
+    A.curtail(t,t0);
     SEXP out;
     PROTECT(out = serial(A));
     SET_ATTR(out,install("class"),mkString("gpgen"));

src/genealogy.h

@@ -491,9 +491,10 @@ class genealogy_t : public nodeseq_t {
     return *this;
   };
   //! curtail the genealogy by removing nodes
-  //! with times later than tnew
-  void curtail (slate_t tnew) {
-    if (!empty()) {
+  //! with times later than tnew and/or earlier than troot
+  void curtail (slate_t tnew, slate_t troot) {
+    if (tnew < troot) troot = tnew;
+    if (!empty() && tnew < time()) {
       node_t *p = back();
       while (!empty() && p->slate > tnew) {
         ball_t *b;
@@ -504,7 +505,7 @@ class genealogy_t : public nodeseq_t {
             p->erase(b); delete b;
             break;
           case green: case blue: // #nocov
-            assert(0);       // #nocov
+            assert(0);           // #nocov
             break;               // #nocov
           }
         }
@@ -523,7 +524,46 @@ class genealogy_t : public nodeseq_t {
       }
     }
     if (tnew < time()) _time = tnew;
-    if (tnew < timezero()) _t0 = tnew;
+    if (!empty() && troot > timezero()) {
+      node_t *p = front();
+      node_t *q;
+      while (!empty() && p->slate <= troot) {
+        ball_t *b;
+        assert(p->holds_own());
+        while (p->size() > 1) {
+          b = p->last_ball();
+          switch (b->color) {
+          case blue:
+            p->erase(b); delete b;
+            break;
+          case black:
+            q = make_node(b->deme());
+            q->slate = troot;
+            b->holder() = q;
+            q->insert(b); p->erase(b);
+            push_back(q);
+            break;
+          case green:
+            q = b->child();
+            if (q->slate < troot) {
+              q->insert(b); p->erase(b);
+              b->holder() = q;
+            } else {
+              node_t *pp = make_node(b->deme());
+              pp->slate = troot;
+              pp->insert(b); p->erase(b);
+              b->holder() = pp;
+              push_back(pp);
+            }
+            break;
+          }
+        }
+        destroy_node(p);
+        if (!empty()) p = front();
+      }
+      sort();
+    }
+    if (troot > timezero()) _t0 = troot;
   };
   //! merge two genealogies:
   //! 1. the node-sequences are merged;

src/init.c

@@ -8,7 +8,7 @@ get_userdata_int_t *get_userdata_int;
 
 SEXP parse_newick (SEXP, SEXP, SEXP);
 SEXP getInfo (SEXP);
-SEXP curtail (SEXP, SEXP);
+SEXP curtail (SEXP, SEXP, SEXP);
 SEXP yaml (SEXP);
 SEXP gendat (SEXP);
 
@@ -34,7 +34,7 @@ static const R_CallMethodDef callMethods[] = {
   METHODS(SIR),
   METHODS(TwoSpecies),
   {"parse_newick", (DL_FUNC) &parse_newick, 3},
-  {"curtail", (DL_FUNC) &curtail, 2},
+  {"curtail", (DL_FUNC) &curtail, 3},
   {"yaml", (DL_FUNC) &yaml, 1},
   {"gendat", (DL_FUNC) &gendat, 1},
   {NULL, NULL, 0}

src/parse.cc

@@ -19,7 +19,7 @@ extern "C" {
     G.parse(x,t0);
     if (!ISNA(tf)) {
       if (G.time() > tf) {
-        G.curtail(tf);
+        G.curtail(tf,t0);
       } else {
         G.time() = tf;
       }

tests/curtail.R

@@ -1,35 +1,60 @@
 png(filename="curtail-%02d.png",res=100,width=6,height=4,units="in")
 
+options(tidyverse.quiet=TRUE,digits=3)
 suppressPackageStartupMessages({
   library(tidyverse)
   library(phylopomp)
 })
 theme_set(theme_bw())
 set.seed(24963811)
-options(digits=3)
 
 runSEIR(time=5,S0=20,omega=2) -> x
 
 plot_grid(
   x |>
     plot(points=TRUE,prune=FALSE,obscure=FALSE,ladderize=FALSE)+
-    geom_vline(xintercept=3.5),
-  x |> lineages(obscure=FALSE,prune=FALSE) |>
+    geom_vline(xintercept=c(1.5,3.5)),
+  x |>
+    lineages(obscure=FALSE,prune=FALSE) |>
     plot()+
     guides(color="none")+
-    geom_vline(xintercept=3.5)+
-    expand_limits(x=5,y=10),
-  x |> curtail(time=3.5) |>
+    geom_vline(xintercept=c(1.5,3.5))+
+    expand_limits(x=c(0,5),y=10),
+  x |>
+    curtail(time=3.5) |>
     plot(points=TRUE,prune=FALSE,obscure=FALSE,ladderize=FALSE)+
     geom_vline(xintercept=3.5)+
-    expand_limits(x=5),
-  x |> curtail(time=3.5) |>
+    expand_limits(x=c(0,5)),
+  x |>
+    curtail(time=3.5) |>
     lineages(obscure=FALSE,prune=FALSE) |>
     plot()+geom_vline(xintercept=3.5)+
     guides(color="none")+
-    expand_limits(x=5,y=10),
+    expand_limits(x=c(0,5),y=10),
+  x |>
+    curtail(troot=1.5) |>
+    plot(points=TRUE,prune=FALSE,obscure=FALSE,ladderize=FALSE)+
+    geom_vline(xintercept=1.5)+
+    expand_limits(x=c(0,5)),
+  x |>
+    curtail(troot=1.5) |>
+    lineages(obscure=FALSE,prune=FALSE) |>
+    plot()+geom_vline(xintercept=1.5)+
+    guides(color="none")+
+    expand_limits(x=c(0,5),y=10),
+  x |>
+    curtail(troot=1.5,time=3.5) |>
+    plot(points=TRUE,prune=FALSE,obscure=FALSE,ladderize=FALSE)+
+    geom_vline(xintercept=1.5)+
+    expand_limits(x=c(0,5)),
+  x |>
+    curtail(troot=1.5,time=3.5) |>
+    lineages(obscure=FALSE,prune=FALSE) |>
+    plot()+geom_vline(xintercept=1.5)+
+    guides(color="none")+
+    expand_limits(x=c(0,5),y=10),
   ncol=2,
-  align="hv",axis="tblr"
+  align="v",axis="bl"
 )
 
 x |>
@@ -47,6 +72,21 @@ plot_grid(
   x |>
     curtail(time=0.5) |>
     diagram(prune=FALSE,obscure=FALSE),
+  x |>
+    curtail(troot=4.5) |>
+    diagram(prune=FALSE,obscure=FALSE),
+  x |>
+    curtail(troot=0.5,time=1) |>
+    diagram(prune=FALSE,obscure=FALSE),
+  ncol=1
+)
+
+runSEIR(time=3,S0=5,omega=2) -> x
+plot_grid(
+  x |> curtail() |> diagram(prune=FALSE,obscure=FALSE),
+  x |> curtail(troot=1) |> diagram(prune=FALSE,obscure=FALSE),
+  x |> curtail(time=2) |> diagram(prune=FALSE,obscure=FALSE),
+  x |> curtail(time=2,troot=1) |> diagram(prune=FALSE,obscure=FALSE),
   ncol=1
 )
 

tests/curtail.Rout.save

@@ -19,36 +19,61 @@ Type 'q()' to quit R.
 
 > png(filename="curtail-%02d.png",res=100,width=6,height=4,units="in")
 > 
+> options(tidyverse.quiet=TRUE,digits=3)
 > suppressPackageStartupMessages({
 +   library(tidyverse)
 +   library(phylopomp)
 + })
 > theme_set(theme_bw())
 > set.seed(24963811)
-> options(digits=3)
 > 
 > runSEIR(time=5,S0=20,omega=2) -> x
 > 
 > plot_grid(
 +   x |>
 +     plot(points=TRUE,prune=FALSE,obscure=FALSE,ladderize=FALSE)+
-+     geom_vline(xintercept=3.5),
-+   x |> lineages(obscure=FALSE,prune=FALSE) |>
++     geom_vline(xintercept=c(1.5,3.5)),
++   x |>
++     lineages(obscure=FALSE,prune=FALSE) |>
 +     plot()+
 +     guides(color="none")+
-+     geom_vline(xintercept=3.5)+
-+     expand_limits(x=5,y=10),
-+   x |> curtail(time=3.5) |>
++     geom_vline(xintercept=c(1.5,3.5))+
++     expand_limits(x=c(0,5),y=10),
++   x |>
++     curtail(time=3.5) |>
 +     plot(points=TRUE,prune=FALSE,obscure=FALSE,ladderize=FALSE)+
 +     geom_vline(xintercept=3.5)+
-+     expand_limits(x=5),
-+   x |> curtail(time=3.5) |>
++     expand_limits(x=c(0,5)),
++   x |>
++     curtail(time=3.5) |>
 +     lineages(obscure=FALSE,prune=FALSE) |>
 +     plot()+geom_vline(xintercept=3.5)+
 +     guides(color="none")+
-+     expand_limits(x=5,y=10),
++     expand_limits(x=c(0,5),y=10),
++   x |>
++     curtail(troot=1.5) |>
++     plot(points=TRUE,prune=FALSE,obscure=FALSE,ladderize=FALSE)+
++     geom_vline(xintercept=1.5)+
++     expand_limits(x=c(0,5)),
++   x |>
++     curtail(troot=1.5) |>
++     lineages(obscure=FALSE,prune=FALSE) |>
++     plot()+geom_vline(xintercept=1.5)+
++     guides(color="none")+
++     expand_limits(x=c(0,5),y=10),
++   x |>
++     curtail(troot=1.5,time=3.5) |>
++     plot(points=TRUE,prune=FALSE,obscure=FALSE,ladderize=FALSE)+
++     geom_vline(xintercept=1.5)+
++     expand_limits(x=c(0,5)),
++   x |>
++     curtail(troot=1.5,time=3.5) |>
++     lineages(obscure=FALSE,prune=FALSE) |>
++     plot()+geom_vline(xintercept=1.5)+
++     guides(color="none")+
++     expand_limits(x=c(0,5),y=10),
 +   ncol=2,
-+   align="hv",axis="tblr"
++   align="v",axis="bl"
 + )
 > 
 > x |>
@@ -78,6 +103,21 @@ $newick
 +   x |>
 +     curtail(time=0.5) |>
 +     diagram(prune=FALSE,obscure=FALSE),
++   x |>
++     curtail(troot=4.5) |>
++     diagram(prune=FALSE,obscure=FALSE),
++   x |>
++     curtail(troot=0.5,time=1) |>
++     diagram(prune=FALSE,obscure=FALSE),
++   ncol=1
++ )
+> 
+> runSEIR(time=3,S0=5,omega=2) -> x
+> plot_grid(
++   x |> curtail() |> diagram(prune=FALSE,obscure=FALSE),
++   x |> curtail(troot=1) |> diagram(prune=FALSE,obscure=FALSE),
++   x |> curtail(time=2) |> diagram(prune=FALSE,obscure=FALSE),
++   x |> curtail(time=2,troot=1) |> diagram(prune=FALSE,obscure=FALSE),
 +   ncol=1
 + )
 >