phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
nodeseq.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 // NODE SEQUENCE CLASS
3 
4 #ifndef _NODESEQ_H_
5 #define _NODESEQ_H_
6 
7 #include <list>
8 #include <unordered_map>
9 #include <string>
10 #include <cstring>
11 #include "node.h"
12 #include "internal.h"
13 
14 typedef typename std::list<node_t*>::const_iterator node_it;
15 typedef typename std::list<node_t*>::iterator node_nit;
16 typedef typename std::list<node_t*>::const_reverse_iterator node_rev_it;
17 
19 class nodeseq_t : public std::list<node_t*> {
20 
21 private:
22 
24  void clean (void) {
25  for (node_it i = begin(); i != end(); i++) delete *i;
26  clear();
27  };
28 
29 public:
30 
32  ~nodeseq_t (void) {
33  clean();
34  };
35 
36 public:
37 
38  // SERIALIZATION
40  size_t bytesize (void) const {
41  size_t s = sizeof(size_t);
42  for (node_it i = begin(); i != end(); i++)
43  s += (*i)->bytesize();
44  return s;
45  };
47  friend raw_t* operator>> (const nodeseq_t& G, raw_t* o) {
48  size_t nnode = G.size();
49  memcpy(o,&nnode,sizeof(size_t)); o += sizeof(size_t);
50  for (node_it i = G.begin(); i != G.end(); i++) {
51  o = (**i >> o);
52  }
53  return o;
54  };
56  friend raw_t* operator>> (raw_t* o, nodeseq_t& G) {
57  G.clean();
58  std::unordered_map<name_t,node_t*> node_names;
59  std::unordered_map<name_t,ball_t*> ball_names;
60  size_t nnode = 0;
61  memcpy(&nnode,o,sizeof(size_t)); o += sizeof(size_t);
62  node_names.reserve(nnode);
63  ball_names.reserve(nnode);
64  for (size_t i = 0; i < nnode; i++) {
65  node_t *p = new node_t();
66  o = (o >> *p);
67  G.push_back(p);
68  node_names.insert({p->uniq,p});
69  }
70  for (node_it i = G.begin(); i != G.end(); i++) {
71  (*i)->repair_owners(node_names,&ball_names);
72  }
73  G.repair_owners(ball_names);
74  G.trace_lineages();
75  return o;
76  };
77 
78 private:
79 
82  void repair_owners (std::unordered_map<name_t,ball_t*>& names) {
83  std::unordered_map<name_t,ball_t*>::const_iterator n;
84  for (node_it i = begin(); i != end(); i++) {
85  node_t *p = *i;
86  n = names.find(p->uniq);
87  assert(n != names.end());
88  ball_t *b = n->second;
89  p->green_ball() = b;
90  }
91  };
92 
93 private:
94 
97  static bool compare (node_t* p, node_t* q) {
98  return (p->slate < q->slate) ||
99  ((p->slate == q->slate) && (p->uniq < q->uniq));
100  };
101 
102 public:
103 
106  merge(other,compare);
107  return *this;
108  };
109 
111  void sort (void) {
112  std::list<node_t*>::sort(compare);
113  };
114 
115 private:
116 
118  slate_t dawn (void) const {
119  return (empty()) ? R_NaReal : front()->slate;
120  };
122  slate_t dusk (void) const {
123  return (empty()) ? R_NaReal : back()->slate;
124  }
125 
126 public:
127 
129  pocket_t* colored (color_t col) const {
130  pocket_t *p = new pocket_t;
131  for (node_it i = begin(); i != end(); i++) {
132  for (ball_it j = (*i)->begin(); j != (*i)->end(); j++) {
133  if ((*j)->is(col)) p->insert(*j);
134  }
135  }
136  return p;
137  };
139  size_t ntime (slate_t t) const {
140  size_t count = 1;
141  for (node_it i = begin(); i != end(); i++) {
142  if (t < (*i)->slate) {
143  t = (*i)->slate;
144  count++;
145  }
146  }
147  return count;
148  };
150  size_t length (void) const {
151  size_t count = 0;
152  for (node_it i = begin(); i != end(); i++) count++;
153  return count;
154  };
156  node_t *position (int n) {
157  int i = 0;
158  node_it k = cbegin();
159  while (i < n && k != cend()) {
160  i++; k++;
161  }
162  assert(k != cend());
163  return *k;
164  };
165 
166 public:
167 
169  void swap (ball_t *a, ball_t *b) {
170  node_t *p = a->holder();
171  node_t *q = b->holder();
172  if (p != q) {
173  p->erase(a); q->insert(a); a->holder() = q;
174  q->erase(b); p->insert(b); b->holder() = p;
175  }
176  };
179  void add (node_t *p, ball_t *a) {
180  swap(a,p->green_ball());
181  p->deme() = a->deme();
182  push_back(p);
183  };
187  void drop (ball_t *a) {
188  assert(a->is(black));
189  node_t *p = a->holder();
190  if (p->size() > 1) {
191  p->erase(a);
192  delete a;
193  if (p->dead_root()) { // remove isolated root
194  destroy_node(p);
195  }
196  } else {
197  swap(a,p->green_ball());
198  destroy_node(p);
199  drop(a); // recurse
200  }
201  };
203  void destroy_node (node_t *p) {
204  assert(p->dead_root());
205  remove(p);
206  delete p;
207  };
210  void comb (void) {
211  for (node_rev_it i = crbegin(); i != crend(); i++) {
212  if ((*i)->size() == 1 && (*i)->holds(green)) {
213  swap((*i)->last_ball(),(*i)->green_ball());
214  }
215  }
216  node_nit j = begin();
217  while (j != end()) {
218  if ((*j)->dead_root()) {
219  destroy_node(*(j++));
220  } else {
221  j++;
222  }
223  }
224  };
225 
226 private:
227 
231  void trace_lineage (ball_t *b, name_t u) {
232  node_t *p = b->holder();
233  while (p->lineage() == null_lineage) {
234  p->lineage() = u;
235  p = p->parent();
236  }
237  };
238 
239 public:
240 
244  void trace_lineages (void) {
245  // we trace each lineage in turn.
246  // because we move from early to late,
247  // the order is guaranteed to be valid.
248  name_t u = 0;
249  for (node_it i = begin(); i != end(); i++) {
250  node_t *p = *i;
251  for (ball_it j = p->begin(); j != p->end(); j++) {
252  ball_t *b = *j;
253  if (b->color==blue) {
254  trace_lineage(b,u);
255  u++;
256  }
257  }
258  }
259  };
260 
261 public:
262 
264  std::string describe (void) const {
265  std::string o = "";
266  for (node_it p = begin(); p != end(); p++) {
267  o += (*p)->describe();
268  }
269  return o;
270  };
272  virtual std::string yaml (std::string tab = "") const {
273  std::string o = "";
274  std::string t = tab + " ";
275  for (node_it p = begin(); p != end(); p++) {
276  o += tab + "- " + (*p)->yaml(t);
277  }
278  return o;
279  };
281  SEXP structure (void) const {
282  SEXP Nodes;
283  PROTECT(Nodes = NEW_LIST(size()));
284  int k = 0;
285  for (node_it i = begin(); i != end(); i++) {
286  SET_ELEMENT(Nodes,k++,(*i)->structure());
287  }
288  UNPROTECT(1);
289  return Nodes;
290  };
292  std::string newick (slate_t t) const {
293  slate_t te = dawn();
294  std::string o = "";
295  for (node_it i = begin(); i != end(); i++) {
296  if ((*i)->is_root()) {
297  o += (*i)->newick(t,te) + ";";
298  }
299  }
300  return o;
301  };
302 };
303 
304 #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
bool is(color_t c) const
is a given ball of the given color?
Definition: ball.h:117
node_t * holder(void) const
in whose pocket do I lie?
Definition: ball.h:109
color_t color
Definition: ball.h:38
Encodes a genealogical node.
Definition: node.h:23
name_t lineage(void) const
view lineage
Definition: node.h:104
bool dead_root(void) const
Definition: node.h:124
node_t * parent(void) const
Definition: node.h:115
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
A sequence of nodes.
Definition: nodeseq.h:19
slate_t dawn(void) const
Earliest time in the sequence.
Definition: nodeseq.h:118
std::string describe(void) const
human-readable info
Definition: nodeseq.h:264
void clean(void)
clean up: delete all nodes, reset globals
Definition: nodeseq.h:24
void trace_lineage(ball_t *b, name_t u)
Definition: nodeseq.h:231
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
node_t * position(int n)
traverse to nth node, retrieve pointer
Definition: nodeseq.h:156
void comb(void)
Definition: nodeseq.h:210
void sort(void)
order nodes in order of increasing time
Definition: nodeseq.h:111
nodeseq_t & operator+=(nodeseq_t &other)
merge two node sequences
Definition: nodeseq.h:105
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 repair_owners(std::unordered_map< name_t, ball_t * > &names)
Definition: nodeseq.h:82
static bool compare(node_t *p, node_t *q)
Definition: nodeseq.h:97
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
void trace_lineages(void)
Definition: nodeseq.h:244
slate_t dusk(void) const
Latest time in the sequence.
Definition: nodeseq.h:122
pocket_t * colored(color_t col) const
Get all balls of a color.
Definition: nodeseq.h:129
~nodeseq_t(void)
destructor
Definition: nodeseq.h:32
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
friend raw_t * operator>>(const nodeseq_t &G, raw_t *o)
binary serialization
Definition: nodeseq.h:47
A pocket is a set of balls.
Definition: pocket.h:32
Rbyte raw_t
Definition: internal.h:43
size_t name_t
Definition: internal.h:45
double slate_t
Definition: internal.h:44
#define n
Definition: lbdp_pomp.c:8
static const name_t null_lineage
Definition: node.h:13
std::list< node_t * >::const_reverse_iterator node_rev_it
Definition: nodeseq.h:16
std::list< node_t * >::const_iterator node_it
Definition: nodeseq.h:14
std::list< node_t * >::iterator node_nit
Definition: nodeseq.h:15
std::set< ball_t *, ball_order >::const_iterator ball_it
Definition: pocket.h:26