diff --git a/R/RcppExports.R b/R/RcppExports.R
index 0ffcae2d1a8a9ac27a9a4c51de4c725631b1a273..4365d9f9ac378e96d6e0d583178d9e84e1a8c471 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -1,7 +1,7 @@
# 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(`_GSSA_CPPSSA`, K, tmax, Nmax, No, nrolls, MES, b0, d0, intervalTime, result_path, result_path_details, print_debug)
+CPPSSA <- function(K, tmax, Nmax, No, nrolls, n0, MES, b, d, v, u, a, yI, y0, crit, intervalTime, result_path = "", result_path_details = "", print_debug = FALSE) {
+ .Call(`_GSSA_CPPSSA`, K, tmax, Nmax, No, nrolls, n0, MES, b, d, v, u, a, yI, y0, crit, intervalTime, result_path, result_path_details, print_debug)
}
diff --git a/R/main.R b/R/main.R
index 860f087e91cc8a2c1f9e382518ce031da5c9a27a..ba7278be63c17a86436efc51473b3c2277e0bd21 100644
--- a/R/main.R
+++ b/R/main.R
@@ -4,13 +4,20 @@
#' @param Nmax A Integer :
#' @param No A Integer :
#' @param nrolls A Integer :
+#' @param n0 A Vector :
#' @param ES A Matrix :
-#' @param b0 A Numeric :
-#' @param d0 A Numeric :
+#' @param b A Vector :
+#' @param d A Vector :
+#' @param v A Matrix :
+#' @param u A Matrix :
+#' @param a A Matrix :
+#' @param yI A Vector :
+#' @param y0 A Vector :
+#' @param crit A Function :
#' @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)
+RSSA <- function(K, tmax, Nmax, No, nrolls, n0, ES, b, d, v, u, a, yI, y0, crit = NULL, intervalTime = 0.1, result_path = "", result_path_details = "", print_debug = FALSE) {
+ CPPSSA(K, tmax, Nmax, No, nrolls, n0, ES, b, d, v, u, a, yI, y0, crit, intervalTime, result_path, result_path_details, print_debug)
}
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 718993ee9d294dedb0d384dc4665b3e7ffdd3f21..6a6b721073cf90bac514ccde58f3c644c9953c22 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -6,8 +6,8 @@
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 _GSSA_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) {
+NumericMatrix CPPSSA(const unsigned long K, const unsigned long tmax, const unsigned long Nmax, const unsigned long No, unsigned long nrolls, NumericVector n0, NumericMatrix MES, NumericVector b, NumericVector d, NumericMatrix v, NumericMatrix u, NumericMatrix a, NumericVector yI, NumericVector y0, Nullable<Function> crit, double intervalTime, const std::string result_path, const std::string result_path_details, bool print_debug);
+RcppExport SEXP _GSSA_CPPSSA(SEXP KSEXP, SEXP tmaxSEXP, SEXP NmaxSEXP, SEXP NoSEXP, SEXP nrollsSEXP, SEXP n0SEXP, SEXP MESSEXP, SEXP bSEXP, SEXP dSEXP, SEXP vSEXP, SEXP uSEXP, SEXP aSEXP, SEXP yISEXP, SEXP y0SEXP, SEXP critSEXP, SEXP intervalTimeSEXP, SEXP result_pathSEXP, SEXP result_path_detailsSEXP, SEXP print_debugSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
@@ -16,20 +16,27 @@ BEGIN_RCPP
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< NumericVector >::type n0(n0SEXP);
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< NumericVector >::type b(bSEXP);
+ Rcpp::traits::input_parameter< NumericVector >::type d(dSEXP);
+ Rcpp::traits::input_parameter< NumericMatrix >::type v(vSEXP);
+ Rcpp::traits::input_parameter< NumericMatrix >::type u(uSEXP);
+ Rcpp::traits::input_parameter< NumericMatrix >::type a(aSEXP);
+ Rcpp::traits::input_parameter< NumericVector >::type yI(yISEXP);
+ Rcpp::traits::input_parameter< NumericVector >::type y0(y0SEXP);
+ Rcpp::traits::input_parameter< Nullable<Function> >::type crit(critSEXP);
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));
+ rcpp_result_gen = Rcpp::wrap(CPPSSA(K, tmax, Nmax, No, nrolls, n0, MES, b, d, v, u, a, yI, y0, crit, intervalTime, result_path, result_path_details, print_debug));
return rcpp_result_gen;
END_RCPP
}
static const R_CallMethodDef CallEntries[] = {
- {"_GSSA_CPPSSA", (DL_FUNC) &_GSSA_CPPSSA, 12},
+ {"_GSSA_CPPSSA", (DL_FUNC) &_GSSA_CPPSSA, 19},
{NULL, NULL, 0}
};
diff --git a/src/main.cpp b/src/main.cpp
index f733b00e7d58f240b832d33b5122533d23ba9c52..ba42ba442e237ea0f5c3f78352c285d4970b5981 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -11,9 +11,9 @@
#include <Rcpp.h>
using namespace Rcpp;
-double get_value(const std::vector<double>& v, unsigned long i, unsigned long j, unsigned long k) {
+/*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;
@@ -23,7 +23,6 @@ unsigned long get_number_pop(const std::vector<unsigned long>& npop) {
return n;
}
-
class Types {
public:
@@ -161,9 +160,9 @@ std::vector<unsigned long> ceffect(unsigned long K, Types& e) {
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) {
+double crate(unsigned long K, Types& e, std::vector<unsigned long>& npop, NumericVector& b, NumericVector& d,
+ NumericMatrix& v, NumericMatrix& u, NumericMatrix& a,
+ NumericVector& yI, NumericVector& y0) {
Types e2 = e.sign();
@@ -172,11 +171,11 @@ double crate(unsigned long K, Types& e, std::vector<unsigned long>& npop, std::v
} 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];
+ return u( (-e.e1)-1, e.e2-1) * 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];
+ return v( e.e1-1, (-e.e2)-1) * 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 a( (-e.e1)-1, (-e.e2)-1) * npop[(-e.e1)-1] * npop[(-e.e2)-1];
}
return 0.0;
@@ -186,8 +185,9 @@ double crate(unsigned long K, Types& e, std::vector<unsigned long>& npop, std::v
// [[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) {
+ unsigned long nrolls, NumericVector n0, NumericMatrix MES, NumericVector b, NumericVector d, NumericMatrix v,
+ NumericMatrix u, NumericMatrix a, NumericVector yI, NumericVector y0, Nullable<Function> crit, 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;
@@ -196,16 +196,22 @@ NumericMatrix CPPSSA(const unsigned long K, const unsigned long tmax, const unsi
std::random_device rd;
std::mt19937 gen(rd());
- std::vector<unsigned long> n0(K, static_cast<unsigned long>(No / K));
+ std::vector<unsigned long> cppn0;
+ for(unsigned long i = 0; i < n0.size(); ++i) {
+ cppn0.push_back(n0[i]);
+ }
+
+ //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<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);
@@ -238,7 +244,7 @@ NumericMatrix CPPSSA(const unsigned long K, const unsigned long tmax, const unsi
unsigned long Nt = 0;
for (short i : check) {
- Nt += n0[i - 1];
+ Nt += cppn0[i - 1];
}
if(print_debug) {
@@ -254,14 +260,15 @@ NumericMatrix CPPSSA(const unsigned long K, const unsigned long tmax, const unsi
}
std::vector<unsigned long>npop;
- std::copy(n0.begin(), n0.end(), std::back_inserter(npop));
+
+ std::copy(cppn0.begin(), cppn0.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);
+ Population p(0.0, cppn0);
tNs.push_back(p);
std::vector<double> rs(ES.size(), 0);
@@ -272,7 +279,26 @@ NumericMatrix CPPSSA(const unsigned long K, const unsigned long tmax, const unsi
double limit = intervalLimit;
- while((t < tmax) && (0 < Nt) && (Nt < Nmax)) {
+
+ while( true ) {
+ //while((t < tmax) && (0 < Nt) && (Nt < Nmax)) {
+
+ bool stop = true;
+
+ if(crit.isNotNull()) {
+ Function f = crit.get();
+ stop = Rcpp::as<bool>(f(t, tmax, Nmax, npop));
+ } else {
+ if((t < tmax) && (0 < Nt) && (Nt < Nmax)) {
+ stop = true;
+ } else {
+ stop = false;
+ }
+ }
+
+ if( ! stop ) {
+ break;
+ }
if(it < nt) {
it++;
@@ -367,30 +393,3 @@ NumericMatrix CPPSSA(const unsigned long K, const unsigned long tmax, const unsi
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;
-}
-