genealogy.h

Go to the documentation of this file.

// -*- C++ -*-
// GENEALOGY class
 
#ifndef _GENEALOGY_H_
#define _GENEALOGY_H_
 
#include <utility>
#include <stdexcept>
#include <vector>
 
#include "nodeseq.h"
#include "internal.h"
 
static const size_t MEMORY_MAX = (1<<28); // 256MB


class genealogy_t : public nodeseq_t {
 
private:
 
  // GENEALOGY member data:
  // - a counter of serial numbers
  // - an initial time
  // - the current time
  // - a sequence of nodes

  name_t _unique;
  slate_t _t0;
  slate_t _time;
  size_t _ndeme;
 
  const static name_t magic = 1123581321;
 
private:

  name_t unique (void) {
    name_t u = _unique;
    _unique++;
    return u;
  };

  void clean (void) {
    _unique = 0;
    _ndeme = 0;
    _t0 = _time = R_NaReal;
  };
 
public:

  size_t ndeme (void) const {
    return _ndeme;
  };

  size_t& ndeme (void) {
    return _ndeme;
  };
 
public:
 
  // SERIALIZATION
  size_t bytesize (void) const {
    return 3*sizeof(name_t) +
      2*sizeof(slate_t) + nodeseq_t::bytesize();
  };

  friend raw_t* operator>> (const genealogy_t& G, raw_t* o) {
    name_t A[3]; A[0] = magic; A[1] = G._unique; A[2] = name_t(G.ndeme());
    slate_t B[2]; B[0] = G.timezero(); B[1] = G.time();
    memcpy(o,A,sizeof(A)); o += sizeof(A);
    memcpy(o,B,sizeof(B)); o += sizeof(B);
    return reinterpret_cast<const nodeseq_t&>(G) >> o;
  };

  friend raw_t* operator>> (raw_t* o, genealogy_t& G) {
    G.clean();
    name_t A[3];
    slate_t B[2];
    memcpy(A,o,sizeof(A)); o += sizeof(A);
    memcpy(B,o,sizeof(B)); o += sizeof(B);
    if (A[0] != magic)
      err("in %s: corrupted genealogy serialization.",__func__);
    G._unique = A[1]; G.ndeme() = size_t(A[2]);
    G.timezero() = B[0]; G.time() = B[1];
    return o >> reinterpret_cast<nodeseq_t&>(G);
  };
 
public:
  // CONSTRUCTORS
  genealogy_t (double t0 = R_NaReal, size_t ndeme = 0) {
    clean();
    _time = _t0 = slate_t(t0);
    _ndeme = ndeme;
  };

  genealogy_t (raw_t *o) {
    o >> *this;
  };

  genealogy_t (SEXP o) {
    if (LENGTH(o)==0)
      err("in %s: cannot deserialize a NULL.",__func__);
    PROTECT(o = AS_RAW(o));
    RAW(o) >> *this;
    UNPROTECT(1);
  };

  genealogy_t (const genealogy_t& G) {
    raw_t *o = new raw_t[G.bytesize()];
    G >> o;
    o >> *this;
    delete[] o;
  };

  genealogy_t& operator= (const genealogy_t& G) {
    clean();
    raw_t *o = new raw_t[G.bytesize()];
    G >> o;
    o >> *this;
    delete[] o;
    return *this;
  };

  genealogy_t (genealogy_t&&) = default;
  genealogy_t& operator= (genealogy_t&&) = default;
  ~genealogy_t (void) {
    clean();
  };

  slate_t& time (void) {
    return _time;
  };

  slate_t time (void) const {
    return _time;
  };

  slate_t& timezero (void) {
    return _t0;
  };

  slate_t timezero (void) const {
    return _t0;
  };
 
public:

  void lineage_count (double *tout, int *deme,
                      int *ell, int *sat, int *etype) const;
  SEXP lineage_count (void) const;

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

  size_t nsample (void) const {
    size_t n = 0;
    for (const node_t *p : *this) {
      if (p->holds(blue)) n++;
    }
    return n;
  };

  size_t nroot (void) const {
    size_t n = 0;
    for (const node_t *p : *this) {
      if (p->is_root()) n++;
    }
    return n;
  };
 
public:

  string_t yaml (string_t tab = "") const;
  SEXP structure (void) const;
  string_t newick (bool extended = true) const;
 
public:

  void valid (void) const {};
  bool check_genealogy_size (size_t grace = 0) const {
    static size_t maxq = MEMORY_MAX/(sizeof(node_t)+2*sizeof(ball_t));
    bool ok = true;
    if (size() > maxq+grace) {
      err("maximum genealogy size exceeded!"); // #nocov
    } else if (size() > maxq) {
      ok = false;               // #nocov
    }
    return ok;
  };
 
private:

  node_t* make_node (name_t d = undeme) {
    check_genealogy_size(0);
    name_t u = unique();
    node_t *p = new node_t(u,_time);
    ball_t *g = new ball_t(p,u,green,d);
    p->green_ball() = g;
    p->insert(g);
    return p;
  };
 
public:

  ball_t* birth (ball_t* a, slate_t t, name_t d) {
    time() = t;
    node_t *p = make_node(a->deme());
    ball_t *b = new ball_t (p,p->uniq,black,d);
    p->insert(b);
    p->slate = time();
    add(p,a);
    return b;
  };

  ball_t* birth (node_t* p, name_t d) {
    ball_t *b = new ball_t(p,unique(),black,d);
    p->insert(b);
    return b;
  };

  void death (ball_t *a, slate_t t) {
    time() = t;
    drop(a);
  };

  ball_t* graft (slate_t t, name_t d) {
    time() = t;
    node_t *p = make_node();
    ball_t *b = new ball_t (p,p->uniq,black,d);
    p->insert(b);
    p->slate = timezero();
    push_front(p);
    return b;
  };

  void sample (ball_t* a, slate_t t) {
    time() = t;
    node_t *p = make_node(a->deme());
    ball_t *b = new ball_t (p,p->uniq,blue,a->deme());
    p->insert(b);
    p->slate = time();
    add(p,a);
  };

  void sample_death (ball_t* a, slate_t t) {
    time() = t;
    node_t *p = make_node(a->deme());
    ball_t *b = new ball_t (p,p->uniq,blue,a->deme());
    p->insert(b);
    p->slate = time();
    add(p,a);
    drop(a);
  };

  void migrate (ball_t* a, slate_t t, name_t d = 0) {
    time() = t;
    node_t *p = make_node(a->deme());
    p->slate = time();
    add(p,a);
    a->deme() = d;
  };

  void sample_migrate (ball_t* a, slate_t t, name_t d = 0) {
    time() = t;
    node_t *p = make_node(a->deme());
    ball_t *b = new ball_t (p,p->uniq,blue,a->deme());
    p->insert(b);
    p->slate = time();
    add(p,a);
    a->deme() = d;
  };

  genealogy_t& prune (void) {
    pocket_t *blacks = colored(black);
    while (!blacks->empty()) {
      ball_t *b = *(blacks->begin());
      blacks->erase(b);
      drop(b);
    }
    delete blacks;
    return *this;
  };

  genealogy_t& obscure (void) {
    // erase deme information from black balls.
    pocket_t *blacks = colored(black);
    while (!blacks->empty()) {
      ball_t *a = *(blacks->begin());
      a->deme() = undeme;
      blacks->erase(a);
    }
    delete blacks;
    // erase deme information from nodes.
    for (node_t *p : *this) {
      p->deme() = undeme;
    }
    // drop superfluous nodes (holding just one ball).
    comb();
    ndeme() = 0;
    return *this;
  };

  void curtail (slate_t tnew, slate_t troot);
 
private:
 
  void reuniqify (name_t shift);
 
public:

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

  genealogy_t& insert_zlb (void) {
    for (node_t *p : *this) {
      if (p->holds(green) && p->holds(blue)) {
        assert(!p->holds(black)); // genealogy should have already been pruned
        ball_t *b = p->last_ball();
        assert(b->is(blue));
        node_t *q = make_node(p->deme());
        q->slate = p->slate;
        swap(q->green_ball(),b);
        push_back(q);
      }
    }
    sort();
    return *this;
  };
 
private:
  void clip_zlb (void);
  void cap_tips (void);
  void cap_roots (void) {
    node_nit j = begin();
    while (j != end()) {
      if ((*j)->is_root() && (*j)->slate > timezero()) {
        node_t *q = make_node();
        q->slate = timezero();
        attach(q,*j);
        push_front(q);
      }
      j++;
    }
    sort();
  };

  node_t *scan_branch_label (string_t::const_iterator b, string_t::const_iterator e, node_t *p, slate_t bl);
 
public:

  genealogy_t& parse (const string_t&);
 
  void time_rescale (slate_t scale, slate_t origin = 0) {
    timezero() = scale*(timezero()-origin);
    for (node_t *p : *this ) {
      p->slate = scale*(p->slate-origin);
    }
    time() = scale*(time()-origin);
  };

  friend SEXP cblv (genealogy_t&);

  genealogy_t& parse_cblv (const double *, const double *, int, double);
 
private:

  std::pair<std::vector<slate_t>, std::vector<slate_t>> cblv (void) const;
 
};
 
#endif