Skip to content

Commit

Permalink
Modify the R and Rcpp code generation to support the R library Simula…
Browse files Browse the repository at this point in the history 
…tionModels. The new versions will perform automatic conversion of outputs and inputs, so that the user is working in the same units as defined in the model description but the calculations are done using only reference units (i.e. no scaling prefixes).
  • Loading branch information
@AleMorales
AleMorales committed Feb 25, 2015
1 parent f611da5 commit 7fc4d97
Showing 1 changed file with 123 additions and 105 deletions.
228 changes: 123 additions & 105 deletions src/codegeneration.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -283,20 +283,22 @@ function create_return_line_R(states, observed)
return_line = "return(list(c("
for i in states
if i != states[end]
return_line = return_line * i * ","
return_line *= i * ","
else
return_line = return_line * i
return_line *= i
end
end
return_line = return_line * "), c("
return_line *= ")"
length(observed) > 0 && (return_line *= ", c(")
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 * ")))"
length(observed) > 0 && (return_line *= ")")
return_line = return_line * "))"
end

# Create the function in R on the Equations section of the model Dict
Expand Down Expand Up @@ -325,7 +327,7 @@ end

function generate_code_R!(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[]
observed = ASCIIString[]
for (key,val) in ode_model.Equations
val.Exported && push!(observed, key)
end
Expand All @@ -334,6 +336,10 @@ function generate_code_R!(ode_model::OdeSource; unit_analysis = true, name = "au
names_derivatives[i] = "d_"*names_derivatives[i]*"_dt"
end
deleteat!(observed, findin(observed, names_derivatives))
coef_observed = OrderedDict{ASCIIString, Float64}()
for i in observed
coef_observed[i] = ode_model.Equations[i].Units.f
end

# Sort the equations
sorted_model = sort_equations(ode_model)
Expand Down Expand Up @@ -362,74 +368,73 @@ function generate_code_R!(ode_model::OdeSource; unit_analysis = true, name = "au
model_function = create_function_R!(sorted_model,observed)

# Create the default arguments
named_states = OrderedDict{String, Any}()
named_states = OrderedDict{ASCIIString, Float64}()
coef_states = OrderedDict{ASCIIString, Float64}()
for (key,val) in sorted_model.States
named_states[key] = val.Value * val.Units.f
coef_states[key] = val.Units.f
end
named_parameters = OrderedDict{String, Any}()
named_parameters = OrderedDict{ASCIIString, Float64}()
coef_parameters = OrderedDict{ASCIIString, Float64}()
for (key,val) in sorted_model.Parameters
named_parameters[key] = val.Value * val.Units.f
coef_parameters[key] = val.Units.f
end
forcings = OrderedDict{String,Any}()
forcings = OrderedDict{ASCIIString,Any}()
coef_forcings = OrderedDict{ASCIIString,Float64}()
c = 0
for (key,value) in sorted_model.Forcings
c += 1
forcings[key] = (float(value.Time), float(value.Value)*value.Units.f)
coef_forcings[key] = value.Units.f
end
write_model_R!(named_states,named_parameters, forcings, observed,
write_model_R!(named_states,coef_states, named_parameters, coef_parameters,
forcings, coef_forcings, observed,coef_observed,
model_function, name, file)
nothing
end


function write_model_R!(States::OrderedDict{String, Any},
Parameters::OrderedDict{String, Any},
Forcings::OrderedDict{String, Any},
Observed::Array{String, 1},
Model::String,
name::String,
file::String)
function write_model_R!(States::OrderedDict{ASCIIString, Float64},
Coef_states::OrderedDict{ASCIIString, Float64},
Parameters::OrderedDict{ASCIIString, Float64},
Coef_parameters::OrderedDict{ASCIIString, Float64},
Forcings::OrderedDict{ASCIIString, Any},
Coef_forcings::OrderedDict{ASCIIString, Float64},
Observed::Array{ASCIIString, 1},
Coef_observed::OrderedDict{ASCIIString, Float64},
Model::ASCIIString,
name::ASCIIString,
file::ASCIIString)
f = open("$(file).R","w")
println(f, "library(simecol); library(RcppSundials)")
println(f, "$name <- new(\"odeModel\",")
println(f, "main = $Model, ")
println(f, "library(SimulationModels); library(RcppSundials)")
println(f, "$name <- ODEmodel\$new(")
transformed_states = string(States)[2:(end-1)]
transformed_states = replace(transformed_states, "=>", "=")
println(f, "init = c($transformed_states),")
units = replace(string(Coef_states)[2:(end-1)], "=>", "=")
println(f, "States = list(Values = c($transformed_states), Coefs = c($units)),")
transformed_parameters = string(Parameters)[2:(end-1)]
transformed_parameters = replace(transformed_parameters, "=>", "=")
println(f, "parms = c($transformed_parameters),")
println(f, "inputs = list(")
units = replace(string(Coef_parameters)[2:(end-1)], "=>", "=")
println(f, "Parameters = list(Values = c($transformed_parameters), Coefs = c($units)),")
println(f, "Forcings = list(Values = 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)]
units = replace(string(Coef_forcings)[2:(end-1)], "=>", "=")
println(f, "), Coefs = c($units)),")
println(f, "Time = 0:1,")
names_observed = string(Observed)[13:(end-1)]
units = replace(string(Coef_observed)[2:(end-1)], "=>", "=")
println(f, "Observed = list(names = c($names_observed), Coefs = c($units)),")
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 = func, jacobian = function(t, states, parameters, forcings) {0}))
colnames(sim) = c(\"time\", $names_states, $names_observed)
class(sim) = c(\"deSolve\",\"matrix\")
sim
})
test = sim($name)
Settings = list(rtol = 1e-6,atol = 1e-10, maxsteps = 1e5, maxord = 5, hini = 0,
hmin = 0, hmax = 0, maxerr = 12, maxnonlin = 12,
maxconvfail = 12, method = "bdf", maxtime = 0, jacobian = 0),
""")
println(f, "model = $Model)")
close(f)
nothing
end
Expand Down Expand Up @@ -469,15 +474,37 @@ function create_return_line_Rcpp(states, observed)
return_line = return_line * "};\n array<vector<double>,2> output{derivatives, observed};\n return output;"
end

# Substitute the annoying x^y for pow(x,y)
function substitute_power(ex::Expr)
new_ex = deepcopy(ex)
for i in 1:length(new_ex.args)
if isa(new_ex.args[i], Expr)
new_ex.args[i] = substitute_power(new_ex.args[i])
elseif new_ex.args[i] == :^
new_ex.args[i] = :pow
end
end
return new_ex
end

substitute_power(ex::Any) = ex

# 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"
if(ismatch(r"params|states|forcs", string(rhs.Expr)))
number = match(r"(?<=\[)[\d]+(?=\])",string(rhs.Expr))
new_expr = replace(string(rhs.Expr), "[$(number.match)]", "[$(number.match)-1]")
code *= "const double $(new_expr);\n"
else
code *= "const double $(rhs.Expr);\n"
end
else
code *= "const double $lhs" * " = " * replace(string(rhs.Expr), ".*", "*") * ";\n"
mod_expr = substitute_power(rhs.Expr)
code *= "const double $lhs" * " = " * replace(string(mod_expr), ".*", "*") * ";\n"
end
end
end
Expand All @@ -492,35 +519,27 @@ function create_function_Rcpp!(model::OdeSorted, observed, name)
code = replace(code, "1.0 *", "")
up_boiler_plate =
"""
#define ARMA_DONT_USE_CXX11
#include <RcppArmadillo.h>
#include <Rcpp.h>
#include <array>
#include <vector>
#include <math.h>
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
const vector<double>& params, const vector<double>& forcs) { \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[]
observed = ASCIIString[]
for (key,val) in ode_model.Equations
val.Exported && push!(observed, key)
end
Expand All @@ -529,6 +548,10 @@ function generate_code_Rcpp!(ode_model::OdeSource; unit_analysis = true, name =
names_derivatives[i] = "d_"*names_derivatives[i]*"_dt"
end
deleteat!(observed, findin(observed, names_derivatives))
coef_observed = OrderedDict{ASCIIString, Float64}()
for i in observed
coef_observed[i] = ode_model.Equations[i].Units.f
end

# Sort the equations
sorted_model = sort_equations(ode_model)
Expand Down Expand Up @@ -557,92 +580,87 @@ function generate_code_Rcpp!(ode_model::OdeSource; unit_analysis = true, name =
model_function = create_function_Rcpp!(sorted_model,observed, name)

# Create the default arguments
named_states = OrderedDict{String, Any}()
named_states = OrderedDict{ASCIIString, Float64}()
coef_states = OrderedDict{ASCIIString, Float64}()
for (key,val) in sorted_model.States
named_states[key] = val.Value * val.Units.f
coef_states[key] = val.Units.f
end
named_parameters = OrderedDict{String, Any}()
named_parameters = OrderedDict{ASCIIString, Float64}()
coef_parameters = OrderedDict{ASCIIString, Float64}()
for (key,val) in sorted_model.Parameters
named_parameters[key] = val.Value * val.Units.f
coef_parameters[key] = val.Units.f
end
forcings = OrderedDict{String,Any}()
forcings = OrderedDict{ASCIIString,Any}()
coef_forcings = OrderedDict{ASCIIString,Float64}()
c = 0
for (key,value) in sorted_model.Forcings
c += 1
forcings[key] = (float(value.Time), float(value.Value)*value.Units.f)
coef_forcings[key] = value.Units.f
end
write_model_Rcpp!(named_states,named_parameters, forcings, observed,
write_model_Rcpp!(named_states,coef_states, named_parameters, coef_parameters,
forcings, coef_forcings, observed,coef_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)
function write_model_Rcpp!(States::OrderedDict{ASCIIString, Float64},
Coef_states::OrderedDict{ASCIIString, Float64},
Parameters::OrderedDict{ASCIIString, Float64},
Coef_parameters::OrderedDict{ASCIIString, Float64},
Forcings::OrderedDict{ASCIIString, Any},
Coef_forcings::OrderedDict{ASCIIString, Float64},
Observed::Array{ASCIIString, 1},
Coef_observed::OrderedDict{ASCIIString, Float64},
Model::ASCIIString,
name::ASCIIString,
file::ASCIIString)
# 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)
"""library(SimulationModels); 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
$(name)_pointer = getNativeSymbolInfo(name = "$name",PACKAGE = "$file")\$address
""")
println(f, "$name <- new(\"odeModel\",")
println(f, "main = function() {NULL}")
println(f, "$name <- ODEmodel\$new(")
transformed_states = string(States)[2:(end-1)]
transformed_states = replace(transformed_states, "=>", "=")
println(f, "init = c($transformed_states),")
units = replace(string(Coef_states)[2:(end-1)], "=>", "=")
println(f, "States = list(Values = c($transformed_states), Coefs = c($units)),")
transformed_parameters = string(Parameters)[2:(end-1)]
transformed_parameters = replace(transformed_parameters, "=>", "=")
println(f, "parms = c($transformed_parameters),")
println(f, "inputs = list(")
units = replace(string(Coef_parameters)[2:(end-1)], "=>", "=")
println(f, "Parameters = list(Values = c($transformed_parameters), Coefs = c($units)),")
println(f, "Forcings = list(Values = 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)
""")
units = replace(string(Coef_forcings)[2:(end-1)], "=>", "=")
println(f, "), Coefs = c($units)),")
println(f, "Time = 0:1,")
names_observed = string(Observed)[13:(end-1)]
units = replace(string(Coef_observed)[2:(end-1)], "=>", "=")
println(f, "Observed = list(names = c($names_observed), Coefs = c($units)),")
println(f, """
Settings = list(rtol = 1e-6,atol = 1e-10, maxsteps = 1e5, maxord = 5, hini = 0,
hmin = 0, hmax = 0, maxerr = 12, maxnonlin = 12,
maxconvfail = 12, method = "bdf", maxtime = 0, jacobian = 0),
""")
println(f, "model = $(name)_pointer)")
close(f)
nothing
end
Expand Down

0 comments on commit 7fc4d97

Please sign in to comment.