phylopomp
Phylodynamics for POMPs
Loading...
Searching...
No Matches
curtail.cc
Go to the documentation of this file.
1// Bare genealogy
2
3#include "genealogy.h"
4#include "generics.h"
5#include "internal.h"
6
9void
11(slate_t tnew, slate_t troot)
12{
13 if (tnew < troot) troot = tnew;
14 if (!empty() && tnew < time()) {
15 node_t *p = back();
16 while (!empty() && p->slate > tnew) {
17 ball_t *b;
18 while (p->size() > 1) {
19 b = p->last_ball();
20 switch (b->color) {
21 case black:
22 p->erase(b); delete b;
23 break;
24 case green: case blue: // #nocov
25 assert(0); // #nocov
26 break; // #nocov
27 }
28 }
29 b = p->last_ball();
30 switch (b->color) {
31 case blue:
32 b->color = black;
33 case black:
34 b->deme() = p->deme();
35 swap(b,p->green_ball());
36 case green:
37 destroy_node(p);
38 break;
39 }
40 if (!empty()) p = back();
41 }
42 }
43 time() = tnew;
44 if (!empty() && troot > timezero()) {
45 node_t *p = front();
46 node_t *q;
47 while (!empty() && p->slate < troot) {
48 ball_t *b;
49 assert(p->holds_own());
50 while (p->size() > 1) {
51 b = p->last_ball();
52 switch (b->color) {
53 case blue:
54 p->erase(b); delete b;
55 break;
56 case black:
57 q = make_node(b->deme());
58 q->slate = troot;
59 move(b,p,q); push_back(q);
60 break;
61 case green:
62 q = b->child();
63 if (q == p) {
64 b = p->first_ball();
65 q = b->child();
66 }
67 if (q->slate < troot) {
68 move(b,p,q);
69 } else {
70 node_t *pp = make_node(b->deme());
71 pp->slate = troot;
72 move(b,p,pp); push_back(pp);
73 }
74 break;
75 }
76 }
77 destroy_node(p);
78 if (!empty()) p = front();
79 }
80 sort();
81 }
82 if (troot > timezero()) timezero() = troot;
83}
84
85extern "C" {
86
88 SEXP curtail (SEXP State, SEXP Time, SEXP Troot) {
89 genealogy_t A = State;
90 double t, t0;
91 t = *REAL(AS_NUMERIC(Time));
92 t0 = *REAL(AS_NUMERIC(Troot));
93 if (ISNA(t)) t = A.time();
94 if (ISNA(t0)) t0 = A.timezero();
95 A.curtail(t,t0);
96 SEXP out;
97 PROTECT(out = serial(A));
98 SET_ATTR(out,install("class"),mkString("gpgen"));
99 UNPROTECT(1);
100 return out;
101 }
102
103}
@ green
Definition ball.h:12
@ black
Definition ball.h:12
@ blue
Definition ball.h:12
Balls function as pointers.
Definition ball.h:27
name_t deme(void) const
view deme
Definition ball.h:84
color_t color
Definition ball.h:36
node_t * child(void) const
a child is the owner of a green ball
Definition ball.h:102
Encodes a genealogy.
Definition genealogy.h:19
slate_t & timezero(void)
view/set zero time.
Definition genealogy.h:153
slate_t & time(void)
view/set current time.
Definition genealogy.h:145
node_t * make_node(name_t d=undeme)
Definition genealogy.h:229
void curtail(slate_t tnew, slate_t troot)
Definition curtail.cc:11
Encodes a genealogical node.
Definition node.h:21
name_t deme(void) const
view deme
Definition node.h:96
ball_t * green_ball(void) const
pointer to my green ball
Definition node.h:88
slate_t slate
Definition node.h:33
bool holds_own(void) const
Definition node.h:118
void sort(void)
order nodes in order of increasing time
Definition nodeseq.h:103
void destroy_node(node_t *p)
remove a dead root node
Definition nodeseq.h:208
void swap(ball_t *a, ball_t *b)
swap balls a and b, wherever they lie
Definition nodeseq.h:164
void move(ball_t *b, node_t *p, node_t *q)
move ball b from p to q
Definition nodeseq.h:159
ball_t * last_ball(void) const
retrieve the last ball
Definition pocket.h:126
ball_t * first_ball(void) const
retrieve the first ball
Definition pocket.h:122
SEXP curtail(SEXP State, SEXP Time, SEXP Troot)
curtail the given genealogy
Definition curtail.cc:88
SEXP serial(const TYPE &X)
binary serialization
Definition generics.h:33
double slate_t
Definition internal.h:53