phylopomp
Phylodynamics for POMPs
Loading...
Searching...
No Matches
genealogy_t Class Reference

Encodes a genealogy. More...

#include <genealogy.h>

Inheritance diagram for genealogy_t:
Collaboration diagram for genealogy_t:

Public Member Functions

size_t ndeme (void) const
 number of demes
 
size_t & ndeme (void)
 number of demes
 
size_t bytesize (void) const
 size of serialized binary form
 
 genealogy_t (double t0=R_NaReal, size_t ndeme=0)
 
 genealogy_t (raw_t *o)
 constructor from serialized binary form
 
 genealogy_t (SEXP o)
 constructor from RAW SEXP (containing binary serialization)
 
 genealogy_t (const genealogy_t &G)
 copy constructor
 
genealogy_toperator= (const genealogy_t &G)
 copy assignment operator
 
 genealogy_t (genealogy_t &&)=default
 move constructor
 
genealogy_toperator= (genealogy_t &&)=default
 move assignment operator
 
 ~genealogy_t (void)
 destructor
 
slate_ttime (void)
 view/set current time.
 
slate_t time (void) const
 view current time.
 
slate_ttimezero (void)
 view/set zero time.
 
slate_t timezero (void) const
 get zero time.
 
void lineage_count (double *tout, int *deme, int *ell, int *sat, int *etype) const
 
SEXP lineage_count (void) const
 lineage count and saturation
 
void gendat (double *tout, int *anc, int *lin, int *sat, int *type, int *deme, int *index, int *child) const
 genealogy information in list format
 
SEXP gendat (void) const
 genealogy information in list format
 
size_t nsample (void) const
 number of samples
 
size_t nroot (void) const
 number of roots
 
string_t yaml (string_t tab="") const
 human/machine-readable info
 
SEXP structure (void) const
 R list description.
 
string_t newick (bool extended=true) const
 put genealogy at current time into Newick format.
 
void valid (void) const
 check the validity of the genealogy.
 
bool check_genealogy_size (size_t grace=0) const
 check the size of the genealogy (to prevent memory exhaustion).
 
ball_tbirth (ball_t *a, slate_t t, name_t d)
 birth into deme d
 
ball_tbirth (node_t *p, name_t d)
 birth of second or subsequent sibling into deme d
 
void death (ball_t *a, slate_t t)
 death
 
ball_tgraft (slate_t t, name_t d)
 graft a new lineage into deme d
 
void sample (ball_t *a, slate_t t)
 insert a sample node
 
void sample_death (ball_t *a, slate_t t)
 insert a sample node and simultaneously terminate the lineage
 
void migrate (ball_t *a, slate_t t, name_t d=0)
 movement into deme d
 
void sample_migrate (ball_t *a, slate_t t, name_t d=0)
 insert a sample node and simultaneously migrate the lineage
 
genealogy_tprune (void)
 prune the tree (drop all black balls)
 
genealogy_tobscure (void)
 erase all deme information
 
void curtail (slate_t tnew, slate_t troot)
 
genealogy_toperator+= (const genealogy_t &other)
 merges two genealogies, adjusting time, t0, and ndeme as needed
 
void insert_zlb (void)
 insert zero-length branches for samples where needed
 
void drop_zlb (void)
 drop all zero-length branches
 
genealogy_tparse (const string_t &s)
 Parse a Newick string and create the indicated genealogy.
 
- Public Member Functions inherited from nodeseq_t
 ~nodeseq_t (void)
 destructor
 
size_t bytesize (void) const
 size of serialized binary form
 
void sort (void)
 order nodes in order of increasing time
 
pocket_tcolored (color_t col) const
 Get all balls of a color.
 
size_t ntime (slate_t t) const
 Number of distinct timepoints.
 
size_t length (void) const
 Number of nodes in the sequence.
 
node_tposition (int n)
 traverse to nth node, retrieve pointer
 
void move (ball_t *b, node_t *p, node_t *q)
 move ball b from p to q
 
void swap (ball_t *a, ball_t *b)
 swap balls a and b, wherever they lie
 
void attach (node_t *p, node_t *q)
 
void detach (node_t *p)
 
void add (node_t *p, ball_t *a)
 
void drop (ball_t *a)
 
void destroy_node (node_t *p)
 remove a dead root node
 
void weed (void)
 drop all dead roots
 
void comb (void)
 
void trace_lineages (void)
 
string_t yaml (string_t tab="") const
 human/machine-readable info
 
SEXP structure (void) const
 R list description.
 
string_t newick (slate_t t, bool showdeme, bool extended) const
 put genealogy at time t into Newick format.
 

Private Member Functions

name_t unique (void)
 get the next unique name
 
void clean (void)
 clean up
 
node_tmake_node (name_t d=undeme)
 
void reuniqify (name_t shift)
 shifts names to avoid overlap
 
void repair_tips (void)
 
void repair_roots (void)
 roots are added at zero time if needed
 
node_tscan_branch (string_t::const_iterator b, string_t::const_iterator e)
 

Private Attributes

name_t _unique
 The next unique name.
 
slate_t _t0
 The initial time.
 
slate_t _time
 The current time.
 
size_t _ndeme
 The number of demes (excluding the undeme).
 

Static Private Attributes

static const name_t magic = 1123581321
 

Friends

raw_toperator>> (const genealogy_t &G, raw_t *o)
 binary serialization
 
raw_toperator>> (raw_t *o, genealogy_t &G)
 binary deserialization
 

Additional Inherited Members

- Static Public Member Functions inherited from nodeseq_t
static bool compare (node_t *p, node_t *q)
 

Detailed Description

Encodes a genealogy.

A genealogy consists of a sequence of nodes and the current time.

Definition at line 19 of file genealogy.h.

Constructor & Destructor Documentation

◆ genealogy_t() [1/5]

genealogy_t::genealogy_t ( double t0 = R_NaReal,
size_t ndeme = 0 )
inline

basic constructor for genealogy class

  • t0 = initial time
  • ndeme = number of demes (excluding the undeme)

Definition at line 102 of file genealogy.h.

102 {
103 clean();
104 _time = _t0 = slate_t(t0);
105 _ndeme = ndeme;
106 };
void clean(void)
clean up
Definition genealogy.h:50
size_t ndeme(void) const
number of demes
Definition genealogy.h:59
slate_t _time
The current time.
Definition genealogy.h:34
size_t _ndeme
The number of demes (excluding the undeme).
Definition genealogy.h:36
slate_t _t0
The initial time.
Definition genealogy.h:32
double slate_t
Definition internal.h:53
Here is the call graph for this function:
Here is the caller graph for this function:

◆ genealogy_t() [2/5]

genealogy_t::genealogy_t ( raw_t * o)
inline

constructor from serialized binary form

Definition at line 108 of file genealogy.h.

108 {
109 o >> *this;
110 };

◆ genealogy_t() [3/5]

genealogy_t::genealogy_t ( SEXP o)
inline

constructor from RAW SEXP (containing binary serialization)

Definition at line 112 of file genealogy.h.

112 {
113 if (LENGTH(o)==0)
114 err("in %s: cannot deserialize a NULL.",__func__);
115 PROTECT(o = AS_RAW(o));
116 RAW(o) >> *this;
117 UNPROTECT(1);
118 };
#define err(...)
Definition internal.h:18

◆ genealogy_t() [4/5]

genealogy_t::genealogy_t ( const genealogy_t & G)
inline

copy constructor

Definition at line 120 of file genealogy.h.

120 {
121 raw_t *o = new raw_t[G.bytesize()];
122 G >> o;
123 o >> *this;
124 delete[] o;
125 };
size_t bytesize(void) const
size of serialized binary form
Definition genealogy.h:71
Rbyte raw_t
Definition internal.h:52
Here is the call graph for this function:

◆ genealogy_t() [5/5]

genealogy_t::genealogy_t ( genealogy_t && )
default

move constructor

Here is the call graph for this function:

◆ ~genealogy_t()

genealogy_t::~genealogy_t ( void )
inline

destructor

Definition at line 140 of file genealogy.h.

140 {
141 clean();
142 };
Here is the call graph for this function:

Member Function Documentation

◆ birth() [1/2]

ball_t * genealogy_t::birth ( ball_t * a,
slate_t t,
name_t d )
inline

birth into deme d

Definition at line 242 of file genealogy.h.

242 {
243 time() = t;
244 node_t *p = make_node(a->deme());
245 ball_t *b = new ball_t (p,p->uniq,black,d);
246 p->insert(b);
247 p->slate = time();
248 add(p,a);
249 return b;
250 };
@ black
Definition ball.h:12
name_t deme(void) const
view deme
Definition ball.h:84
slate_t & time(void)
view/set current time.
Definition genealogy.h:145
node_t * make_node(name_t d=undeme)
Definition genealogy.h:229
name_t uniq
Definition node.h:32
void insert(ball_t *a)
insert a ball into the pocket of a node
Definition node.h:156
slate_t slate
Definition node.h:33
void add(node_t *p, ball_t *a)
Definition nodeseq.h:184
Here is the call graph for this function:

◆ birth() [2/2]

ball_t * genealogy_t::birth ( node_t * p,
name_t d )
inline

birth of second or subsequent sibling into deme d

Definition at line 252 of file genealogy.h.

252 {
253 ball_t *b = new ball_t(p,unique(),black,d);
254 p->insert(b);
255 return b;
256 };
name_t unique(void)
get the next unique name
Definition genealogy.h:43
Here is the call graph for this function:

◆ bytesize()

size_t genealogy_t::bytesize ( void ) const
inline

size of serialized binary form

Definition at line 71 of file genealogy.h.

71 {
72 return 3*sizeof(name_t) +
73 2*sizeof(slate_t) + nodeseq_t::bytesize();
74 };
size_t bytesize(void) const
size of serialized binary form
Definition nodeseq.h:38
size_t name_t
Definition internal.h:54
Here is the call graph for this function:
Here is the caller graph for this function:

◆ check_genealogy_size()

bool genealogy_t::check_genealogy_size ( size_t grace = 0) const
inline

check the size of the genealogy (to prevent memory exhaustion).

Definition at line 214 of file genealogy.h.

214 {
215 static size_t maxq = MEMORY_MAX/(sizeof(node_t)+2*sizeof(ball_t));
216 bool ok = true;
217 if (size() > maxq+grace) {
218 err("maximum genealogy size exceeded!"); // #nocov
219 } else if (size() > maxq) {
220 ok = false; // #nocov
221 }
222 return ok;
223 };
static const size_t MEMORY_MAX
Definition genealogy.h:13
Here is the caller graph for this function:

◆ clean()

void genealogy_t::clean ( void )
inlineprivate

clean up

Definition at line 50 of file genealogy.h.

50 {
51 _unique = 0;
52 _ndeme = 0;
53 _t0 = _time = R_NaReal;
54 };
name_t _unique
The next unique name.
Definition genealogy.h:30
Here is the caller graph for this function:

◆ curtail()

void genealogy_t::curtail ( slate_t tnew,
slate_t troot )

curtail the genealogy by removing nodes with times later than tnew and/or earlier than troot

Definition at line 10 of file curtail.cc.

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}
@ green
Definition ball.h:12
@ blue
Definition ball.h:12
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
slate_t & timezero(void)
view/set zero time.
Definition genealogy.h:153
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
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
Here is the call graph for this function:
Here is the caller graph for this function:

◆ death()

void genealogy_t::death ( ball_t * a,
slate_t t )
inline

death

Definition at line 258 of file genealogy.h.

258 {
259 time() = t;
260 drop(a);
261 };
void drop(ball_t *a)
Definition nodeseq.h:192
Here is the call graph for this function:

◆ drop_zlb()

void genealogy_t::drop_zlb ( void )
inline

drop all zero-length branches

Definition at line 373 of file genealogy.h.

373 {
374 for (node_t *p : *this) {
375 if (!p->holds_own() &&
376 p->slate == p->parent()->slate &&
377 p->deme() == p->parent()->deme()) {
378 while (!p->empty()) {
379 ball_t *b = p->last_ball();
380 p->erase(b); p->parent()->insert(b);
381 }
382 detach(p);
383 }
384 }
385 };
node_t * parent(void) const
Definition node.h:115
void detach(node_t *p)
Definition nodeseq.h:179
Here is the call graph for this function:
Here is the caller graph for this function:

◆ gendat() [1/2]

void genealogy_t::gendat ( double * tout,
int * anc,
int * lin,
int * sat,
int * type,
int * deme,
int * index,
int * child ) const

genealogy information in list format

Definition at line 8 of file gendat.cc.

11 {
12 int m, n, k;
13 node_it i, j;
14 for (k = 0, n = 0, i = begin(); i != end(); i++, n++) {
15 node_t *p = *i;
16 assert(!p->holds(black)); // tree should be pruned first
17 tout[n] = p->slate;
18 deme[n] = p->deme();
19 if (p->is_root()) {
20 type[n] = 0; // root node
21 } else if (p->holds(blue)) {
22 type[n] = 1; // sample node
23 deme[n] = p->ball(blue)->deme();
24 } else {
25 type[n] = 2; // internal node
26 }
27 lin[n] = p->lineage(); // 0-based indexing
28 sat[n] = p->nchildren();
29 index[n] = k;
30 k += sat[n];
31 child[n] = NA_INTEGER;
32 if (p->is_root()) {
33 anc[n] = n; // 0-based indexing
34 } else {
35 for (m = 0, j = begin(); j != i; j++, m++) {
36 node_t *q = *j;
37 if (p->parent()->uniq == q->uniq) {
38 anc[n] = m;
39 break;
40 }
41 }
42 }
43 }
44 tout[n] = time();
45 for (k = 0, n = 0, i = begin(); i != end(); i++, n++) {
46 node_t *p = *i;
47 j = i; j++;
48 for (m = n+1; j != end(); m++, j++) {
49 node_t *q = *j;
50 if (p->uniq == q->parent()->uniq) {
51 child[k++] = m;
52 }
53 }
54 }
55}
name_t lineage(void) const
view lineage
Definition node.h:104
bool is_root(void) const
Definition node.h:121
int nchildren(void) const
number of descendants
Definition node.h:131
ball_t * ball(const color_t c) const
retrieve the first ball of the specified color.
Definition pocket.h:130
bool holds(ball_t *b) const
does this node hold the given ball?
Definition pocket.h:109
static const int deme
Definition lbdp.cc:7
#define n
Definition lbdp_pomp.c:9
std::list< node_t * >::const_iterator node_it
Definition nodeseq.h:12
Here is the call graph for this function:
Here is the caller graph for this function:

◆ gendat() [2/2]

SEXP genealogy_t::gendat ( void ) const

genealogy information in list format

Definition at line 59 of file gendat.cc.

60 {
61 SEXP t0, tout, anc, lin, sat, type, deme, index, child, ns, nr, nn;
62 SEXP out, outn;
63 size_t n = length();
64 PROTECT(t0 = NEW_NUMERIC(1));
65 PROTECT(tout = NEW_NUMERIC(n+1));
66 PROTECT(type = NEW_INTEGER(n));
67 PROTECT(deme = NEW_INTEGER(n));
68 PROTECT(lin = NEW_INTEGER(n));
69 PROTECT(sat = NEW_INTEGER(n));
70 PROTECT(index = NEW_INTEGER(n));
71 PROTECT(child = NEW_INTEGER(n));
72 PROTECT(anc = NEW_INTEGER(n));
73 PROTECT(ns = NEW_INTEGER(1));
74 PROTECT(nr = NEW_INTEGER(1));
75 PROTECT(nn = NEW_INTEGER(1));
76 PROTECT(out = NEW_LIST(12));
77 PROTECT(outn = NEW_CHARACTER(12));
78 set_list_elem(out,outn,t0,"t0",0);
79 set_list_elem(out,outn,tout,"nodetime",1);
80 set_list_elem(out,outn,type,"nodetype",2);
81 set_list_elem(out,outn,deme,"deme",3);
82 set_list_elem(out,outn,lin,"lineage",4);
83 set_list_elem(out,outn,sat,"saturation",5);
84 set_list_elem(out,outn,index,"index",6);
85 set_list_elem(out,outn,child,"child",7);
86 set_list_elem(out,outn,anc,"ancestor",8);
87 set_list_elem(out,outn,ns,"nsample",9);
88 set_list_elem(out,outn,nr,"nroot",10);
89 set_list_elem(out,outn,nn,"nnode",11);
90 SET_NAMES(out,outn);
91 gendat(REAL(tout),INTEGER(anc),INTEGER(lin),INTEGER(sat),
92 INTEGER(type),INTEGER(deme),INTEGER(index),INTEGER(child));
93 *REAL(t0) = double(timezero()); // zero-time
94 *INTEGER(ns) = nsample(); // number of samples
95 *INTEGER(nr) = nroot(); // number of roots
96 *INTEGER(nn) = length(); // number of nodes
97 UNPROTECT(14);
98 return out;
99}
size_t nsample(void) const
number of samples
Definition genealogy.h:183
SEXP gendat(void) const
genealogy information in list format
Definition gendat.cc:60
size_t nroot(void) const
number of roots
Definition genealogy.h:192
size_t length(void) const
Number of nodes in the sequence.
Definition nodeseq.h:142
static int set_list_elem(SEXP list, SEXP names, SEXP element, const char *name, int pos)
Definition internal.h:76
Here is the call graph for this function:
Here is the caller graph for this function:

◆ graft()

ball_t * genealogy_t::graft ( slate_t t,
name_t d )
inline

graft a new lineage into deme d

Definition at line 263 of file genealogy.h.

263 {
264 time() = t;
265 node_t *p = make_node();
266 ball_t *b = new ball_t (p,p->uniq,black,d);
267 p->insert(b);
268 p->slate = timezero();
269 push_front(p);
270 return b;
271 };
Here is the call graph for this function:

◆ insert_zlb()

void genealogy_t::insert_zlb ( void )
inline

insert zero-length branches for samples where needed

Definition at line 358 of file genealogy.h.

358 {
359 for (node_t *p : *this) {
360 if (p->holds(green) && p->holds(blue)) {
361 assert(!p->holds(black)); // genealogy should have already been pruned
362 ball_t *b = p->last_ball();
363 assert(b->is(blue));
364 node_t *q = make_node(p->deme());
365 q->slate = p->slate;
366 swap(q->green_ball(),b);
367 push_back(q);
368 }
369 }
370 sort();
371 };
bool is(color_t c) const
is a given ball of the given color?
Definition ball.h:115
Here is the call graph for this function:
Here is the caller graph for this function:

◆ lineage_count() [1/2]

void genealogy_t::lineage_count ( double * tout,
int * deme,
int * ell,
int * sat,
int * etype ) const

lineage count, saturation, and event-type. types are:

  • 0 = non-event
  • -1 = root
  • 1 = sample
  • 2 = non-sample node
  • 3 = end of interval

Definition at line 39 of file lineages.cc.

42{
43 size_t nd = ndeme()+1;
44 slate_t tcur = timezero();
45 for (size_t j = 0; j < nd; j++) {
46 tout[j] = tcur;
47 deme[j] = j;
48 sat[j] = ell[j] = 0;
49 etype[j] = 0;
50 }
51 for (const node_t *p : *this) {
52 if (tcur < p->slate) {
53 tout += nd; ell += nd; sat += nd;
54 deme += nd; etype += nd;
55 tcur = p->slate;
56 for (size_t j = 0; j < nd; j++) {
57 tout[j] = tcur;
58 deme[j] = j;
59 ell[j] = (ell-nd)[j];
60 sat[j] = 0;
61 etype[j] = 0;
62 }
63 }
64 p->lineage_incr(ell,sat,etype);
65 }
66 tout += nd; ell += nd; sat += nd;
67 deme += nd; etype += nd;
68 tcur = time();
69 for (size_t j = 0; j < nd; j++) {
70 tout[j] = tcur;
71 sat[j] = ell[j] = 0;
72 deme[j] = j;
73 etype[j] = 3;
74 }
75}
void lineage_incr(int *incr, int *sat, int *etype) const
Definition lineages.cc:14
#define ell
Definition lbdp_pomp.c:11
Here is the call graph for this function:

◆ lineage_count() [2/2]

SEXP genealogy_t::lineage_count ( void ) const

lineage count and saturation

Definition at line 79 of file lineages.cc.

81{
82 SEXP tout, deme, ell, sat, etype, out, outn;
83 int nt = ntime(timezero())+1;
84 int nl = (ndeme()+1)*nt;
85 PROTECT(tout = NEW_NUMERIC(nl));
86 PROTECT(deme = NEW_INTEGER(nl));
87 PROTECT(ell = NEW_INTEGER(nl));
88 PROTECT(sat = NEW_INTEGER(nl));
89 PROTECT(etype = NEW_INTEGER(nl));
90 PROTECT(out = NEW_LIST(5));
91 PROTECT(outn = NEW_CHARACTER(5));
92 set_list_elem(out,outn,tout,"time",0);
93 set_list_elem(out,outn,deme,"deme",1);
94 set_list_elem(out,outn,ell,"lineages",2);
95 set_list_elem(out,outn,sat,"saturation",3);
96 set_list_elem(out,outn,etype,"event_type",4);
97 SET_NAMES(out,outn);
98 lineage_count(REAL(tout),INTEGER(deme),INTEGER(ell),
99 INTEGER(sat),INTEGER(etype));
100 UNPROTECT(7);
101 return out;
102}
SEXP lineage_count(void) const
lineage count and saturation
Definition lineages.cc:80
size_t ntime(slate_t t) const
Number of distinct timepoints.
Definition nodeseq.h:131
Here is the call graph for this function:
Here is the caller graph for this function:

◆ make_node()

node_t * genealogy_t::make_node ( name_t d = undeme)
inlineprivate

create a node holding its own green ball. this does not insert into the genealogy.

Definition at line 229 of file genealogy.h.

229 {
231 name_t u = unique();
232 node_t *p = new node_t(u,_time);
233 ball_t *g = new ball_t(p,u,green,d);
234 p->green_ball() = g;
235 p->insert(g);
236 return p;
237 };
bool check_genealogy_size(size_t grace=0) const
check the size of the genealogy (to prevent memory exhaustion).
Definition genealogy.h:214
Here is the call graph for this function:
Here is the caller graph for this function:

◆ migrate()

void genealogy_t::migrate ( ball_t * a,
slate_t t,
name_t d = 0 )
inline

movement into deme d

Definition at line 292 of file genealogy.h.

292 {
293 time() = t;
294 node_t *p = make_node(a->deme());
295 p->slate = time();
296 add(p,a);
297 a->deme() = d;
298 };
Here is the call graph for this function:

◆ ndeme() [1/2]

size_t & genealogy_t::ndeme ( void )
inline

number of demes

Definition at line 63 of file genealogy.h.

63 {
64 return _ndeme;
65 };

◆ ndeme() [2/2]

size_t genealogy_t::ndeme ( void ) const
inline

number of demes

Definition at line 59 of file genealogy.h.

59 {
60 return _ndeme;
61 };
Here is the caller graph for this function:

◆ newick()

string_t genealogy_t::newick ( bool extended = true) const

put genealogy at current time into Newick format.

Definition at line 93 of file newick.cc.

95{
96 return nodeseq_t::newick(time(),(ndeme() > 0),extended);
97}
string_t newick(slate_t t, bool showdeme, bool extended) const
put genealogy at time t into Newick format.
Definition newick.cc:79
Here is the call graph for this function:
Here is the caller graph for this function:

◆ nroot()

size_t genealogy_t::nroot ( void ) const
inline

number of roots

Definition at line 192 of file genealogy.h.

192 {
193 size_t n = 0;
194 for (const node_t *p : *this) {
195 if (p->holds_own()) n++;
196 }
197 return n;
198 };
Here is the call graph for this function:
Here is the caller graph for this function:

◆ nsample()

size_t genealogy_t::nsample ( void ) const
inline

number of samples

Definition at line 183 of file genealogy.h.

183 {
184 size_t n = 0;
185 for (const node_t *p : *this) {
186 if (p->holds(blue)) n++;
187 }
188 return n;
189 };
Here is the call graph for this function:
Here is the caller graph for this function:

◆ obscure()

genealogy_t & genealogy_t::obscure ( void )
inline

erase all deme information

Definition at line 321 of file genealogy.h.

321 {
322 // erase deme information from black balls.
323 pocket_t *blacks = colored(black);
324 while (!blacks->empty()) {
325 ball_t *a = *(blacks->begin());
326 a->deme() = undeme;
327 blacks->erase(a);
328 }
329 delete blacks;
330 // erase deme information from nodes.
331 for (node_t *p : *this) {
332 p->deme() = undeme;
333 }
334 // drop superfluous nodes (holding just one ball).
335 comb();
336 ndeme() = 0;
337 return *this;
338 };
static const name_t undeme
Definition ball.h:15
void comb(void)
Definition nodeseq.h:226
pocket_t * colored(color_t col) const
Get all balls of a color.
Definition nodeseq.h:121
Here is the call graph for this function:
Here is the caller graph for this function:

◆ operator+=()

genealogy_t & genealogy_t::operator+= ( const genealogy_t & other)

merges two genealogies, adjusting time, t0, and ndeme as needed

merge two genealogies:

  1. the node-sequences are merged;
  2. the root time advances as necessary;
  3. the current time retreats as necessary;
  4. the unique-name stack advances as necessary.

Definition at line 29 of file sum.cc.

31{
32 genealogy_t G = other;
33 slate_t t = time();
34 slate_t t0 = timezero();
35 t0 = (t0 > G.timezero()) ? t0 : G.timezero();
36 t = (t < G.time()) ? t : G.time();
38 merge(G,compare);
39 timezero() = (timezero() < G.timezero()) ? timezero() : G.timezero();
40 time() = (time() > G.time()) ? time() : G.time();
41 ndeme() = (ndeme() > G.ndeme()) ? ndeme() : G.ndeme();
42 _unique = G._unique;
43 curtail(t,t0);
44 return *this;
45}
void reuniqify(name_t shift)
shifts names to avoid overlap
Definition sum.cc:21
genealogy_t(double t0=R_NaReal, size_t ndeme=0)
Definition genealogy.h:102
void curtail(slate_t tnew, slate_t troot)
Definition curtail.cc:11
static bool compare(node_t *p, node_t *q)
Definition nodeseq.h:95
Here is the call graph for this function:

◆ operator=() [1/2]

genealogy_t & genealogy_t::operator= ( const genealogy_t & G)
inline

copy assignment operator

Definition at line 127 of file genealogy.h.

127 {
128 clean();
129 raw_t *o = new raw_t[G.bytesize()];
130 G >> o;
131 o >> *this;
132 delete[] o;
133 return *this;
134 };
Here is the call graph for this function:

◆ operator=() [2/2]

genealogy_t & genealogy_t::operator= ( genealogy_t && )
default

move assignment operator

Here is the call graph for this function:

◆ parse()

genealogy_t & genealogy_t::parse ( const string_t & s)

Parse a Newick string and create the indicated genealogy.

Definition at line 130 of file parse.cc.

132{
133 node_t *p = 0, *q;
134 slate_t tf = timezero();
135 string_t::const_reverse_iterator b = s.crbegin(), e = b;
136 bool open = false; // branch-string reading-frame open?
137 int stack = 0, sqstack = 0;
138 if (!s.empty() && *b != ';')
139 err("in '%s': invalid Newick format: no final semicolon.",__func__);
140 while (b != s.crend()) {
141 switch (*b) {
142 case ';': // root
143 if (stack != 0)
144 err("in '%s': invalid Newick: unbalanced parentheses.",__func__);
145 p = make_node();
146 p->slate = timezero();
147 push_front(p);
148 b++;
149 e = b;
150 open = true;
151 break;
152 case ')': // internal node
153 if (open) {
154 q = scan_branch(b.base(),e.base());
155 if (q->holds(black))
156 err("in '%s': 'type=extant' on internal node.",__func__);
157 q->slate += p->slate;
158 attach(p,q);
159 push_back(q);
160 tf = (q->slate > tf) ? q->slate : tf;
161 ndeme() = (ndeme() > q->deme()) ? ndeme() : q->deme();
162 p = q;
163 } else {
164 err("in '%s': invalid Newick: missing comma or semicolon.",__func__);
165 }
166 stack++;
167 b++;
168 e = b;
169 open = true;
170 break;
171 case '(': // tip node, eldest sister
172 if (open) {
173 q = scan_branch(b.base(),e.base());
174 q->slate += p->slate;
175 attach(p,q);
176 push_back(q);
177 tf = (q->slate > tf) ? q->slate : tf;
178 ndeme() = (ndeme() > q->deme()) ? ndeme() : q->deme();
179 }
180 p = p->parent();
181 b++;
182 stack--;
183 e = b;
184 open = false;
185 break;
186 case ',': // tip node, younger sister
187 if (stack <= 0)
188 err("in '%s': invalid Newick string: misplaced comma or unbalanced parentheses.",__func__);
189 if (open) {
190 q = scan_branch(b.base(),e.base());
191 q->slate += p->slate;
192 attach(p,q);
193 push_back(q);
194 tf = (q->slate > tf) ? q->slate : tf;
195 ndeme() = (ndeme() > q->deme()) ? ndeme() : q->deme();
196 }
197 b++;
198 e = b;
199 open = true;
200 break;
201 case ']': // skip metadata
202 sqstack++;
203 while (b != s.crend() && sqstack > 0) {
204 b++;
205 if (*b == ']') sqstack++;
206 if (*b == '[') sqstack--;
207 }
208 if (sqstack != 0)
209 err("in '%s': invalid Newick format: unbalanced square brackets.",__func__);
210 if (b != s.crend()) b++;
211 break;
212 default:
213 b++;
214 break;
215 }
216 }
217 if (stack != 0)
218 err("in '%s': invalid Newick format: unbalanced parentheses.",__func__);
219 if (open) {
220 q = scan_branch(b.base(),e.base());
221 q->slate += p->slate;
222 attach(p,q);
223 push_back(q);
224 tf = (q->slate > tf) ? q->slate : tf;
225 ndeme() = (ndeme() > q->deme()) ? ndeme() : q->deme();
226 }
227 time() = tf;
228 sort(); repair_tips(); drop_zlb(); weed();
229 return *this;
230}
void drop_zlb(void)
drop all zero-length branches
Definition genealogy.h:373
void repair_tips(void)
Definition parse.cc:9
node_t * scan_branch(string_t::const_iterator b, string_t::const_iterator e)
Definition parse.cc:101
void weed(void)
drop all dead roots
Definition nodeseq.h:214
void attach(node_t *p, node_t *q)
Definition nodeseq.h:174
Here is the call graph for this function:
Here is the caller graph for this function:

◆ prune()

genealogy_t & genealogy_t::prune ( void )
inline

prune the tree (drop all black balls)

Definition at line 310 of file genealogy.h.

310 {
311 pocket_t *blacks = colored(black);
312 while (!blacks->empty()) {
313 ball_t *b = *(blacks->begin());
314 blacks->erase(b);
315 drop(b);
316 }
317 delete blacks;
318 return *this;
319 };
Here is the call graph for this function:
Here is the caller graph for this function:

◆ repair_roots()

void genealogy_t::repair_roots ( void )
inlineprivate

roots are added at zero time if needed

Definition at line 393 of file genealogy.h.

393 {
394 node_nit j = begin();
395 while (j != end()) {
396 if ((*j)->holds_own() && (*j)->slate > timezero()) {
397 node_t *q = make_node();
398 q->slate = timezero();
399 attach(q,*j);
400 push_front(q);
401 }
402 j++;
403 }
404 sort();
405 };
std::list< node_t * >::iterator node_nit
Definition nodeseq.h:13
Here is the call graph for this function:

◆ repair_tips()

void genealogy_t::repair_tips ( void )
private
  • tips without descendants are reclassified as samples.
  • tips with black balls are swapped out.

Definition at line 8 of file parse.cc.

10{
11 for (node_t *p : *this) {
12 if (p->empty()) {
13 ball_t *b = new ball_t(p,p->uniq,blue,p->deme());
14 p->insert(b);
15 } else if (p->holds(black))
16 swap(p->last_ball(),p->green_ball());
17 }
18}
Here is the call graph for this function:
Here is the caller graph for this function:

◆ reuniqify()

void genealogy_t::reuniqify ( name_t shift)
private

shifts names to avoid overlap

Definition at line 20 of file sum.cc.

22{
23 this->_unique += shift;
24 for (node_t *p : *this) p->reuniqify(shift);
25}
void reuniqify(name_t shift)
shifts name to avoid overlap
Definition sum.cc:10
Here is the call graph for this function:
Here is the caller graph for this function:

◆ sample()

void genealogy_t::sample ( ball_t * a,
slate_t t )
inline

insert a sample node

Definition at line 273 of file genealogy.h.

273 {
274 time() = t;
275 node_t *p = make_node(a->deme());
276 ball_t *b = new ball_t (p,p->uniq,blue,a->deme());
277 p->insert(b);
278 p->slate = time();
279 add(p,a);
280 };
Here is the call graph for this function:

◆ sample_death()

void genealogy_t::sample_death ( ball_t * a,
slate_t t )
inline

insert a sample node and simultaneously terminate the lineage

Definition at line 282 of file genealogy.h.

282 {
283 time() = t;
284 node_t *p = make_node(a->deme());
285 ball_t *b = new ball_t (p,p->uniq,blue,a->deme());
286 p->insert(b);
287 p->slate = time();
288 add(p,a);
289 drop(a);
290 };
Here is the call graph for this function:

◆ sample_migrate()

void genealogy_t::sample_migrate ( ball_t * a,
slate_t t,
name_t d = 0 )
inline

insert a sample node and simultaneously migrate the lineage

Definition at line 300 of file genealogy.h.

300 {
301 time() = t;
302 node_t *p = make_node(a->deme());
303 ball_t *b = new ball_t (p,p->uniq,blue,a->deme());
304 p->insert(b);
305 p->slate = time();
306 add(p,a);
307 a->deme() = d;
308 };
Here is the call graph for this function:

◆ scan_branch()

node_t * genealogy_t::scan_branch ( string_t::const_iterator b,
string_t::const_iterator e )
private

Scan the branch (label+branch-length) string. This has format [&&PhyloPOMP deme=d type=s]s:f

Scan the branch string. This has format s[&&PhyloPOMP deme=d type=s]s:f

Definition at line 100 of file parse.cc.

103{
104 name_t deme = 0;
105 color_t col = green;
106 slate_t bl = 0;
107 if (b != e) {
108 std::smatch m;
109 if (std::regex_match(b,e,m,std::regex("^.*?\\[&&PhyloPOMP.+?deme=(\\w+).*?\\].*$")))
110 deme = scan_name(m[1].str());
111 if (std::regex_match(b,e,m,std::regex("^.*?\\[&&PhyloPOMP.+?type=(\\w+).*?\\].*$")))
112 col = scan_color(m[1].str());
113 if (std::regex_match(b,e,m,std::regex("^.*:([-+]?[0-9]*\\.?[0-9]+(?:[eE][-+]?[0-9]+)?)$")))
114 bl = scan_slate(m[1].str());
115 else
116 warn("in '%s': in branch-string '%s': no branch-length detected: assuming zero branch length.",
117 __func__,string_t(b,e).c_str());
118 }
119 node_t *q = make_node(deme);
120 if (col != green) {
121 ball_t *b = new ball_t(q,q->uniq,col,deme);
122 q->insert(b);
123 }
124 q->slate = bl;
125 return q;
126}
color_t
BALL COLORS.
Definition ball.h:12
#define warn(...)
Definition internal.h:19
static color_t scan_color(const std::string &s)
simple function for scanning the color
Definition parse.cc:78
static name_t scan_name(const string_t &s)
Definition parse.cc:52
static slate_t scan_slate(const string_t &s)
Definition parse.cc:25
Here is the call graph for this function:
Here is the caller graph for this function:

◆ structure()

SEXP genealogy_t::structure ( void ) const

R list description.

Definition at line 84 of file structure.cc.

86{
87 SEXP O, On, T0, Time, Nodes, Ndeme;
88 PROTECT(O = NEW_LIST(4));
89 PROTECT(On = NEW_CHARACTER(4));
90 PROTECT(Time = NEW_NUMERIC(1));
91 *REAL(Time) = double(time());
92 PROTECT(T0 = NEW_NUMERIC(1));
93 *REAL(T0) = double(timezero());
94 PROTECT(Ndeme = NEW_INTEGER(1));
95 *INTEGER(Ndeme) = int(ndeme());
96 PROTECT(Nodes = nodeseq_t::structure());
97 set_list_elem(O,On,Time,"time",0);
98 set_list_elem(O,On,T0,"t0",1);
99 set_list_elem(O,On,Ndeme,"ndeme",2);
100 set_list_elem(O,On,Nodes,"nodes",3);
101 SET_NAMES(O,On);
102 UNPROTECT(6);
103 return O;
104}
SEXP structure(void) const
R list description.
Definition structure.cc:71
Here is the call graph for this function:

◆ time() [1/2]

slate_t & genealogy_t::time ( void )
inline

view/set current time.

Definition at line 145 of file genealogy.h.

145 {
146 return _time;
147 };
Here is the caller graph for this function:

◆ time() [2/2]

slate_t genealogy_t::time ( void ) const
inline

view current time.

Definition at line 149 of file genealogy.h.

149 {
150 return _time;
151 };

◆ timezero() [1/2]

slate_t & genealogy_t::timezero ( void )
inline

view/set zero time.

Definition at line 153 of file genealogy.h.

153 {
154 return _t0;
155 };
Here is the caller graph for this function:

◆ timezero() [2/2]

slate_t genealogy_t::timezero ( void ) const
inline

get zero time.

Definition at line 157 of file genealogy.h.

157 {
158 return _t0;
159 };

◆ unique()

name_t genealogy_t::unique ( void )
inlineprivate

get the next unique name

Definition at line 43 of file genealogy.h.

43 {
44 name_t u = _unique;
45 _unique++;
46 return u;
47 };
Here is the caller graph for this function:

◆ valid()

void genealogy_t::valid ( void ) const
inline

check the validity of the genealogy.

Definition at line 212 of file genealogy.h.

212{};

◆ yaml()

string_t genealogy_t::yaml ( string_t tab = "") const

human/machine-readable info

Definition at line 63 of file yaml.cc.

65{
66 string_t o;
67 string_t t = tab + " ";
68 o = tab + "t0: " + std::to_string(timezero()) + "\n"
69 + tab + "time: " + std::to_string(time()) + "\n"
70 + tab + "ndeme: " + std::to_string(ndeme()) + "\n"
71 + tab + "nodes:\n" + nodeseq_t::yaml(tab);
72 return o;
73}
string_t yaml(string_t tab="") const
human/machine-readable info
Definition yaml.cc:52
Here is the call graph for this function:
Here is the caller graph for this function:

Friends And Related Symbol Documentation

◆ operator>> [1/2]

raw_t * operator>> ( const genealogy_t & G,
raw_t * o )
friend

binary serialization

Definition at line 76 of file genealogy.h.

76 {
77 name_t A[3]; A[0] = magic; A[1] = G._unique; A[2] = name_t(G.ndeme());
78 slate_t B[2]; B[0] = G.timezero(); B[1] = G.time();
79 memcpy(o,A,sizeof(A)); o += sizeof(A);
80 memcpy(o,B,sizeof(B)); o += sizeof(B);
81 return reinterpret_cast<const nodeseq_t&>(G) >> o;
82 };
static const name_t magic
Definition genealogy.h:38

◆ operator>> [2/2]

raw_t * operator>> ( raw_t * o,
genealogy_t & G )
friend

binary deserialization

Definition at line 84 of file genealogy.h.

84 {
85 G.clean();
86 name_t A[3];
87 slate_t B[2];
88 memcpy(A,o,sizeof(A)); o += sizeof(A);
89 memcpy(B,o,sizeof(B)); o += sizeof(B);
90 if (A[0] != magic)
91 err("in %s: corrupted genealogy serialization.",__func__);
92 G._unique = A[1]; G.ndeme() = size_t(A[2]);
93 G.timezero() = B[0]; G.time() = B[1];
94 return o >> reinterpret_cast<nodeseq_t&>(G);
95 };

Field Documentation

◆ _ndeme

size_t genealogy_t::_ndeme
private

The number of demes (excluding the undeme).

Definition at line 36 of file genealogy.h.

◆ _t0

slate_t genealogy_t::_t0
private

The initial time.

Definition at line 32 of file genealogy.h.

◆ _time

slate_t genealogy_t::_time
private

The current time.

Definition at line 34 of file genealogy.h.

◆ _unique

name_t genealogy_t::_unique
private

The next unique name.

Definition at line 30 of file genealogy.h.

◆ magic

const name_t genealogy_t::magic = 1123581321
staticprivate

Definition at line 38 of file genealogy.h.


The documentation for this class was generated from the following files: