popul_proc.h

Go to the documentation of this file.

// -*- C++ -*-
// POPUL_PROC class
 
#ifndef _POPUL_PROC_H_
#define _POPUL_PROC_H_
 
#include <string>
#include <cstring>
 
#include "internal.h"


template <class STATE, class PARAMETERS, size_t NEVENT, size_t NDEME = 1>
class popul_proc_t  {
 
protected:
  // TYPE DEFINITIONS
  typedef STATE state_t;
  typedef PARAMETERS parameters_t;
  static const size_t nevent = NEVENT;
  static const size_t ndeme = NDEME;
 
protected:
  // MEMBER DATA
  slate_t next;                 // time of next event
  size_t event;                 // mark of next event
  slate_t current;              // current time
  state_t state;                // current state
  parameters_t params;          // model parameters
 
private:
 
  void clean (void) {};         // memory cleanup
 
public:
 
  // SERIALIZATION
  size_t bytesize (void) const {
    return 2*sizeof(slate_t) + sizeof(size_t)
      + sizeof(state_t) + sizeof(parameters_t);
  };
  size_t bytesize (void) const  {…}

  friend raw_t* operator>> (const popul_proc_t &X, raw_t *o) {
    slate_t A[2]; A[0] = X.current; A[1] = X.next;
    memcpy(o,A,sizeof(A)); o += sizeof(A);
    memcpy(o,&X.event,sizeof(size_t)); o += sizeof(size_t);
    memcpy(o,&X.state,sizeof(state_t)); o += sizeof(state_t);
    memcpy(o,&X.params,sizeof(parameters_t)); o += sizeof(parameters_t);
    return o;
  };
  friend raw_t* operator>> (const popul_proc_t &X, raw_t *o) {…}

  friend raw_t* operator>> (raw_t *o, popul_proc_t &X) {
    X.clean();
    slate_t A[2];
    memcpy(A,o,sizeof(A)); o += sizeof(A);
    X.current = A[0]; X.next = A[1];
    memcpy(&X.event,o,sizeof(size_t)); o += sizeof(size_t);
    memcpy(&X.state,o,sizeof(state_t)); o += sizeof(state_t);
    memcpy(&X.params,o,sizeof(parameters_t)); o += sizeof(parameters_t);
    return o;
  };
  friend raw_t* operator>> (raw_t *o, popul_proc_t &X) {…}
 
public:
  // CONSTRUCTORS, ETC.
  popul_proc_t (double t0 = 0) {
    clean();
    next = current = slate_t(t0);
    event = 0;
  };
  popul_proc_t (double t0 = 0) {…}

  popul_proc_t (raw_t *o) {
    o >> *this;
  };
  popul_proc_t (raw_t *o) {…}

  popul_proc_t (const popul_proc_t & X) {
    raw_t *o = new raw_t[X.bytesize()];
    X >> o;
    o >> *this;
    delete[] o;
  };
  popul_proc_t (const popul_proc_t & X) {…}

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

  popul_proc_t (popul_proc_t &&) = delete;
  popul_proc_t & operator= (popul_proc_t &&) = delete;
  ~popul_proc_t (void) {
    clean();
  };
  ~popul_proc_t (void) {…}
 
public:
 
  // INFORMATION EXTRACTORS
  slate_t time (void) const {
    return current;
  };
  slate_t time (void) const  {…}
 
public:
 
  virtual void valid (void) const {};
 
public:

  void update_params (double*, int);
  void update_IVPs (double*, int);
  double event_rates (double *rate, int n) const;
  virtual void rinit (void) = 0;
  virtual void jump (int e) = 0;
  std::string yaml (std::string tab) const;
 
public:

  void update_clocks (void) {
    double rate[nevent];
    double total_rate = event_rates(rate,nevent);
    if (total_rate > 0) {
      next = current+rexp(1/total_rate);
    } else {
      next = R_PosInf;
    }
    double u = runif(0,total_rate);
    event = 0;
    while (u > rate[event] && event < nevent) {
      if (rate[event] < 0)
        err("in '%s' (%s line %d): invalid negative rate[%zd]=%lg", // #nocov
            __func__,__FILE__,__LINE__,event,rate[event]); // #nocov
      u -= rate[event];
      event++;
    }
    assert(event < nevent);
  };
  void update_clocks (void) {…}

  int play (double tfin) {
    int count = 0;
    if (current > tfin)
      err("cannot simulate backward! (current t=%lg, requested t=%lg)",current,tfin);
    while (next < tfin) {
      current = next;
      jump(event);
      update_clocks();
      count++;
      R_CheckUserInterrupt();
    }
    if (next > tfin) current = tfin; // relies on Markov property
    return count;
  };
  int play (double tfin) {…}
 
};
class popul_proc_t {…};
 
#define PARAM_SET(X) if (!ISNA(p[m])) params.X = p[m];  \
  m++;
#define PARAM_SET(X) if (!ISNA(p[m])) params.X = p[m];  \ …
#define RATE_CALC(X) total += rate[m++] = (X);
#define YAML_PARAM(X) (t + #X + ": " + std::to_string(params.X) + "\n")
#define YAML_STATE(X) (t + #X + ": " + std::to_string(state.X) + "\n")
 
#endif