Add a first prototype for C++ version of code. This would only work i…

…nside an R package (cannot be loaded via R CMD SHLIB). Also, exponents are not corrected.

src/codegeneration.jl

@@ -411,7 +411,11 @@ function write_model_R!(States::OrderedDict{String, Any},
  names_observed = string(Observed)[8:(end-1)]
  names_states = string(collect(keys(States)))[8:(end-1)]
  println(f, """
-solver = function(y, times, func, parms, inputs) {
+solver = function(y, times, func, parms, inputs, ...) {
+settings = list(...)
+if(!is.null(settings["atol"])) atol = settings["atol"] else atol = 1e-6
+if(!is.null(settings["rtol"])) rtol = settings["rtol"] else rtol = 1e-6
+if(!is.null(settings["method"])) method = settings["method"] else method = "bdf"
 sim = as.matrix(cvode_R(times = times, states = y,
  parameters = parms,
  forcings_data = inputs,
@@ -437,3 +441,216 @@ function generate_code_R!(source::String; unit_analysis = false,name = "autogene
  ode_model = convert_reaction_model(reaction_model)
  generate_code_R!(ode_model, unit_analysis = unit_analysis, name = name, file = file, jacobian = jacobian, sensitivities = sensitivities)
 end
+
+####################################################################################
+####################################################################################
+############################ RC++ CODE GENERATION ################################
+####################################################################################
+####################################################################################
+
+# Create a return line with an array containing the derivatives and observed functions (for the STL version of the model)
+function create_return_line_Rcpp(states, observed)
+ return_line = "vector<double> derivatives{"
+ for i in states
+ if i != states[end]
+ return_line = return_line * i * ","
+ else
+ return_line = return_line * i
+ end
+ end
+ return_line = return_line * "};\n vector<double> observed{"
+ for i in observed
+ if i != observed[end]
+ return_line = return_line * i * ", "
+ else
+ return_line = return_line * i
+ end
+ end
+ return_line = return_line * "};\n array<vector<double>,2> output{derivatives, observed};\n return output;"
+end
+
+# Create the function in C++ using STL with the RcppSundials API
+function create_function_Rcpp!(model::OdeSorted, observed, name)
+ code = ""
+ for level in 1:length(model.SortedEquations)
+ for (lhs, rhs) in model.SortedEquations[level]
+ if level == 1
+ code *= "const double $(rhs.Expr);\n"
+ else
+ code *= "const double $lhs" * " = " * replace(string(rhs.Expr), ".*", "*") * ";\n"
+ end
+ end
+ end
+ # Determine what the time derivatives are
+ time_derivatives = String[]
+ for i in collect(keys(model.States))
+ push!(time_derivatives, "d_"*i*"_dt")
+ end
+ # Return line
+ return_line = create_return_line_Rcpp(time_derivatives,observed)
+ # Return the output
+ code = replace(code, "1.0 *", "")
+ up_boiler_plate =
+"""
+#define ARMA_DONT_USE_CXX11
+#include <RcppArmadillo.h>
+#include <Rcpp.h>
+#include <array>
+#include <vector>
+using namespace std;
+using namespace Rcpp;
+
+extern "C" {
+ array<vector<double>, 2> $(name)(const double& t, const vector<double>& states,
+ const vector<double>& parameters, const vector<double>& forcings) { \n
+"""
+ low_boiler_plate =
+"""
+\n
+}
+arma::mat $(name)_jacobian(const double& t, const vector<double>& states,
+ const vector<double>& parameters, const vector<double>& forcings) {
+ arma::mat output = arma::eye(states.size(), states.size());
+ return output;
+}
+};
+"""
+ return paste("\n",up_boiler_plate, code, return_line,low_boiler_plate)
+end
+
+function generate_code_Rcpp!(ode_model::OdeSource; unit_analysis = true, name = "autogenerated_model", file = "autogenerated_model", jacobian = false, sensitivities = false)
+ # Generate the observed variables (everything that is exported but it is not a time derivative)
+ observed = String[]
+ for (key,val) in ode_model.Equations
+ val.Exported && push!(observed, key)
+ end
+ names_derivatives = collect(keys(ode_model.States))
+ for i in 1:length(names_derivatives)
+ names_derivatives[i] = "d_"*names_derivatives[i]*"_dt"
+ end
+ deleteat!(observed, findin(observed, names_derivatives))
+
+ # Sort the equations
+ sorted_model = sort_equations(ode_model)
+
+ # Created compressed model (only if Jacobian or Sensitivities are required!)
+ jacobian_function = "() -> ()"
+ sensitivity_function = "() -> ()"
+ sensitivity_jacobian_function = "() -> ()"
+ if jacobian || sensitivities
+ compressed_model = compress_model(sorted_model, level = 2)
+ jacobian && (jacobian_function = generate_jacobian_function(compressed_model, name))
+ sensitivities && ((sensitivity_function, sensitivity_jacobian_function) = generate_extended_system(compressed_model, name))
+ end
+
+ # Check the units
+ unit_analysis && check_units(sorted_model)
+
+ # Go through the equations and substitute * by ×
+ for i in 1:length(sorted_model.SortedEquations)
+ for (key,val) in sorted_model.SortedEquations[i]
+ sorted_model.SortedEquations[i][key].Expr = sub_product(sorted_model.SortedEquations[i][key].Expr)
+ end
+ end
+
+ # Generate the rhs function
+ model_function = create_function_Rcpp!(sorted_model,observed, name)
+
+ # Create the default arguments
+ named_states = OrderedDict{String, Any}()
+ for (key,val) in sorted_model.States
+ named_states[key] = val.Value * val.Units.f
+ end
+ named_parameters = OrderedDict{String, Any}()
+ for (key,val) in sorted_model.Parameters
+ named_parameters[key] = val.Value * val.Units.f
+ end
+ forcings = OrderedDict{String,Any}()
+ c = 0
+ for (key,value) in sorted_model.Forcings
+ c += 1
+ forcings[key] = (float(value.Time), float(value.Value)*value.Units.f)
+ end
+ write_model_Rcpp!(named_states,named_parameters, forcings, observed,
+ model_function, name, file)
+ nothing
+end
+
+
+function write_model_Rcpp!(States::OrderedDict{String, Any},
+ Parameters::OrderedDict{String, Any},
+ Forcings::OrderedDict{String, Any},
+ Observed::Array{String, 1},
+ Model::String,
+ name::String,
+ file::String)
+ # Create the C++ file
+ f = open("$(file).cpp","w")
+ println(f, Model)
+ close(f)
+ f = open("$(file).R","w")
+ println(f,
+"""library(simecol); library(RcppSundials)
+
+ # These are the pointers to the C++ functions
+ # Note that in the current version the Jacobian is a dummy function (no content)
+ system("R CMD SHLIB $(file).cpp -o $(file).so")
+ dyn.load("$(file).so")
+ # so simulations must always be run with the option Jacobian 0
+ $(name)_pointer = getNativeSymbolInfo(name = "$name",
+ PACKAGE = "$file")\$address
+ $(name)jacobian_pointer = getNativeSymbolInfo(name = "$(name)_jacobian",
+ PACKAGE = "$file")\$address
+
+""")
+ println(f, "$name <- new(\"odeModel\",")
+ println(f, "main = function() {NULL}")
+ transformed_states = string(States)[2:(end-1)]
+ transformed_states = replace(transformed_states, "=>", "=")
+ println(f, "init = c($transformed_states),")
+ transformed_parameters = string(Parameters)[2:(end-1)]
+ transformed_parameters = replace(transformed_parameters, "=>", "=")
+ println(f, "parms = c($transformed_parameters),")
+ println(f, "inputs = list(")
+ forcs = ""
+ for (key,val) in Forcings
+ forcs *= "$key = cbind(c($(string(val[1])[2:(end-1)])),c($(string(val[2])[2:(end-1)]))),"
+ end
+ forcs = forcs[1:(end-1)]
+ println(f, forcs)
+ println(f, "),")
+ println(f, "times = 0:1,")
+ names_observed = string(Observed)[8:(end-1)]
+ names_states = string(collect(keys(States)))[8:(end-1)]
+ println(f,
+"""
+solver = function(y, times, func, parms, inputs, ...) {
+settings = list(...)
+if(!is.null(settings["atol"])) atol = settings["atol"] else atol = 1e-6
+if(!is.null(settings["rtol"])) rtol = settings["rtol"] else rtol = 1e-6
+if(!is.null(settings["method"])) method = settings["method"] else method = "bdf"
+sim = as.matrix(cvode_R(times = times, states = y,
+ parameters = parms,
+ forcings_data = inputs,
+ settings = list(rtol = 1e-6,
+ atol = 1e-6, maxsteps = 1e3, maxord = 5, hini = 1e-3,
+ hmin = 0, hmax = 100, maxerr = 5, maxnonlin = 10,
+ maxconvfail = 10, method = "bdf", maxtime = 1, jacobian = 0),
+ model = $(name)_pointer, jacobian = $(name)jacobian_pointer))
+colnames(sim) = c(\"time\", $names_states, $names_observed)
+class(sim) = c(\"deSolve\",\"matrix\")
+sim
+})
+test = sim($name)
+""")
+ close(f)
+ nothing
+end
+
+
+function generate_code_Rcpp!(source::String; unit_analysis = false,name = "autogenerated_model", file = "autogenerated_model", jacobian = false, sensitivities = false)
+ parsed_model = process_file(source)
+ reaction_model = convert_master_equation(parsed_model)
+ ode_model = convert_reaction_model(reaction_model)
+ generate_code_Rcpp!(ode_model, unit_analysis = unit_analysis, name = name, file = file, jacobian = jacobian, sensitivities = sensitivities)
+end