Add source

.gitignore

+#RStudio
+*.Rproj
+*.Rbuildignore
+.Rproj.user
+.Rhistory
+.RData
+#C++
+# Prerequisites
+*.d
+
+# Compiled Object files
+*.slo
+*.lo
+*.o
+*.obj
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Compiled Dynamic libraries
+*.so
+*.dylib
+*.dll
+
+# Fortran module files
+*.mod
+*.smod
+
+# Compiled Static libraries
+*.lai
+*.la
+*.a
+*.lib
+
+# Executables
+*.exe
+*.out
+*.app

DESCRIPTION

+Package: GillepsieSSA
+Type: Package
+Title: What the Package Does (Title Case)
+Version: 0.1.0
+Author: Who wrote it
+Maintainer: The package maintainer <yourself@somewhere.net>
+Description: More about what it does (maybe more than one line)
+    Use four spaces when indenting paragraphs within the Description.
+License: What license is it under?
+Encoding: UTF-8
+LazyData: true
+Imports: Rcpp (>= 1.0.0)
+LinkingTo: Rcpp
\ No newline at end of file

NAMESPACE

+useDynLib(GillepsieSSA)
+exportPattern("^[^\\CPP]")
+importFrom(Rcpp, evalCpp)

R/RcppExports.R

+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+CPPSSA <- function(K, tmax, Nmax, No, nrolls, MES, b0, d0, intervalTime, result_path = "", result_path_details = "", print_debug = FALSE) {
+    .Call('_GillepsieSSA_CPPSSA', PACKAGE = 'GillepsieSSA', K, tmax, Nmax, No, nrolls, MES, b0, d0, intervalTime, result_path, result_path_details, print_debug)
+}
+

R/main.R

+#' Gillepsie SSA
+#' @param K A Integer :
+#' @param tmax A Integer :
+#' @param Nmax A Integer :
+#' @param No A Integer :
+#' @param nrolls A Integer :
+#' @param ES A Matrix :
+#' @param b0 A Numeric :
+#' @param d0 A Numeric :
+#' @param intervalTime A Numeric : The interval time for the result
+#' @param result_path A Integer : The path of results file (optional)
+#' @param result_path_details A Integer : The path of results details files (optional)
+#' @param print_debug A Logical : more verbose
+RSSA <- function(K, tmax, Nmax, No, nrolls, ES, b0, d0, intervalTime = 0.1, result_path = "", result_path_details = "", print_debug = FALSE) {
+  CPPSSA(K, tmax, Nmax, No, nrolls, ES, b0, d0, intervalTime, result_path, result_path_details, print_debug)
+}

src/Makevars

+PKG_CXXFLAGS= -DCOMPATIBILITYRCPP
+## This is a C++11 package
+CXX_STD = CXX11
+#PKG_CPPFLAGS=-I$(LIB_GSL)/include -I../inst/include
+#PKG_LIBS=-L$(LIB_GSL)/lib -lgsl -lgslcblas $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()")
+

src/Makevars.win

+# -fno-omit-frame-pointer to overcome a (compiler?) bug, in case one would recompile
+# with -O3 (or -UNDEBUG?) for windows 32bit:
+# see https://github.com/RcppCore/RcppEigen/issues/34
+PKG_CXXFLAGS= -DCOMPATIBILITYRCPP -fno-omit-frame-pointer
+## This is a C++11 package
+CXX_STD = CXX11
+#PKG_CPPFLAGS=-I$(LIB_GSL)/include -I../inst/include
+#PKG_LIBS=-L$(LIB_GSL)/lib -lgsl -lgslcblas $(shell "${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "Rcpp:::LdFlags()")
+

src/RcppExports.cpp

+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+// CPPSSA
+NumericMatrix CPPSSA(const unsigned long K, const unsigned long tmax, const unsigned long Nmax, const unsigned long No, unsigned long nrolls, NumericMatrix MES, double b0, double d0, double intervalTime, const std::string result_path, const std::string result_path_details, bool print_debug);
+RcppExport SEXP _GillepsieSSA_CPPSSA(SEXP KSEXP, SEXP tmaxSEXP, SEXP NmaxSEXP, SEXP NoSEXP, SEXP nrollsSEXP, SEXP MESSEXP, SEXP b0SEXP, SEXP d0SEXP, SEXP intervalTimeSEXP, SEXP result_pathSEXP, SEXP result_path_detailsSEXP, SEXP print_debugSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const unsigned long >::type K(KSEXP);
+    Rcpp::traits::input_parameter< const unsigned long >::type tmax(tmaxSEXP);
+    Rcpp::traits::input_parameter< const unsigned long >::type Nmax(NmaxSEXP);
+    Rcpp::traits::input_parameter< const unsigned long >::type No(NoSEXP);
+    Rcpp::traits::input_parameter< unsigned long >::type nrolls(nrollsSEXP);
+    Rcpp::traits::input_parameter< NumericMatrix >::type MES(MESSEXP);
+    Rcpp::traits::input_parameter< double >::type b0(b0SEXP);
+    Rcpp::traits::input_parameter< double >::type d0(d0SEXP);
+    Rcpp::traits::input_parameter< double >::type intervalTime(intervalTimeSEXP);
+    Rcpp::traits::input_parameter< const std::string >::type result_path(result_pathSEXP);
+    Rcpp::traits::input_parameter< const std::string >::type result_path_details(result_path_detailsSEXP);
+    Rcpp::traits::input_parameter< bool >::type print_debug(print_debugSEXP);
+    rcpp_result_gen = Rcpp::wrap(CPPSSA(K, tmax, Nmax, No, nrolls, MES, b0, d0, intervalTime, result_path, result_path_details, print_debug));
+    return rcpp_result_gen;
+END_RCPP
+}
+
+static const R_CallMethodDef CallEntries[] = {
+    {"_GillepsieSSA_CPPSSA", (DL_FUNC) &_GillepsieSSA_CPPSSA, 12},
+    {NULL, NULL, 0}
+};
+
+RcppExport void R_init_GillepsieSSA(DllInfo *dll) {
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+    R_useDynamicSymbols(dll, FALSE);
+}

src/main.cpp

+#include <iostream>
+#include <vector>
+#include <math.h>
+#include <numeric>
+#include <random>
+#include <fstream>
+#include <iterator>
+#include <string>
+#include <algorithm>
+
+#include <Rcpp.h>
+using namespace Rcpp;
+
+double get_value(const std::vector<double>& v, unsigned long i, unsigned long j, unsigned long k) {
+  return v[i*k+j];
+}
+
+unsigned long get_number_pop(const std::vector<unsigned long>& npop) {
+  unsigned long n = 0;
+  for (short i : npop) {
+    n += i;
+  }
+  return n;
+}
+
+
+class Types {
+
+public:
+
+  Types() : e1(0), e2(0), e3(0) {}
+
+  Types(short e1, short e2, short e3) : e1(e1), e2(e2), e3(e3) {}
+
+  void print() {
+    std::cout << "e{" << e1 << ", " << e2 << ", " << e3 << "}" << std::endl;
+  }
+
+  Types sign() {
+    Types t(0, 0, 0);
+
+    if(this->e1 >= 1) {
+      t.e1 = 1;
+    } else if(this->e1 <= -1) {
+      t.e1 = -1;
+    }
+
+    if(this->e2 >= 1) {
+      t.e2 = 1;
+    } else if(this->e2 <= -1) {
+      t.e2 = -1;
+    }
+
+    if(this->e3 >= 1) {
+      t.e3 = 1;
+    } else if(this->e3 <= -1) {
+      t.e3 = -1;
+    }
+
+    return t;
+  }
+
+  short int e1 = 0;
+  short int e2 = 0;
+  short int e3 = 0;
+
+  bool operator==(const Types& b){
+    if(this->e1 != b.e1){
+      return false;
+    }
+
+    if(this->e2 != b.e2){
+      return false;
+    }
+
+    return this->e3 == b.e3;
+  }
+
+};
+
+Types B ( 1,  0, 0);
+Types D ( -1,  0, 0);
+Types T (-1,  1, 0);
+Types I ( 1, -1, 0);
+Types F (-1, -1, 1);
+
+class Population {
+
+public:
+
+  Population() : time(0.0) {}
+  Population(double time, const std::vector<unsigned long>& n0) : time(time) {
+    std::copy(n0.begin(), n0.end(), std::back_inserter(npop));
+  }
+
+  double time;
+  std::vector<unsigned long> npop;
+
+  unsigned long size() {
+    return get_number_pop(this->npop);
+  }
+
+  std::string generate() const {
+
+    unsigned long index = 0;
+    std::string result = "";
+    for (unsigned long val : npop) {
+      if(index == 0) {
+        result += std::to_string(val);
+      } else {
+        result += ";" + std::to_string(val);
+      }
+      index++;
+    }
+    return result;
+  }
+
+};
+
+std::vector<unsigned long> ceffect(unsigned long K, Types& e) {
+  Types e2 = e.sign();
+
+  std::vector<unsigned long> v(K, 0);
+
+  if(e2 == B) {
+    v[e.e1 - 1] = 1;
+  } else if(e2 == D) {
+    v[(-e.e1)-1] = -1;
+  } else if(e2 == T) {
+    v[e.e2-1] = 1;
+
+    std::vector<unsigned long> v2(K, 0);
+    v2[(-e.e1)-1] = 1;
+
+    for (unsigned long i = 0; i < K; i++) {
+      v[i] = v[i] - v2[i];
+    }
+  } else if(e2 == I) {
+    v[e.e1-1] = 1;
+
+    std::vector<unsigned long> v2(K, 0);
+    v2[(-e.e2)-1] = 1;
+
+    for (unsigned long i = 0; i < K; i++) {
+      v[i] = v[i] - v2[i];
+    }
+  } else if(e2 == F) {
+    v[e.e3-1] = 1;
+
+    std::vector<unsigned long> v2(K, 0);
+    v2[(-e.e1)-1] = 1;
+
+    std::vector<unsigned long> v3(K, 0);
+    v3[(-e.e2)-1] = 1;
+
+    for (unsigned long i = 0; i < K; i++) {
+      v[i] = v[i] - v2[i] - v3[i];
+    }
+  }
+
+  return v;
+}
+
+double crate(unsigned long K, Types& e, std::vector<unsigned long>& npop, std::vector<double>& b, std::vector<double>& d,
+             std::vector<double>& v, std::vector<double>& u, std::vector<double>& a,
+             std::vector<unsigned long>& yI, std::vector<unsigned long>& y0) {
+
+  Types e2 = e.sign();
+
+  if(e2 == B) {
+    return b[e.e1-1] * npop[e.e1-1] + yI[e.e1-1];
+  } else if(e2 == D) {
+    return d[(-e.e1)-1] * npop[(-e.e1)-1] + y0[(-e.e1)-1];
+  } else if(e2 == T) {
+    return get_value(u, (-e.e1)-1, e.e2-1, K) * npop[(-e.e1)-1];
+  } else if(e2 == I) {
+    return get_value(v, e.e1-1, (-e.e2)-1, K) * npop[e.e1-1] * npop[(-e.e2)-1];
+  } else if(e2 == F) {
+    return get_value(a, (-e.e1)-1,(-e.e2)-1, K) * npop[(-e.e1)-1] * npop[(-e.e2)-1];
+  }
+
+  return 0.0;
+}
+
+// std::vector<Population>
+
+// [[Rcpp::export]]
+NumericMatrix CPPSSA(const unsigned long K, const unsigned long tmax, const unsigned long Nmax, const unsigned long No,
+                               unsigned long nrolls, NumericMatrix MES, double b0, double d0, double intervalTime, const std::string result_path = "", const std::string result_path_details = "",
+                               bool print_debug = false) {
+
+  if(print_debug) {
+    std::cout << "Start simulation" << std::endl;
+  }
+
+  std::random_device rd;
+  std::mt19937 gen(rd());
+
+  std::vector<unsigned long> n0(K, static_cast<unsigned long>(No / K));
+
+  std::vector<double> b(K, b0);
+  std::vector<double> d(K, d0);
+  std::vector<double> v(K*K, 0);
+  std::vector<double> u(K*K, 0);
+  std::vector<double> a(K*K, 0);
+
+  std::vector<unsigned long> yI(K, 0);
+  std::vector<unsigned long> y0(K, 0);
+
+  std::vector<unsigned long> check(K, 0);
+  std::iota(check.begin(), check.end(), 1);
+
+  std::vector<Types> ES(MES.nrow() * MES.ncol());
+
+  for (unsigned long i = 0; i < MES.nrow(); i++) {
+    Types t(MES(i, 0), MES(i, 1), MES(i, 2));
+    ES[i] = t;
+  }
+
+  /*
+  unsigned long rows_ES = 2*K;
+
+  std::vector<Types> ES(rows_ES);
+
+  for (unsigned long i = 0; i < K; i++) {
+    Types t1(static_cast<unsigned long>(i+1), 0, 0);
+    Types t2(static_cast<unsigned long>(-(i+1)), 0, 0);
+    ES[i] = t1;
+    ES[i+K] = t2;
+  }*/
+
+  std::uniform_real_distribution<double> dist (0, ES.size());
+
+  double t = 0.0;
+
+  unsigned long it = 0;
+
+  unsigned long Nt = 0;
+
+  for (short i : check) {
+    Nt += n0[i - 1];
+  }
+
+  if(print_debug) {
+    std::cout << "t = " << static_cast<int>(t) << "s -> npop = " << Nt << std::endl;
+  }
+
+  std::vector<double> ts(nrolls, 0.0);
+
+  std::exponential_distribution<double> distribution(1);
+
+  for (unsigned long i=0; i<nrolls; ++i) {
+    ts[i] =  distribution(gen);
+  }
+
+  std::vector<unsigned long>npop;
+  std::copy(n0.begin(), n0.end(), std::back_inserter(npop));
+
+  unsigned long Ne = ES.size();
+
+  unsigned long nt = ts.size();
+
+  std::vector<Population> tNs(1);
+  Population p(0.0, n0);
+  tNs.push_back(p);
+
+  std::vector<double> rs(ES.size(), 0);
+
+  unsigned long ie = 0;
+
+  double intervalLimit = 1;
+  double limit = intervalLimit;
+
+
+  while((t < tmax) && (0 < Nt < Nmax)) {
+
+    if(it < nt) {
+      it++;
+    } else {
+      it = 0;
+    }
+
+    for(unsigned long i = 0; i < ES.size(); i++) {
+      rs[i] = crate(K, ES[i], npop, b, d, v, u, a, yI, y0);
+    }
+
+    std::discrete_distribution<> disc(rs.begin(), rs.end());
+
+    ie = static_cast<unsigned long>(disc(gen));
+
+    std::vector<unsigned long> veffect = ceffect(K, ES[ie]);
+
+    for(unsigned long i = 0; i < veffect.size(); i++) {
+      npop[i] += veffect[i];
+    }
+
+    double sum = 0.0;
+
+    for (double val : rs) {
+      sum += val;
+    }
+
+    t = t + ts[it] / sum;
+
+    Nt = 0;
+    for(unsigned long i = 0; i < npop.size(); i++) {
+      Nt += npop[i];
+    }
+
+    if(t > limit) {
+      if(print_debug) {
+        std::cout << "t = " << static_cast<unsigned long>(t) << "s -> npop = " << Nt << std::endl;
+      }
+      limit += intervalLimit;
+    }
+
+    Population cp(t, npop);
+    tNs.push_back(cp);
+  }
+
+  std::ofstream output_file(result_path);
+  std::ofstream output_details_file(result_path_details);
+
+  double currentTime = 0.0;
+  unsigned long currentNumber = 0;
+  double sumVal = 0.0;
+
+  output_file << 0 << ";" << No << std::endl;
+  output_details_file << 0 << ";" << tNs[1].generate() << std::endl;
+
+  std::vector<Population> tmp_tNs;
+
+  for (unsigned long i=0; i<tNs.size(); ++i) {
+    Population p = tNs[i];
+    sumVal += p.size();
+    if(p.time > currentTime ) {
+      tmp_tNs.push_back(p);
+      currentTime += intervalTime;
+      output_file << currentTime << ";" << static_cast<int>(sumVal / currentNumber) << std::endl;
+      output_details_file << currentTime << ";" << p.generate() << std::endl;
+      currentNumber = 0;
+      sumVal = 0;
+
+    }
+    currentNumber++;
+  }
+
+  if(print_debug) {
+    std::cout << "End simulation" << std::endl;
+  }
+
+  NumericMatrix xx (tmp_tNs.size(), K);
+
+  for (int j=0; j<K; j++) {
+    xx(0, j) = tNs[1].npop[j];
+  }
+
+  for (int i=1; i<tmp_tNs.size(); i++) {
+    for (int j=0; j<K; j++) {
+      xx(i, j) = tmp_tNs[i].npop[j];
+      //std::cout << tmp_tNs[i].generate() << std::endl;
+    }
+
+  }
+
+  return xx;
+}
+
+int main() {
+
+  unsigned long K = 10;
+
+  double b0 = 0.5;
+
+  double d0 = 0.45;
+
+  unsigned long No = static_cast<unsigned long>(pow(10, 3));
+
+  unsigned long tmax = 50;
+
+  unsigned long Nmax = 50 * No;
+
+  unsigned long nrolls = static_cast<unsigned long>(pow(10, 6));
+
+  double intervalTime = 0.1;
+
+  bool print_debug = true;
+
+  //std::vector<Population> out1 = CPPSSA(K, tmax, Nmax, No, nrolls, b0, d0, intervalTime, "/home/jimmy/jimmy/project/Gillespie_SSA/output/out1.txt", "/home/jimmy/jimmy/project/Gillespie_SSA/output/out1_details.txt", print_debug);
+
+  //std::vector<Population> out2 = CPPSSA(K, tmax, Nmax, No, nrolls, b0, d0, intervalTime, "/home/jimmy/jimmy/project/Gillespie_SSA/output/out2.txt", "/home/jimmy/jimmy/project/Gillespie_SSA/output/out2_details.txt", print_debug);
+
+  return 0;
+}
+