phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
genealogy.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 // GENEALOGY class
3 
4 #ifndef _GENEALOGY_H_
5 #define _GENEALOGY_H_
6 
7 #include <string>
8 #include <cstring>
9 #include <utility>
10 #include <stdexcept>
11 
12 #include "nodeseq.h"
13 #include "inventory.h"
14 #include "internal.h"
15 
16 static const size_t MEMORY_MAX = (1<<28); // 256MB
17 
19 
22 class genealogy_t : public nodeseq_t {
23 
24 private:
25 
26  // GENEALOGY member data:
27  // - a counter of serial numbers
28  // - an initial time
29  // - the current time
30  // - a sequence of nodes
31 
39  size_t _ndeme;
40 
41  const static name_t magic = 1123581321;
42 
43 private:
44 
46  name_t unique (void) {
47  name_t u = _unique;
48  _unique++;
49  return u;
50  };
51 
53  void clean (void) {
54  _unique = 0;
55  _t0 = _time = R_NaReal;
56  };
57 
58 public:
59 
61  size_t ndeme (void) const {
62  return _ndeme;
63  };
65  size_t& ndeme (void) {
66  return _ndeme;
67  };
68 
69 public:
70 
71  // SERIALIZATION
73  size_t bytesize (void) const {
74  return 3*sizeof(name_t) +
75  2*sizeof(slate_t) + nodeseq_t::bytesize();
76  };
78  friend raw_t* operator>> (const genealogy_t& G, raw_t* o) {
79  name_t A[3]; A[0] = magic; A[1] = G._unique; A[2] = name_t(G._ndeme);
80  slate_t B[2]; B[0] = G._t0; B[1] = G._time;
81  memcpy(o,A,sizeof(A)); o += sizeof(A);
82  memcpy(o,B,sizeof(B)); o += sizeof(B);
83  return reinterpret_cast<const nodeseq_t&>(G) >> o;
84  };
86  friend raw_t* operator>> (raw_t* o, genealogy_t& G) {
87  G.clean();
88  name_t A[3];
89  slate_t B[2];
90  memcpy(A,o,sizeof(A)); o += sizeof(A);
91  memcpy(B,o,sizeof(B)); o += sizeof(B);
92  if (A[0] != magic)
93  err("in %s (%s line %d) corrupted genealogy serialization.",
94  __func__,__FILE__,__LINE__);
95  G._unique = A[1]; G._ndeme = size_t(A[2]);
96  G._t0 = B[0]; G._time = B[1];
97  return o >> reinterpret_cast<nodeseq_t&>(G);
98  };
99 
100 
101 public:
102  // CONSTRUCTORS
105  genealogy_t (double t0 = 0, name_t u = 0, size_t nd = 1) {
106  clean();
107  _ndeme = nd;
108  _unique = u;
109  _time = _t0 = slate_t(t0);
110  };
113  o >> *this;
114  };
116  genealogy_t (SEXP o) {
117  if (LENGTH(o)==0)
118  err("in %s (%s line %d): cannot deserialize a NULL.",
119  __func__,__FILE__,__LINE__);
120  PROTECT(o = AS_RAW(o));
121  RAW(o) >> *this;
122  UNPROTECT(1);
123  };
126  raw_t *o = new raw_t[G.bytesize()];
127  G >> o;
128  o >> *this;
129  delete[] o;
130  };
133  clean();
134  raw_t *o = new raw_t[G.bytesize()];
135  G >> o;
136  o >> *this;
137  delete[] o;
138  return *this;
139  };
141  genealogy_t (genealogy_t&&) = default;
145  ~genealogy_t (void) {
146  clean();
147  };
148 
150  slate_t& time (void) {
151  return _time;
152  };
154  slate_t time (void) const {
155  return _time;
156  };
158  slate_t timezero (void) const {
159  return _t0;
160  };
161 
162 public:
163 
171  void lineage_count (double *tout, int *deme,
172  int *ell, int *sat, int *etype) const {
173  slate_t tcur = timezero();
174  for (size_t j = 0; j < _ndeme; j++) {
175  tout[j] = tcur;
176  deme[j] = j+1;
177  sat[j] = ell[j] = 0;
178  etype[j] = 0;
179  }
180  for (node_it i = begin(); i != end(); i++) {
181  node_t *p = *i;
182  if (tcur < p->slate) {
183  tout += _ndeme; ell += _ndeme; sat += _ndeme;
184  deme += _ndeme; etype += _ndeme;
185  tcur = p->slate;
186  for (size_t j = 0; j < _ndeme; j++) {
187  tout[j] = tcur;
188  deme[j] = j+1;
189  ell[j] = (ell-_ndeme)[j];
190  sat[j] = 0;
191  etype[j] = 0;
192  }
193  }
194  p->lineage_incr(ell,sat,etype);
195  }
196  tout += _ndeme; ell += _ndeme; sat += _ndeme;
197  deme += _ndeme; etype += _ndeme;
198  tcur = time();
199  for (size_t j = 0; j < _ndeme; j++) {
200  tout[j] = tcur;
201  sat[j] = ell[j] = 0;
202  deme[j] = j+1;
203  etype[j] = 3;
204  }
205  };
207  SEXP lineage_count (void) const {
208  SEXP tout, deme, ell, sat, etype, out, outn;
209  int nt = ntime(timezero())+1;
210  int nl = _ndeme*nt;
211  PROTECT(tout = NEW_NUMERIC(nl));
212  PROTECT(deme = NEW_INTEGER(nl));
213  PROTECT(ell = NEW_INTEGER(nl));
214  PROTECT(sat = NEW_INTEGER(nl));
215  PROTECT(etype = NEW_INTEGER(nl));
216  PROTECT(out = NEW_LIST(5));
217  PROTECT(outn = NEW_CHARACTER(5));
218  set_list_elem(out,outn,tout,"time",0);
219  set_list_elem(out,outn,deme,"deme",1);
220  set_list_elem(out,outn,ell,"lineages",2);
221  set_list_elem(out,outn,sat,"saturation",3);
222  set_list_elem(out,outn,etype,"event_type",4);
223  SET_NAMES(out,outn);
224  lineage_count(REAL(tout),INTEGER(deme),INTEGER(ell),
225  INTEGER(sat),INTEGER(etype));
226  UNPROTECT(7);
227  return out;
228  };
229 
231  void gendat (double *tout, int *anc, int *lin,
232  int *sat, int *type,
233  int *index, int *child) const {
234  int m, n, k;
235  node_it i, j;
236  for (k = 0, n = 0, i = begin(); i != end(); i++, n++) {
237  node_t *p = *i;
238  assert(!p->holds(black)); // tree should be pruned first
239  tout[n] = p->slate;
240  if (p->is_root()) {
241  type[n] = 0; // root node
242  } else if (p->holds(blue)) {
243  type[n] = 1; // sample node
244  } else {
245  type[n] = 2; // internal node
246  }
247  lin[n] = p->lineage(); // 0-based indexing
248  sat[n] = p->nchildren();
249  index[n] = k;
250  k += sat[n];
251  child[n] = NA_INTEGER;
252  if (p->is_root()) {
253  anc[n] = n; // 0-based indexing
254  } else {
255  for (m = 0, j = begin(); j != i; j++, m++) {
256  node_t *q = *j;
257  if (p->parent()->uniq == q->uniq) {
258  anc[n] = m;
259  break;
260  }
261  }
262  }
263  }
264  tout[n] = time();
265  for (k = 0, n = 0, i = begin(); i != end(); i++, n++) {
266  node_t *p = *i;
267  j = i; j++;
268  for (m = n+1; j != end(); m++, j++) {
269  node_t *q = *j;
270  if (p->uniq == q->parent()->uniq) {
271  child[k++] = m;
272  }
273  }
274  }
275  };
277  SEXP gendat (void) const {
278  SEXP tout, anc, lin, sat, type, index, child, ns, nn;
279  SEXP out, outn;
280  size_t n = length();
281  PROTECT(tout = NEW_NUMERIC(n+1));
282  PROTECT(type = NEW_INTEGER(n));
283  PROTECT(lin = NEW_INTEGER(n));
284  PROTECT(sat = NEW_INTEGER(n));
285  PROTECT(index = NEW_INTEGER(n));
286  PROTECT(child = NEW_INTEGER(n));
287  PROTECT(anc = NEW_INTEGER(n));
288  PROTECT(ns = NEW_INTEGER(1));
289  PROTECT(nn = NEW_INTEGER(1));
290  PROTECT(out = NEW_LIST(9));
291  PROTECT(outn = NEW_CHARACTER(9));
292  set_list_elem(out,outn,tout,"nodetime",0);
293  set_list_elem(out,outn,type,"nodetype",1);
294  set_list_elem(out,outn,lin,"lineage",2);
295  set_list_elem(out,outn,sat,"saturation",3);
296  set_list_elem(out,outn,index,"index",4);
297  set_list_elem(out,outn,child,"child",5);
298  set_list_elem(out,outn,anc,"ancestor",6);
299  set_list_elem(out,outn,ns,"nsample",7);
300  set_list_elem(out,outn,nn,"nnode",8);
301  SET_NAMES(out,outn);
302  gendat(REAL(tout),INTEGER(anc),INTEGER(lin),INTEGER(sat),
303  INTEGER(type),INTEGER(index),INTEGER(child));
304  *INTEGER(ns) = nsample(); // number of samples
305  *INTEGER(nn) = length(); // number of nodes
306  UNPROTECT(11);
307  return out;
308  };
309 
311  size_t nsample (void) const {
312  size_t n = 0;
313  for (node_it k = cbegin(); k != cend(); k++) {
314  if ((*k)->holds(blue)) n++;
315  }
316  return n;
317  };
318 
319 public:
320 
322  SEXP structure (void) const {
323  SEXP O, On, T0, Time, Nodes;
324  PROTECT(O = NEW_LIST(3));
325  PROTECT(On = NEW_CHARACTER(3));
326  PROTECT(Time = NEW_NUMERIC(1));
327  *REAL(Time) = double(time());
328  PROTECT(T0 = NEW_NUMERIC(1));
329  *REAL(T0) = double(timezero());
330  PROTECT(Nodes = nodeseq_t::structure());
331  set_list_elem(O,On,Time,"time",0);
332  set_list_elem(O,On,T0,"t0",1);
333  set_list_elem(O,On,Nodes,"nodes",2);
334  SET_NAMES(O,On);
335  UNPROTECT(5);
336  return O;
337  };
338 
339 public:
340 
342  std::string describe (void) const {
343  std::string o = "t0 = " + std::to_string(double(timezero()))
344  + "\ntime = " + std::to_string(double(time())) + "\n"
346  return o;
347  };
348 
349 public:
350 
352  virtual std::string yaml (std::string tab = "") const {
353  std::string o;
354  std::string t = tab + " ";
355  o = tab + "t0: " + std::to_string(timezero()) + "\n"
356  + tab + "time: " + std::to_string(time()) + "\n"
357  + tab + "nodes:\n" + nodeseq_t::yaml(tab);
358  return o;
359  };
360 
361 public:
362 
364  std::string newick (void) const {
365  return nodeseq_t::newick(time());
366  };
367 
368 public:
369 
371  void valid (void) const {};
373  bool check_genealogy_size (size_t grace = 0) const {
374  static size_t maxq = MEMORY_MAX/(sizeof(node_t)+2*sizeof(ball_t));
375  bool ok = true;
376  if (size() > maxq+grace) {
377  err("maximum genealogy size exceeded!"); // #nocov
378  } else if (size() > maxq) {
379  ok = false; // #nocov
380  }
381  return ok;
382  };
383 
384 public:
385 
390  name_t u = unique();
391  node_t *p = new node_t(u,_time);
392  ball_t *g = new ball_t(p,u,green,d);
393  p->green_ball() = g;
394  p->insert(g);
395  return p;
396  };
397 
398 public:
399 
402  time() = t;
403  node_t *p = make_node(a->deme());
404  ball_t *b = new ball_t (p,p->uniq,black,d);
405  p->insert(b);
406  p->slate = time();
407  add(p,a);
408  return b;
409  };
412  ball_t *b = new ball_t(p,unique(),black,d);
413  p->insert(b);
414  return b;
415  };
417  void death (ball_t *a, slate_t t) {
418  time() = t;
419  drop(a);
420  };
423  time() = t;
424  node_t *p = make_node(d);
425  ball_t *b = new ball_t (p,p->uniq,black,d);
426  p->insert(b);
427  p->slate = timezero();
428  push_front(p);
429  return b;
430  };
432  void sample (ball_t* a, slate_t t) {
433  time() = t;
434  node_t *p = make_node(a->deme());
435  ball_t *b = new ball_t (p,p->uniq,blue,a->deme());
436  p->insert(b);
437  p->slate = time();
438  add(p,a);
439  };
441  void sample_death (ball_t* a, slate_t t) {
442  time() = t;
443  node_t *p = make_node(a->deme());
444  ball_t *b = new ball_t (p,p->uniq,blue,a->deme());
445  p->insert(b);
446  p->slate = time();
447  add(p,a);
448  drop(a);
449  };
451  ball_t* migrate (ball_t* a, slate_t t, name_t d = 0) {
452  time() = t;
453  node_t *p = make_node(a->deme());
454  p->slate = time();
455  add(p,a);
456  a->deme() = d;
457  return a;
458  };
461  std::pair<node_it, node_it> extant (void) const {
462  return std::pair<node_it,node_it>(cbegin(),cend());
463  };
465  genealogy_t& prune (void) {
466  pocket_t *blacks = colored(black);
467  while (!blacks->empty()) {
468  ball_t *b = *(blacks->begin());
469  blacks->erase(b);
470  drop(b);
471  }
472  delete blacks;
473  return *this;
474  };
477  // erase deme information from black balls.
478  pocket_t *blacks = colored(black);
479  while (!blacks->empty()) {
480  ball_t *a = *(blacks->begin());
481  a->deme() = 0;
482  blacks->erase(a);
483  }
484  delete blacks;
485  // erase deme information from nodes.
486  for (node_it i = begin(); i != end(); i++) {
487  (*i)->deme() = 0;
488  }
489  // drop superfluous nodes (holding just one ball).
490  comb();
491  _ndeme = 1;
492  return *this;
493  };
496  void curtail (slate_t tnew) {
497  if (!empty()) {
498  node_t *p = back();
499  while (!empty() && p->slate > tnew) {
500  ball_t *b;
501  while (p->size() > 1) {
502  b = p->last_ball();
503  switch (b->color) {
504  case black:
505  p->erase(b); delete b;
506  break;
507  case green: case blue: // #nocov
508  assert(0); // #nocov
509  break; // #nocov
510  }
511  }
512  b = p->last_ball();
513  switch (b->color) {
514  case blue:
515  b->color = black;
516  case black:
517  b->deme() = p->deme();
518  swap(b,p->green_ball());
519  case green:
520  destroy_node(p);
521  break;
522  }
523  if (!empty()) p = back();
524  }
525  }
526  if (tnew < time()) _time = tnew;
527  if (tnew < timezero()) _t0 = tnew;
528  };
535  reinterpret_cast<nodeseq_t&>(*this) += reinterpret_cast<nodeseq_t&>(G);
536  _t0 = (_t0 > G._t0) ? G._t0 : _t0;
537  _time = (_time < G._time) ? G._time : _time;
538  _unique = (_unique < G._unique) ? G._unique : _unique;
539  _ndeme = (_ndeme < G.ndeme()) ? G.ndeme() : _ndeme;
540  return *this;
541  };
542 
543 private:
544 
545  size_t scan_color (const std::string& s, color_t* col) const {
546  return 1;
547  };
550  size_t scan_label (const std::string& s, color_t* col,
551  name_t *deme, slate_t *time) const {
552  size_t n = s.size();
553  if (n < 1)
554  err("in '%s' (%s line %d): invalid Newick format: empty label.",__func__,__FILE__,__LINE__);
555  size_t sz, i = 1;
556  switch (s[0]) {
557  case 'o':
558  *col = black;
559  break;
560  case 'b':
561  *col = blue;
562  break;
563  case 'g': case 'm':
564  *col = green;
565  break;
566  default:
567  err("in '%s' (%s line %d): invalid Newick label: expected one of 'b','g','m', or 'o', got '%c'.",
568  __func__,__FILE__,__LINE__,s[0]);
569  break;
570  }
571  while (i < n && s[i] == '_') i++;
572  if (i == n)
573  err("in '%s': invalid Newick format: premature termination.",__func__);
574  if (s[i] == '(' || s[i] == ')' || s[i] == ',' || s[i] == ';')
575  err("in '%s' (%s line %d): invalid Newick format.",__func__,__FILE__,__LINE__);
576  if (s[i] == ':') {
577  *deme = 0;
578  } else {
579  try {
580  *deme = name_t(stoi(s.substr(i),&sz));
581  i += sz;
582  }
583  catch (const std::invalid_argument& e) {
584  err("in '%s' (%s line %d): invalid Newick format: deme should be indicated with an integer.",
585  __func__,__FILE__,__LINE__);
586  }
587  catch (const std::out_of_range& e) {
588  err("in '%s': invalid Newick format: deme out of range.",__func__);
589  }
590  catch (const std::exception& e) {
591  err("in '%s': parsing deme label: %s.",__func__,e.what());
592  }
593  catch (...) {
594  err("in '%s': other deme-parsing error.",__func__);
595  }
596  }
597  // skip to branch length
598  while (i < n && s[i] != ':' &&
599  s[i] != '(' && s[i] != ')' && s[i] != ',' && s[i] != ';') i++;
600  if (i == n || s[i] != ':')
601  err("in '%s': invalid Newick format: missing or invalid branch length.",__func__);
602  i++;
603  try {
604  *time = slate_t(stod(s.substr(i),&sz));
605  }
606  catch (const std::invalid_argument& e) {
607  err("in '%s': invalid Newick format: branch length should be a non-negative decimal number.",__func__);
608  }
609  catch (const std::out_of_range& e) {
610  err("in '%s': invalid Newick format: branch length out of range.",__func__);
611  }
612  catch (const std::exception& e) {
613  err("in '%s': parsing branch-length: %s.",__func__,e.what());
614  }
615  catch (...) {
616  err("in '%s': other branch-length parsing error.",__func__);
617  }
618  i += sz;
619  return i;
620  };
623  size_t scan_ball (const std::string& s, const slate_t t0, node_t *p) {
624  color_t col;
625  name_t deme;
626  slate_t t;
627  ball_t *b;
628  size_t i = scan_label(s,&col,&deme,&t);
629  t += t0;
630  _time = (_time < t) ? t : _time;
631  _ndeme = (_ndeme <= deme) ? deme+1 : _ndeme;
632  if (col != black) err("in '%s': bad Newick string (1)",__func__);
633  if (p == 0) err("in '%s': bad Newick string (2)",__func__);
634  b = new ball_t(p,unique(),col,deme);
635  p->insert(b);
636  return i;
637  };
639  size_t scan_node (const std::string& s, const slate_t t0, node_t **q) {
640  size_t i;
641  color_t col;
642  name_t deme;
643  slate_t t;
644  name_t u = unique();
645  i = scan_label(s,&col,&deme,&t);
646  t += t0;
647  _ndeme = (_ndeme <= deme) ? deme+1 : _ndeme;
648  _time = (_time < t) ? t : _time;
649  node_t *p = new node_t(u,t);
650  ball_t *g = new ball_t(p,u,green,deme);
651  p->green_ball() = g;
652  p->insert(g);
653  if (col==blue) {
654  ball_t *b = new ball_t(p,p->uniq,blue,deme);
655  p->insert(b);
656  }
657  push_back(p);
658  *q = p;
659  return i;
660  };
663  size_t scan_tree (const std::string& s,
664  const slate_t t0, node_t **root) {
665  size_t n = s.size();
666  size_t i = 1, j = 1, k;
667  size_t stack = 1;
668  while (j < n && stack > 0) {
669  switch (s[j]) {
670  case '(':
671  stack++;
672  break;
673  case ')':
674  stack--;
675  break;
676  default:
677  break;
678  }
679  j++;
680  }
681  if (stack > 0)
682  err("in '%s': premature end of Newick string.",__func__);
683  node_t *p = 0;
684  k = j;
685  k += scan_node(s.substr(j),t0,&p);
686  parse(s.substr(i,j-i-1),p->slate,p);
687  *root = p;
688  return k;
689  };
690 
691 public:
692 
694  genealogy_t& parse (const std::string& s, slate_t t0, node_t *p = 0) {
695  size_t i = 0;
696  size_t n = s.size();
697  while (i < n) {
698  switch (s[i]) {
699  case '(':
700  {
701  genealogy_t G(t0,unique());
702  node_t *q = 0;
703  i += G.scan_tree(s.substr(i),t0,&q);
704  *this += G;
705  if (p != 0) {
706  ball_t *g = q->green_ball();
707  q->erase(g); p->insert(g); g->holder() = p;
708  }
709  }
710  break;
711  case 'b':
712  {
713  node_t *q = 0;
714  i += scan_node(s.substr(i),t0,&q);
715  if (p != 0) {
716  ball_t *g = q->green_ball();
717  q->erase(g); p->insert(g); g->holder() = p;
718  }
719  }
720  break;
721  case 'o':
722  i += scan_ball(s.substr(i),t0,p);
723  break;
724  case ',': case ';': case ' ': case '\n': case '\t':
725  i++;
726  break;
727  case ')':
728  err("in '%s': invalid Newick string: unbalanced parentheses.",__func__);
729  break;
730  default:
731  err("in '%s': invalid Newick string.",__func__);
732  break;
733  }
734  }
735  sort();
736  return *this;
737  };
738 };
739 
740 #endif
color_t
BALL COLORS.
Definition: ball.h:14
@ green
Definition: ball.h:14
@ black
Definition: ball.h:14
@ blue
Definition: ball.h:14
Balls function as pointers.
Definition: ball.h:29
name_t deme(void) const
view deme
Definition: ball.h:86
node_t * holder(void) const
in whose pocket do I lie?
Definition: ball.h:109
color_t color
Definition: ball.h:38
Encodes a genealogy.
Definition: genealogy.h:22
ball_t * birth(ball_t *a, slate_t t, name_t d)
birth into deme d
Definition: genealogy.h:401
slate_t time(void) const
view current time.
Definition: genealogy.h:154
size_t nsample(void) const
number of samples
Definition: genealogy.h:311
~genealogy_t(void)
destructor
Definition: genealogy.h:145
ball_t * birth(node_t *p, name_t d)
birth of second or subsequent sibling into deme d
Definition: genealogy.h:411
SEXP structure(void) const
R list description.
Definition: genealogy.h:322
genealogy_t & operator+=(genealogy_t &G)
Definition: genealogy.h:534
void clean(void)
clean up
Definition: genealogy.h:53
std::string newick(void) const
put genealogy at current time into Newick format.
Definition: genealogy.h:364
genealogy_t & operator=(const genealogy_t &G)
copy assignment operator
Definition: genealogy.h:132
virtual std::string yaml(std::string tab="") const
machine-readable info
Definition: genealogy.h:352
genealogy_t(double t0=0, name_t u=0, size_t nd=1)
Definition: genealogy.h:105
void valid(void) const
check the validity of the genealogy.
Definition: genealogy.h:371
genealogy_t(SEXP o)
constructor from RAW SEXP (containing binary serialization)
Definition: genealogy.h:116
void death(ball_t *a, slate_t t)
death
Definition: genealogy.h:417
size_t scan_ball(const std::string &s, const slate_t t0, node_t *p)
Definition: genealogy.h:623
SEXP lineage_count(void) const
lineage count and saturation
Definition: genealogy.h:207
size_t bytesize(void) const
size of serialized binary form
Definition: genealogy.h:73
genealogy_t & prune(void)
prune the tree (drop all black balls)
Definition: genealogy.h:465
genealogy_t & parse(const std::string &s, slate_t t0, node_t *p=0)
Parse a Newick string and create the indicated genealogy.
Definition: genealogy.h:694
std::string describe(void) const
human-readable info
Definition: genealogy.h:342
std::pair< node_it, node_it > extant(void) const
Definition: genealogy.h:461
size_t ndeme(void) const
number of demes
Definition: genealogy.h:61
genealogy_t(const genealogy_t &G)
copy constructor
Definition: genealogy.h:125
bool check_genealogy_size(size_t grace=0) const
check the size of the genealogy (to prevent memory exhaustion).
Definition: genealogy.h:373
size_t scan_label(const std::string &s, color_t *col, name_t *deme, slate_t *time) const
Definition: genealogy.h:550
slate_t timezero(void) const
get zero time.
Definition: genealogy.h:158
slate_t & time(void)
view/set current time.
Definition: genealogy.h:150
slate_t _time
The current time.
Definition: genealogy.h:37
size_t scan_node(const std::string &s, const slate_t t0, node_t **q)
Scan the Newick string and create the indicated node.
Definition: genealogy.h:639
SEXP gendat(void) const
nodelist in data-frame format
Definition: genealogy.h:277
size_t _ndeme
The number of demes.
Definition: genealogy.h:39
ball_t * migrate(ball_t *a, slate_t t, name_t d=0)
movement into deme d
Definition: genealogy.h:451
node_t * make_node(name_t d)
Definition: genealogy.h:388
void sample_death(ball_t *a, slate_t t)
insert a sample node and simultaneously terminate the lineage
Definition: genealogy.h:441
genealogy_t(raw_t *o)
constructor from serialized binary form
Definition: genealogy.h:112
genealogy_t & obscure(void)
erase all deme information
Definition: genealogy.h:476
ball_t * graft(slate_t t, name_t d)
graft a new lineage into deme d
Definition: genealogy.h:422
size_t scan_tree(const std::string &s, const slate_t t0, node_t **root)
Definition: genealogy.h:663
static const name_t magic
Definition: genealogy.h:41
name_t unique(void)
get the next unique name
Definition: genealogy.h:46
size_t scan_color(const std::string &s, color_t *col) const
Definition: genealogy.h:545
size_t & ndeme(void)
number of demes
Definition: genealogy.h:65
void sample(ball_t *a, slate_t t)
insert a sample node
Definition: genealogy.h:432
void lineage_count(double *tout, int *deme, int *ell, int *sat, int *etype) const
Definition: genealogy.h:171
slate_t _t0
The initial time.
Definition: genealogy.h:35
void gendat(double *tout, int *anc, int *lin, int *sat, int *type, int *index, int *child) const
nodelist in data-frame format
Definition: genealogy.h:231
name_t _unique
The next unique name.
Definition: genealogy.h:33
friend raw_t * operator>>(const genealogy_t &G, raw_t *o)
binary serialization
Definition: genealogy.h:78
genealogy_t(genealogy_t &&)=default
move constructor
void curtail(slate_t tnew)
Definition: genealogy.h:496
Encodes a genealogical node.
Definition: node.h:23
name_t lineage(void) const
view lineage
Definition: node.h:104
node_t * parent(void) const
Definition: node.h:115
bool is_root(void) const
Definition: node.h:121
ball_t * green_ball(void) const
pointer to my green ball
Definition: node.h:88
name_t deme(void) const
view deme
Definition: node.h:96
name_t uniq
Definition: node.h:30
slate_t slate
Definition: node.h:35
void lineage_incr(int *incr, int *sat, int *etype) const
Definition: node.h:151
int nchildren(void) const
number of descendants
Definition: node.h:131
A sequence of nodes.
Definition: nodeseq.h:19
std::string describe(void) const
human-readable info
Definition: nodeseq.h:264
virtual std::string yaml(std::string tab="") const
human- & machine-readable info
Definition: nodeseq.h:272
size_t ntime(slate_t t) const
Number of distinct timepoints.
Definition: nodeseq.h:139
void comb(void)
Definition: nodeseq.h:210
void sort(void)
order nodes in order of increasing time
Definition: nodeseq.h:111
void destroy_node(node_t *p)
remove a dead root node
Definition: nodeseq.h:203
SEXP structure(void) const
R list description.
Definition: nodeseq.h:281
size_t bytesize(void) const
size of serialized binary form
Definition: nodeseq.h:40
void drop(ball_t *a)
Definition: nodeseq.h:187
void swap(ball_t *a, ball_t *b)
swap balls a and b, wherever they lie
Definition: nodeseq.h:169
size_t length(void) const
Number of nodes in the sequence.
Definition: nodeseq.h:150
pocket_t * colored(color_t col) const
Get all balls of a color.
Definition: nodeseq.h:129
std::string newick(slate_t t) const
put genealogy at time t into Newick format.
Definition: nodeseq.h:292
void add(node_t *p, ball_t *a)
Definition: nodeseq.h:179
A pocket is a set of balls.
Definition: pocket.h:32
ball_t * last_ball(void) const
retrieve the last ball
Definition: pocket.h:125
bool holds(ball_t *b) const
does this node hold the given ball?
Definition: pocket.h:112
static const size_t MEMORY_MAX
Definition: genealogy.h:16
Rbyte raw_t
Definition: internal.h:43
static int set_list_elem(SEXP list, SEXP names, SEXP element, const char *name, int pos)
Definition: internal.h:67
size_t name_t
Definition: internal.h:45
#define err(...)
Definition: internal.h:18
double slate_t
Definition: internal.h:44
static int deme
Definition: lbdp.cc:7
#define ell
Definition: lbdp_pomp.c:10
#define n
Definition: lbdp_pomp.c:8
std::list< node_t * >::const_iterator node_it
Definition: nodeseq.h:14