| Title: | A Code Generator for ODE-Based Models |
|---|---|
| Description: | Provides an R6 class and several utility methods to facilitate the implementation of models based on ordinary differential equations. The heart of the package is a code generator that creates compiled 'Fortran' (or 'R') code which can be passed to a numerical solver. There is direct support for solvers contained in packages 'deSolve' and 'rootSolve'. |
| Authors: | David Kneis <[email protected]> |
| Maintainer: | David Kneis <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.8.6 |
| Built: | 2025-02-05 05:57:26 UTC |
| Source: | https://github.com/dkneis/rodeo |
This package provides methods to
import a conceptual ODE-based model stored in tabular form (i.e. as text files or spreadsheets).
generate source code (either R or Fortran) to be passed to an ODE-solver.
visualize and export basic information about a model, e.g. for documentation purposes.
Consult the package vignette for details. The concept of writing an ODE system in tabular/matrix form is nicely introduced, e. g., in the book of Reichert, P., Borchardt, D., Henze, M., Rauch, W., Shanahan, P., Somlyody, L., and Vanrolleghem, P. A. (2001): River water quality model No. 1, IWA publishing, ISBN 9781900222822.
The current source code repository is https://github.com/dkneis/rodeo.
See rodeo-class for the rodeo class
and the corresponding class methods.
Type help(package="rodeo") or see the links below to access the
documentation of non-class methods contained in the package.
buildFromWorkbook Builds and compiles a model
fully specified in a workbook (supports .xlsx and .ods format).
forcingFunctions Generation of forcing functions
in Fortran.
exportDF Export of data frames as TEX or HTML code.
stoiCreate Creates a stoichiometry matrix from a set
of chemical reactions.
stoiCheck Validates a stoichiometry matrix by checking
for conservation of mass.
Useful links:
The function builds a rodeo-based model by importing
all declarations and equations from tables stored as delimited text.
buildFromText( declarations, equations, sep = "\t", dim = 1, set_defaults = TRUE, fortran = FALSE, sources = NULL, ... )buildFromText( declarations, equations, sep = "\t", dim = 1, set_defaults = TRUE, fortran = FALSE, sources = NULL, ... )
declarations |
File path of a delimited text file holding the declaration of state variables, parameters, and functions. See below for details about the expected file contents. |
equations |
File path of a delimited text file holding mathematical expressions of process rates and stoichiometric factors forming the right hand sides of a system of simultaneous ODE. See below for details about the expected file contents. |
sep |
The column delimiter used in the input text files. |
dim |
The number of spatial compartments, possibly in multiple
dimensions. For single-box models without spatial resolution, use
|
set_defaults |
If |
fortran |
Controls the language of code generation. The default
( |
sources |
Only relevant if |
... |
Optional arguments passed to |
An object of class rodeo.
The delimited text files provided as input are parsed by
read.table with header=TRUE and the delimiter
specified by sep. The files must contain the following:
'declarations' Declares the identifiers of state variables,
parameters, and functions used in the model equations. Mandatory columns
are 'type', 'name', 'unit', 'description', and 'default'. Entries in
the type column must be one of 'variable', 'parameter', or 'function'.
If source code is generated for R (fortran=FALSE), any declared
functions must be accessible in the environment where the model is
run. If fortran=TRUE, the functions must be implemented in the
file(s) listed in sources to be included in compilation.
Entries in the 'name' column must be unique, valid identifier names
(character followed by characters, digits, underscore).
Entries in the 'default' column shall be numeric.
'equations' Specifies the model equations. Mandatory columns are 'name', 'unit', 'description', 'rate' plus one column for every state variable of the model. The 'rate' columns holds math expressions for the process rates and columns named after state variables contain the corresponding expressions representing stoichiometric factors. All columns are of type character.
The best way to understand the contents of the input files is to study
the examples in the folder 'models' shipped with the package. Type
system.file("models", package="rodeo") at the R prompt to see
where this folder is installed on your system. A full example is given below.
David Kneis [email protected]
buildFromWorkbook provides similar functionality
# Build and run a SEIR type epidemic model decl <- system.file("models/SEIR_declarations.txt", package="rodeo") eqns <- system.file("models/SEIR_equations.txt", package="rodeo") m <- buildFromText(decl, eqns) x <- m$dynamics(times=0:30, fortran=FALSE) print(head(x))# Build and run a SEIR type epidemic model decl <- system.file("models/SEIR_declarations.txt", package="rodeo") eqns <- system.file("models/SEIR_equations.txt", package="rodeo") m <- buildFromText(decl, eqns) x <- m$dynamics(times=0:30, fortran=FALSE) print(head(x))
The function builds a rodeo-based model by importing
all declarations and equations from a workbook established with common
spreadsheet software.
buildFromWorkbook( workbook, dim = 1, set_defaults = TRUE, fortran = FALSE, sources = NULL, ... )buildFromWorkbook( workbook, dim = 1, set_defaults = TRUE, fortran = FALSE, sources = NULL, ... )
workbook |
File path of the workbook with extension '.xlsx'. See below for the mandatory worksheets that must be present in the workbook. |
dim |
The number of spatial compartments, possibly in multiple
dimensions. For single-box models without spatial resolution, use
|
set_defaults |
If |
fortran |
Controls the language of code generation. The default
( |
sources |
Only relevant if |
... |
Optional arguments passed to |
An object of class rodeo.
The file provided as workbook must contain two sheets:
'declarations' Declares the identifiers of state variables,
parameters, and functions used in the model equations. Mandatory columns
are 'type', 'name', 'unit', 'description', and 'default'. Entries in
the type column must be one of 'variable', 'parameter', or 'function'.
If source code is generated for R (fortran=FALSE), any declared
functions must be accessible in the environment where the model is
run. If fortran=TRUE, the functions must be implemented in the
file(s) listed in sources to be included in compilation.
Entries in the 'name' column must be unique, valid identifier names
(character followed by characters, digits, underscore).
Entries in the 'default' column shall be numeric.
'equations' Specifies the model equations. Mandatory columns are 'name', 'unit', 'description', 'rate' plus one column for every state variable of the model. The 'rate' columns holds math expressions for the process rates and columns named after state variables contain the corresponding expressions representing stoichiometric factors. All columns are of type character.
The best way to understand the contents of a suitable workbook is to study
the examples in the folder 'models' shipped with the package. Type
system.file("models", package="rodeo") at the R prompt to see
where this folder is installed on your system. A full example is given below.
David Kneis [email protected]
buildFromText provides similar functionality
# Build and run a SEIR type epidemic model m <- buildFromWorkbook( system.file("models/SEIR.xlsx", package="rodeo") ) x <- m$dynamics(times=0:30, fortran=FALSE) print(head(x))# Build and run a SEIR type epidemic model m <- buildFromWorkbook( system.file("models/SEIR.xlsx", package="rodeo") ) x <- m$dynamics(times=0:30, fortran=FALSE) print(head(x))
Creates and 'compiles' a function for use with numerical methods from
package deSolve or rootSolve.
sources |
Name(s) of source files(s) where functions appearing in
process rates or stoichiometric factors are implemented. Can be |
fortran |
If |
target |
Name of a 'target environment'. Currently, 'deSolve' is the only supported value. |
lib |
File path to be used for the generated library (without the platform specific extension). Note that any uppercase characters will be converted to lowercase. By default, the file is created in R's temporary folder under a random name. |
reuse |
If |
invisible(NULL)
The expected language of the external code passed in sources
depends on the value of fortran.
If fortran is FALSE, R code is generated and made executable
by eval and parse. Auxiliary code
passed via sources is made available via source.
The created R function is stored in the object.
If fortran is TRUE, the external code passed in
sources must implement a module with the fixed name 'functions'.
This module must contain all user-defined functions referenced in process
rates or stoichiometric factors.
If fortran is TRUE, a shared library is created. The library
is immediately loaded with dyn.load and it is
automatically unloaded with dyn.unload when the
object's finalize method is called.
The name of the library (base name without extension) as well as the name
of the function to compute the derivatives are stored in the object.
These names can be queried with the
libName and libFunc methods, respectively.
Unless a file path is specified via the lib argument, the library is
created in the folder returned by tempdir under a
unique random name.
This method internally calls generate.
data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) # This would trigger compilation assuming that 'functionsCode.f95' contains # a Fortran implementation of all functions; see vignette for full example ## Not run: model$compile(sources="functionsCode.f95") ## End(Not run)data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) # This would trigger compilation assuming that 'functionsCode.f95' contains # a Fortran implementation of all functions; see vignette for full example ## Not run: model$compile(sources="functionsCode.f95") ## End(Not run)
Compute a dynamic solution with the numerical algorithms from package
deSolve.
times |
Times of interest (numeric vector). |
fortran |
Switch between compiled Fortran and R code (logical). Default
is |
proNames |
Assign names to output columns holding the process rates?
Default is |
... |
Auxiliary arguments passed to |
The matrix returned by the integrator (see ode).
This method can only be used after compile has been
called.
The ... argument should not be used to assign values to any of
y, parms, times, func. If fortran is
TRUE it should also not assign values to dllname,
nout, or outnames. All these arguments of
ode get their appropriate values automatically.
Use step for integration over a single time step
with a built-in, Fortran-based solver.
Generates code to include tabular data in a tex document or web site.
exportDF( x, tex = FALSE, colnames = NULL, width = NULL, align = NULL, funHead = NULL, funCell = NULL, lines = TRUE, indent = 2 )exportDF( x, tex = FALSE, colnames = NULL, width = NULL, align = NULL, funHead = NULL, funCell = NULL, lines = TRUE, indent = 2 )
x |
The data frame being exported. |
tex |
Logical. Allows to switch between generation of TEX code and HTML. |
colnames |
Displayed column names. If |
width |
Either |
align |
Either |
funHead |
Either |
funCell |
Like |
lines |
Logical. Switches table borders on/off. |
indent |
Integer. Number of blanks used to indent the generated code. |
A character string (usually needs to be exported to a file).
The functions funHead and funCell are useful to apply
formatting or character replacement. For example, one could use
function(x) {paste0("\\bold{",toupper(x),"}")}
to generate bold, uppercase column names in a TEX table.
David Kneis [email protected]
The xtable packages provides similar functionality with
more sophisticated options. Consider the 'pandoc' software do convert
documents from one markup language to another one. Finally, consider the
latex package 'datatools' for direct inclusion of delimited text files
(e.g. produced by write.table) in tex documents.
# Create example table df <- data.frame(stringsAsFactors=FALSE, name= c("growth", "dead"), unit= c("1/d","1/d"), expression= c("r * N * (1 - N/K)"," d * N")) # Export as TEX: header in bold, 1st colum in italics, last column as math tex <- exportDF(df, tex=TRUE, colnames=c(expression="process rate expression"), width=c(expression=0.5), align=c(expression="p"), funHead=setNames(replicate(ncol(df), function(x){paste0("\\textbf{",x,"}")}),names(df)), funCell=c(name=function(x){paste0("\\textit{",x,"}")}, expression=function(x){paste0("$",x,"$")}) ) cat(tex,"\n") # Export as HTML: non-standard colors are used for all columns tf <- tempfile(fileext=".html") write(x= exportDF(df, tex=FALSE, funHead=setNames(replicate(ncol(df), function(x){paste0("<font color='red'>",x,"</font>")}),names(df)), funCell=setNames(replicate(ncol(df), function(x){paste0("<font color='blue'>",x,"</font>")}),names(df)) ), file=tf) ## Not run: browseURL(tf) file.remove(tf) ## End(Not run)# Create example table df <- data.frame(stringsAsFactors=FALSE, name= c("growth", "dead"), unit= c("1/d","1/d"), expression= c("r * N * (1 - N/K)"," d * N")) # Export as TEX: header in bold, 1st colum in italics, last column as math tex <- exportDF(df, tex=TRUE, colnames=c(expression="process rate expression"), width=c(expression=0.5), align=c(expression="p"), funHead=setNames(replicate(ncol(df), function(x){paste0("\\textbf{",x,"}")}),names(df)), funCell=c(name=function(x){paste0("\\textit{",x,"}")}, expression=function(x){paste0("$",x,"$")}) ) cat(tex,"\n") # Export as HTML: non-standard colors are used for all columns tf <- tempfile(fileext=".html") write(x= exportDF(df, tex=FALSE, funHead=setNames(replicate(ncol(df), function(x){paste0("<font color='red'>",x,"</font>")}),names(df)), funCell=setNames(replicate(ncol(df), function(x){paste0("<font color='blue'>",x,"</font>")}),names(df)) ), file=tf) ## Not run: browseURL(tf) file.remove(tf) ## End(Not run)
rodeo ObjectClean-up method for objects of the rodeo-class.
The method is called implicitly for its side effects when a
rodeo object is destroyed.
At present, the method just unloads the object-specific
shared libraries created with the compile or
initStepper methods.
Generates Fortran code to return the current values of forcing functions based on interpolation in tabulated time series data.
forcingFunctions(x)forcingFunctions(x)
x |
Data frame with colums 'name', 'file', 'column', 'mode', 'default'. See below for expected entries. |
A character string holding generated Fortran code. Must be written to
disk, e.g. using write, prior to compilation.
The fields of the input data frame are interpreted as follows:
name Name of the forcing function as declared in the table
of functions.
file Name of the text file containing the time series data
either as an absolute or relative path.
Time information is expected as numeric values in the first column (e.g.
as number of seconds after some reference date).
The period is used as the decimal character in floating point numbers,
numeric values can also be given in scientific format (e.g. as 0.314e+1).
Allowed column delimiters are blank, tab, or comma. A sequence of white
spaces collapses to a single delimiter but this is not the case for
commas. It is strictly recommended to use a consistent delimiter
character within a particular file.
Blank lines are allowed everywhere in the file, comment lines must start
with a '#'. The first non-blank, non-comment line is interpreted as
column headers and the name of the first column (holding time info)
is essentially ignored).
column Name of the column in file from which data are
to be read.
mode Integer code to control how the interpolation is
performed. Use 0 for constant interpolation with full weight given to the
value at the end of a time interval. Use 1 for constant interpolation
with full weight given to the value at the begin of a time interval. Any
other values (< 0 or > 1) result in linear interpolation with weights
being set automatically.
default Logical. If FALSE, the generated function has
the interface 'f(time)'. If TRUE, the generated function has a
two-argument interface 'f(time, z)'. If the actual argument 'z' is
NaN, the function behaves just like the single-argument version,
i.e. interpolation in tabulated data is performed. If 'z' is not
NaN, the function returns the value of 'z'.
The generated code provides a single module named 'forcings' which defines
as many forcing functions as there are rows in x.
The module 'forcings' needs to be made available to the compiler
(either at the command line or via inclusion in another file with Fortran's
include mechanism). In addition, it must be referenced in the module
'functions' with an appropriate 'use' statement (see example below).
The generated function return scalar values of type double precision. If an error condition is encountered, the return value of a functions equals the largest possible double precision value (generated by Fortran's 'huge' function). In addition, errors trigger calls of the subroutines 'rexit' (at top level) or 'rwarn' (at lower levels). These two functions are available automatically if the Fortran code is compiled using 'R CMD SHLIB'. Otherwise, the two functions need to be defined (see examples below).
In the two-argument version, the second argument is tested against NaN
using 'ISNAN'. This function is not part of the Fortran standard but it is
supported by most compilers, including gfortan. The Fortran 2003 standard
conformal function would be 'IS_IEEE_NAN' which is not yet supported by
compiler versions normally installed with R (March 2016).
David Kneis [email protected]
## Not run: ! Example of a Fortran file to define functions include 'forcings.f95' ! include generated forcings file in compilation module functions use forcings ! make forcings available as functions implicit none contains ! ... any non-forcing functions go here ... end module ## End(Not run) ## Not run: ! Definition of 'rexit' and 'rwarn' for testing of the generated code ! outside of R subroutine rexit (x) character(len=*), intent(in):: x write(*,*) "ERROR: ",trim(adjustl(x)) stop 1 end subroutine subroutine rwarn (x) character(len=*), intent(in):: x write(*,*) "WARNING: ",trim(adjustl(x)) end subroutine ## End(Not run)## Not run: ! Example of a Fortran file to define functions include 'forcings.f95' ! include generated forcings file in compilation module functions use forcings ! make forcings available as functions implicit none contains ! ... any non-forcing functions go here ... end module ## End(Not run) ## Not run: ! Definition of 'rexit' and 'rwarn' for testing of the generated code ! outside of R subroutine rexit (x) character(len=*), intent(in):: x write(*,*) "ERROR: ",trim(adjustl(x)) stop 1 end subroutine subroutine rwarn (x) character(len=*), intent(in):: x write(*,*) "WARNING: ",trim(adjustl(x)) end subroutine ## End(Not run)
Declaration of functions referenced at the ODE's right hand sides of the bacteria growth example model.
A data frame with the following fields:
Name of the function
Unit of the return value
Short description (text)
This is a low-level method to translate the ODE-model specification into a
function that computes process rates and the state variables derivatives
(either in R or Fortran). You probably want to use the high-level method
compile instead, which uses generate internally.
lang |
Character string to select the language of the generated source code. Currently either 'f95' (for Fortran) or 'r' (for R). |
name |
Name for the generated function (character string). It should start with a letter, optionally followed by letters, numbers, or underscores. |
The generated source code as a string. Must be written to
disk, e.g. using write, prior to compilation.
Details of this low-level method may change in the future.
See other methods of the rodeo-class, especially
compile which internally uses this method.
data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) fortranCode <- model$generate(lang="f95") ## Not run: write(fortranCode, file="") ## End(Not run)data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) fortranCode <- model$generate(lang="f95") ## Not run: write(fortranCode, file="") ## End(Not run)
Query values of parameters of a rodeo-based model.
asArray |
Logical. If |
useNames |
Logical. Used to enable/disable element names for the return
vector when |
A numeric vector or array.
The corresponding 'set' method is setPars and examples
can be found there. Use getVars to query the values of
variables rather than parameters.
Query values of state variables of a rodeo-based model.
asArray |
Logical. If |
useNames |
Logical. Used to enable/disable element names for the return
vector when |
A numeric vector or array.
The corresponding 'set' method is setVars and examples
can be found there. Use getPars to query the values of
parameters rather than variables.
rodeo ObjectInitializes an object of the rodeo-class with data frames
holding the specification of an ODE system.
vars |
Declaration of state variables appearing in the ODE system. Data frame with mandatory columns 'name', 'unit', 'description'. |
pars |
Declaration of parameters (i.e. constants) appearing in the ODE
system. Data frame with the same mandatory columns as |
funs |
Declaration of functions being referenced in the ODE
system. Data frame with the same mandatory columns as |
pros |
Declaration of process rates. Data frame with mandatory columns 'name', 'unit', 'description', 'expression'. |
stoi |
Declaration of stoichiometric factors. A data frame with
mandatory columns 'variable', 'process', 'expression', if |
asMatrix |
Logical. Specifies whether stoichiometry information is given in matrix or data base format. |
dim |
An integer vector, specifying the number of boxes in each spatial
dimension. Use |
The method is called implicitly for its side effects when a
rodeo object is instantiated with new.
It has no accessible return value.
The mandatory fields of the input data frames should be of type
character. Additional fields may be present in these data frames and the
contents becomes part of the rodeo object.
The 'expression' fields of pros and stoi (or the contents of
the stoichiometry matrix) should be valid mathematical expressions in R and
Fortran. These can involve the names of declared state variables,
parameters, and functions as well as numeric constants or basic math
operators. Branching or loop constructs are not allowed (but these can
appear inside referenced functions).
There are currently few reserved words that cannot be used as variable,
parameter, function, or process names. The reserved words are 'time',
'left', and 'right'.
Initialization does not assign numeric values to state variables or
parameters. Use the decicated methods setVars and
setPars for that purpose.
See the package vignette for examples.
data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) print(model)data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) print(model)
Initializes rodeo's built-in ODE solver. This method must be called
prior to using step. Only works with Fortran, not plain R.
sources |
Name(s) of fortran source file(s) where a module with the
fixed name 'functions' is implemented. This module must contain all
user-defined functions referenced in process rates or
stoichiometric factors. Can be |
method |
Name of a the solver. Currently, 'rk5' is the only supported value (Runge-Kutta method of Cash and Karp). |
invisible(NULL)
After this method was called, step can be used to
perform the integration.
To perform integration with the solvers from package
deSolve use compile instead of this
method.
Return the name of the library function for use with
deSolve or rootSolve methods.
The name of the function to compute derivatives which is contained in
the library built with compile. This name has to be supplied as
the func argument of the solver methods
in deSolve or rootSolve.
Return the pure name of the shared library for use with
deSolve or rootSolve methods.
The base name of the shared library file created with
compile after stripping of the the platform specific
extension. This name has to be supplied as the dllname argument of
the solver methods in deSolve or
rootSolve.
Declaration of parameters of the bacteria growth example model.
A data frame with the following fields:
Name of the parameter
Unit
Short description (text)
Visualizes the stoichiometry matrix using standard plot methods. The sign of stoichiometric factors is displayed as upward and downward pointing triangles. Also visualized are dependencies of process rates on variables.
box |
A positive integer pointing to a spatial sub-unit of the model. |
time |
Time. The value is ignored in the case of autonomous models. |
cex |
Character expansion factor. |
colPositive |
Color for positive stoichiometric factors. |
colNegative |
Color for negative stoichiometric factors. |
colInteract |
Color used to highlight dependencies. |
colBack |
Color of background. |
colGrid |
Color of a grid. |
lwdGrid |
Grid line width. |
translateVars |
Optional function to recode variable labels. Must take the original vector as argument and return the altered version. |
translatePros |
Optional function to recode process labels. Must take the original vector as argument and return the altered version. |
The values of state variables and parameters must have been set using
the setVars and setPars methods. If the
stoichiometric factors are mathematical expressions involving
function references, these functions must be defined in R (even if the
numerical computations are based on generated Fortran code).
See other methods of the rodeo-class or
stoichiometry for computing the stoichiometric factors only.
Alternative options for displaying stoichiometry information are described
in the package vignette.
data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) model$setVars(c(bac=0.1, sub=0.5)) model$setPars(c(mu=0.8, half=0.1, yield= 0.1, vol=1000, flow=50, sub_in=1)) monod <- function(c,h) {c/(c+h)} model$plotStoichiometry(box=c(1))data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) model$setVars(c(bac=0.1, sub=0.5)) model$setPars(c(mu=0.8, half=0.1, yield= 0.1, vol=1000, flow=50, sub_in=1)) monod <- function(c,h) {c/(c+h)} model$plotStoichiometry(box=c(1))
Definition of processes of the bacteria growth example model.
A data frame with the following fields:
Name of the process
Unit of the rate expression
Short description (text)
Mathematical expression (as a string)
rodeo ClassThis documents the rodeo class to represent an ODE-based
model. See the rodeo-package main page or type
help(package="rodeo") for an introduction to the package of
the same name.
prosTblA data frame with fields 'name', 'unit', 'description', and 'expression' defining the process rates.
stoiTblA data frame with fields 'variable', 'process', and 'expression' reprenting the stoichiometry matrix in data base format.
varsTblA data frame with fields 'name', 'unit', 'description' declaring
the state variables of the model. The declared names become valid
identifiers to be used in the expression fields of prosTbl or stoiTbl.
parsTblA data frame of the same structure as vars declaring the
parameters of the model. The declared names become valid
identifiers to be used in the expression fields of prosTbl or stoiTbl.
funsTblA data frame of the same structure as vars declaring any
functions referenced in the expression fields of prosTbl or stoiTbl.
dimInteger vector specifying the spatial dimensions.
varsNumeric vector, holding values of state variables.
parsNumeric vector, holding values of parameters.
For most of the methods below, a separate help page is available describing its arguments and usage.
initialize Initialization method for new objects.
namesVars, namesPars, namesFuns, namesPros Functions
returning the names of declared state variables, parameters,
functions, and processes, respectively (character vectors). No arguments.
lenVars, lenPars, lenFuns, lenPros Functions
returning the number of declared state variables, parameters, functions
and processes, respectively (integer). No arguments.
getVarsTable, getParsTable, getFunsTable, getProsTable,
getStoiTable Functions returning the data frames used to initialize
the object. No arguments
getDim Returns the spatial dimensions as an integer vector.
No arguments.
compile Compiles a Fortran library for use with
numerical methods from packages deSolve or
rootSolve.
generate Translate the ODE-model specification into a
function that computes process rates and the state variables' derivatives
(either in R or Fortran). Consider to use the high-level method
compile.
setVars Assign values to state variables.
setPars Assign values to parameters.
getVars Returns the values of state variables.
getPars Returns the values of parameters.
stoichiometry Returns the stoichiometry matrix, either
evaluated (numeric) or as text.
plotStoichiometry Plots qualitative stoichiometry
information.
See the rodeo-package main page or type
help(package="rodeo") to find the documentation of any non-class
methods contained in the rodeo package.
## Example using high-level functions: Epidemiological SEIR model # Import model from workbook shipped with this package m <- buildFromWorkbook(system.file("models/SEIR.xlsx", package="rodeo")) # Set parameters and initial state (defaults stored in workbook) p <- m$getParsTable() m$setPars(setNames(p$default, p$name)) v <- m$getVarsTable() m$setVars(setNames(v$default, v$name)) # Run a dynamic simulations and print parts of the result sim <- m$dynamics(time=1:30, fortran=FALSE) print(head(sim)) print(tail(sim)) ### Example using low-level functions: Bacterial growth in bioreactor # Creation of model object (data frames shipped as package data) data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) # Assignment of parameters and initial values model$setPars(c(mu=0.8, half=0.1, yield= 0.1, vol=1000, flow=50, sub_in=1)) model$setVars(c(bac=0.01, sub=0)) # Implementation of functions declared in 'funs' monod <- function(c,h) {c/(c+h)} # Creation of derivatives function in a low-level way; calling # the 'compile' method is a more convenient alternative code <- model$generate(name="derivs", lang="r") derivs <- eval(parse(text=code)) # Explicit call of an integrator from the deSolve package times <- 0:96 out <- deSolve::ode(y=model$getVars(), times=times, func=derivs, parms=model$getPars()) colnames(out) <- c("time", model$namesVars(), model$namesPros()) plot(out)## Example using high-level functions: Epidemiological SEIR model # Import model from workbook shipped with this package m <- buildFromWorkbook(system.file("models/SEIR.xlsx", package="rodeo")) # Set parameters and initial state (defaults stored in workbook) p <- m$getParsTable() m$setPars(setNames(p$default, p$name)) v <- m$getVarsTable() m$setVars(setNames(v$default, v$name)) # Run a dynamic simulations and print parts of the result sim <- m$dynamics(time=1:30, fortran=FALSE) print(head(sim)) print(tail(sim)) ### Example using low-level functions: Bacterial growth in bioreactor # Creation of model object (data frames shipped as package data) data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) # Assignment of parameters and initial values model$setPars(c(mu=0.8, half=0.1, yield= 0.1, vol=1000, flow=50, sub_in=1)) model$setVars(c(bac=0.01, sub=0)) # Implementation of functions declared in 'funs' monod <- function(c,h) {c/(c+h)} # Creation of derivatives function in a low-level way; calling # the 'compile' method is a more convenient alternative code <- model$generate(name="derivs", lang="r") derivs <- eval(parse(text=code)) # Explicit call of an integrator from the deSolve package times <- 0:96 out <- deSolve::ode(y=model$getVars(), times=times, func=derivs, parms=model$getPars()) colnames(out) <- c("time", model$namesVars(), model$namesPros()) plot(out)
Triggers dynamic simulation(s) of a 0D model for one to many scenarios. The individual scenarios can differ by the initial state or the values of parameters. Optionally produces basic plots.
times |
Numeric vector defining the times for which the future state of the system is computed. |
scenarios |
Either |
fortran |
Logical. Passed to the respective argument of
|
plot.vars |
Logical. Plot the dynamics of all state variables? |
plot.pros |
Logical. Plot the dynamics of process rates? |
leg |
Keyword to set the position of the legend (if plots are created). |
mar |
Numeric vector of length 4 to control figure margins. See the
|
... |
Possible optional arguments passed to the numerical solver,
namely |
A data frame with at least three columns and one row for each
time requested via the times argument. The first column
indicates the scenario, the second column holds the time, and further
columns hold the values of state variables and process rates.
If the scenarios argument is used to update initial values and /
or parameters, the following applies: For each parameter not being included
(by name) in the vector for a particular scenario, the default value
will be used (as stored in the workbook from which the model was imported).
The same is true for the initial values of variables. See the examples
below.
David Kneis [email protected]
Look at buildFromWorkbook for how to create a
suitable model object.
# build the model m <- buildFromWorkbook(system.file("models/SEIR.xlsx", package="rodeo")) # run with defaults x <- m$scenarios(times=0:30, scenarios=NULL) # run with updated values x <- m$scenarios(times=0:30, scenarios=list(default=c(t=1, i=0.2, r=0.4), fast=c(t=1, i=0.5, r=0.1)) )# build the model m <- buildFromWorkbook(system.file("models/SEIR.xlsx", package="rodeo")) # run with defaults x <- m$scenarios(times=0:30, scenarios=NULL) # run with updated values x <- m$scenarios(times=0:30, scenarios=list(default=c(t=1, i=0.2, r=0.4), fast=c(t=1, i=0.5, r=0.1)) )
Assign values to parameters of a rodeo-based model.
x |
A numeric vector or array, depending on the model's spatial
dimensions. Consult the help page of the sister method
|
NULL (invisible). The assigned numeric data are stored in the
object and can be accessed by the getPars method.
Look at the notes and examples for the setVars method.
The corresponding 'get' method is getPars. Use
setVars to assign values to variables rather than parameters.
Consult the help page of the latter function for examples.
Assign values to state variables of a rodeo-based model.
x |
A numeric vector or array, depending on the model's spatial dimensions. See details below. |
NULL (invisible). The assigned numeric data are stored in the
object and can be accessed by the getVars method.
For a 0-dimensional model (i.e. a model without spatial resolution),
x must be a numeric vector whose length equals the number of state
variables. The element names of x must match those returned by the
object's namesVars method. See the examples for how to bring the
vector elements into required order.
For models with a spatial resolution, x must be a numeric array of
proper dimensions. The last dimension (cycling slowest) corresponds to the
variables and the first dimension (cycling fastest) corresponds to the
models' highest spatial dimension. Thus, dim(x) must be equal to
c(rev(obj$getDim()), obj$namesVars()) where obj is the object
whose variables are to be assigned. The names of the array's last dimension
must match the return value of obj$namesVars().
In the common 1-dimensional case, this just means that x must be a
matrix with column names matching the return value of
obj$namesVars() and as many rows as given by obj$getDim().
The corresponding 'get' method is getVars. Use
setPars to assign values to parameters rather than variables.
data(vars, pars, funs, pros, stoi) x0 <- c(bac=0.1, sub=0.5) # 0-dimensional model model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) model$setVars(x0) print(model$getVars()) # How to sort vector elements x0 <- c(sub=0.5, bac=0.1) # doesn't match order of variables model$setVars(x0[model$namesVars()]) # 1-dimensional model with 3 boxes nBox <- 3 model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(nBox)) x1 <- matrix(rep(x0, each=nBox), nrow=nBox, ncol=model$lenVars(), dimnames=list(NULL, model$namesVars())) model$setVars(x1) print(model$getVars()) print(model$getVars(asArray=TRUE)) # 2-dimensional model with 3 x 4 boxes model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(3,4)) x2 <- array(rep(x0, each=3*4), dim=c(4,3,model$lenVars()), dimnames=list(dim2=NULL, dim1=NULL, variables=model$namesVars())) model$setVars(x2) print(model$getVars()) print(model$getVars(asArray=TRUE))data(vars, pars, funs, pros, stoi) x0 <- c(bac=0.1, sub=0.5) # 0-dimensional model model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) model$setVars(x0) print(model$getVars()) # How to sort vector elements x0 <- c(sub=0.5, bac=0.1) # doesn't match order of variables model$setVars(x0[model$namesVars()]) # 1-dimensional model with 3 boxes nBox <- 3 model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(nBox)) x1 <- matrix(rep(x0, each=nBox), nrow=nBox, ncol=model$lenVars(), dimnames=list(NULL, model$namesVars())) model$setVars(x1) print(model$getVars()) print(model$getVars(asArray=TRUE)) # 2-dimensional model with 3 x 4 boxes model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(3,4)) x2 <- array(rep(x0, each=3*4), dim=c(4,3,model$lenVars()), dimnames=list(dim2=NULL, dim1=NULL, variables=model$namesVars())) model$setVars(x2) print(model$getVars()) print(model$getVars(asArray=TRUE))
Performs integration over a single time step using a built-in ODE solver. At present, a single solver is implement with limited options. The interface of this method may change when support for other solvers is added.
t0 |
Numeric. Initial time. |
h |
Numeric. Length of time step of interest. |
hmin |
Minimum tolerated internal step size. The default of |
maxsteps |
Maximum tolerated number of sub-steps. |
tol |
Numeric. Relative accuracy requested. This is currently a global value, i.e. one cannot set the accuracy per state variable. |
method |
String. Currently, 'rk5' is the only method implemented. This is a Runge-Kutta Cash-Karp solver adapted from Press et al. (2002), Numerical recipes in Fortran 90. It is designed to handle non-stiff problems only. |
check |
Logical. Can be used to avoid repeated checks of arguments. This may increase performance in repeated calls. |
A named numeric vector holding the values of state variables and process rates in all boxes.
This method can only be used after a call to initStepper
has been made.
Use deSolve for advanced solvers with more
options and capabilities to handle stiff problems.
Definition of the links between simulated processes and state variables in the bacteria growth example model.
A data frame with the following fields:
Name of the state variable
Name of the process
Mathematical expression (as a string)
Validates the stoichiometry matrix by checking for conservation of mass (more typically conservation of moles).
stoiCheck(stoi, comp, env = globalenv(), zero = .Machine$double.eps * 2)stoiCheck(stoi, comp, env = globalenv(), zero = .Machine$double.eps * 2)
stoi |
Stoichiometry matrix either in evaluated
( |
comp |
Matrix defining the elemental composition of compounds.
Column names of |
env |
An environment or list supplying constants, functions, and
operators needed to evaluate expressions in |
zero |
A number close to zero. If the absolute result value of a mass balance computation is less than this, the result is set to 0 (exactly). |
A numeric matrix with the following properties:
There is one row for each process, thus the number and names of rows
are the same as in stoi.
There is one column per checked element, hence column names are
equal to the row names of comp.
A matrix element at position represent the mass balance
for the process in row with respect to the element in column
. A value of zero indicates a close balance. Positive (negative)
values indicate that mass is gained (lost) in the respective process.
David Kneis [email protected]
Use stoiCreate to create a stoichiometry matrix
from a set of reactions in common notation.
# Eq. 1 and 2 are from Soetaert et al. (1996), Geochimica et Cosmochimica # Acta, 60 (6), 1019-1040. 'OM' is organic matter. Constants 'nc' and 'pc' # represent the nitrogen/carbon and phosphorus/carbon ratio, respectively. reactions <- c( oxicDegrad= "OM + O2 -> CO2 + nc * NH3 + pc * H3PO4 + H2O", denitrific= "OM + 0.8*HNO3 -> CO2 + nc*NH3 + 0.4*N2 + pc*H3PO4 + 1.4*H2O", dissPhosp1= "H3PO4 <-> H + H2PO4", dissPhosp2= "H2PO4 <-> H + HPO4" ) # Non-evaluated stoichiometry matrix stoi <- stoiCreate(reactions, toRight="_f", toLeft="_b") # Parameters ('nc' and 'pc' according to Redfield ratio) pars <- list(nc=16/106, pc=1/106) # Elemental composition comp <- rbind( OM= c(C=1, N="nc", P="pc", H="2 + 3*nc + 3*pc"), O2= c(C=0, N=0, P=0, H=0), CO2= c(C=1, N=0, P=0, H=0), NH3= c(C=0, N=1, P=0, H=3), H3PO4= c(C=0, N=0, P=1, H=3), H2PO4= c(C=0, N=0, P=1, H=2), HPO4= c(C=0, N=0, P=1, H=1), H2O= c(C=0, N=0, P=0, H=2), HNO3= c(C=0, N=1, P=0, H=1), N2= c(C=0, N=2, P=0, H=0), H= c(C=0, N=0, P=0, H=1) ) # We need the transposed form comp <- t(comp) # Mass balance check bal <- stoiCheck(stoi, comp=comp, env=pars) print(bal) print(colSums(bal) == 0)# Eq. 1 and 2 are from Soetaert et al. (1996), Geochimica et Cosmochimica # Acta, 60 (6), 1019-1040. 'OM' is organic matter. Constants 'nc' and 'pc' # represent the nitrogen/carbon and phosphorus/carbon ratio, respectively. reactions <- c( oxicDegrad= "OM + O2 -> CO2 + nc * NH3 + pc * H3PO4 + H2O", denitrific= "OM + 0.8*HNO3 -> CO2 + nc*NH3 + 0.4*N2 + pc*H3PO4 + 1.4*H2O", dissPhosp1= "H3PO4 <-> H + H2PO4", dissPhosp2= "H2PO4 <-> H + HPO4" ) # Non-evaluated stoichiometry matrix stoi <- stoiCreate(reactions, toRight="_f", toLeft="_b") # Parameters ('nc' and 'pc' according to Redfield ratio) pars <- list(nc=16/106, pc=1/106) # Elemental composition comp <- rbind( OM= c(C=1, N="nc", P="pc", H="2 + 3*nc + 3*pc"), O2= c(C=0, N=0, P=0, H=0), CO2= c(C=1, N=0, P=0, H=0), NH3= c(C=0, N=1, P=0, H=3), H3PO4= c(C=0, N=0, P=1, H=3), H2PO4= c(C=0, N=0, P=1, H=2), HPO4= c(C=0, N=0, P=1, H=1), H2O= c(C=0, N=0, P=0, H=2), HNO3= c(C=0, N=1, P=0, H=1), N2= c(C=0, N=2, P=0, H=0), H= c(C=0, N=0, P=0, H=1) ) # We need the transposed form comp <- t(comp) # Mass balance check bal <- stoiCheck(stoi, comp=comp, env=pars) print(bal) print(colSums(bal) == 0)
Return and optionally evaluate the mathematical expression appearing in the stoichiometry matrix.
box |
Either |
time |
Time. The value is ignored in the case of autonomous models. |
A matrix of numeric or character type, depending on the value of
box.
If the stoichiometric factors are mathematical expressions involving function references, these functions must be defined in R (even if the numerical computations are based on generated Fortran code).
See other methods of the rodeo-class or
plotStoichiometry for a graphical representation of the
stoichiometric factors only.
data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) model$setPars(c(mu=0.8, half=0.1, yield= 0.1, vol=1000, flow=50, sub_in=1)) model$setVars(c(bac=0.1, sub=0.5)) monod <- function(c,h) {c/(c+h)} print(model$stoichiometry(box=NULL)) print(model$stoichiometry(box=c(1)))data(vars, pars, funs, pros, stoi) model <- rodeo$new(vars, pars, funs, pros, stoi, dim=c(1)) model$setPars(c(mu=0.8, half=0.1, yield= 0.1, vol=1000, flow=50, sub_in=1)) model$setVars(c(bac=0.1, sub=0.5)) monod <- function(c,h) {c/(c+h)} print(model$stoichiometry(box=NULL)) print(model$stoichiometry(box=c(1)))
Creates a stoichiometry matrix from a set of reaction equations.
stoiCreate( reactions, eval = FALSE, env = globalenv(), toRight = "_forward", toLeft = "_backward" )stoiCreate( reactions, eval = FALSE, env = globalenv(), toRight = "_forward", toLeft = "_backward" )
reactions |
A named vector of character strings, each representing a (chemical) reaction. See syntax details below. |
eval |
Logical. If |
env |
Only relevant if |
toRight |
Only relevant for reversible reactions. The passed character
string is appended to the name of the respective element of
|
toLeft |
Like |
A matrix with the following properties:
The number of columns equals the total number of components present
in reactions. The components' names are used as column names.
The number of rows equals the length of reactions plus the
number of reversible reactions. Thus, a single row is created for each
non-reversible reaction but two rows are created for reversible ones.
The latter represent the forward and backward reaction (in that order).
The row names are constructed from the names of reactions, making
use of the suffixes toRight and toLeft in the case of
reversible reactions.
The matrix is filled with the stoichiometric factors extracted from
reactions. Empty elements are set to zero.
The type of the matrix (character or
numeric) depends on the value of eval.
The syntax rules for reaction equations are as follows (see examples):
There must be a left hand side and a right hand side. Sides must be separated by one of the arrows '->', '<-', or '<->' with the latter indicating a reversible reaction.
Names of component(s) must appear at each side of the reaction. These must be legal row/column names in R. If multiple components are consumed or produced, they must be separated by '+'.
Any stoichiometric factors need to appear before the respective component name using '*' as the separating character. Stoichiometric factors being equal to unity can be omitted.
A stoichiometric factor is treated as a mathematical expression. In common cases, it is just a numeric constant. However, the expression can also involve references to variables or functions. If such an expression contains math operators '*' or '+' it needs to be enclosed in parenthesis.
David Kneis [email protected]
Use stoiCheck to validate the mass balance of
the generated matrix.
# EXAMPLE 1: From https://en.wikipedia.org/wiki/Petersen_matrix (July 2016) # reactions <- c( formS= "A + 2 * B -> S", equiES= "E + S <-> ES", decoES= "ES -> E + P" ) stoi <- stoiCreate(reactions, eval=TRUE, toRight="_f", toLeft="_b") print(stoi) # EXAMPLE 2: Decomposition of organic matter (selected equations only) # # Eq. 1 and 2 are from Soetaert et al. (1996), Geochimica et Cosmochimica # Acta, 60 (6), 1019-1040. 'OM' is organic matter. Constants 'nc' and 'pc' # represent the nitrogen/carbon and phosphorus/carbon ratio, respectively. reactions <- c( oxicDegrad= "OM + O2 -> CO2 + nc * NH3 + pc * H3PO4 + H2O", denitrific= "OM + 0.8*HNO3 -> CO2 + nc*NH3 + 0.4*N2 + pc*H3PO4 + 1.4*H2O", dissPhosp1= "H3PO4 <-> H + H2PO4", dissPhosp2= "H2PO4 <-> H + HPO4" ) # Non-evaluated matrix stoi <- stoiCreate(reactions, toRight="_f", toLeft="_b") print(stoi) # Evaluated matrix ('nc' and 'pc' according to Redfield ratio) pars <- list(nc=16/106, pc=1/106) stoi <- stoiCreate(reactions, eval=TRUE, env=pars, toRight="_f", toLeft="_b") print(stoi)# EXAMPLE 1: From https://en.wikipedia.org/wiki/Petersen_matrix (July 2016) # reactions <- c( formS= "A + 2 * B -> S", equiES= "E + S <-> ES", decoES= "ES -> E + P" ) stoi <- stoiCreate(reactions, eval=TRUE, toRight="_f", toLeft="_b") print(stoi) # EXAMPLE 2: Decomposition of organic matter (selected equations only) # # Eq. 1 and 2 are from Soetaert et al. (1996), Geochimica et Cosmochimica # Acta, 60 (6), 1019-1040. 'OM' is organic matter. Constants 'nc' and 'pc' # represent the nitrogen/carbon and phosphorus/carbon ratio, respectively. reactions <- c( oxicDegrad= "OM + O2 -> CO2 + nc * NH3 + pc * H3PO4 + H2O", denitrific= "OM + 0.8*HNO3 -> CO2 + nc*NH3 + 0.4*N2 + pc*H3PO4 + 1.4*H2O", dissPhosp1= "H3PO4 <-> H + H2PO4", dissPhosp2= "H2PO4 <-> H + HPO4" ) # Non-evaluated matrix stoi <- stoiCreate(reactions, toRight="_f", toLeft="_b") print(stoi) # Evaluated matrix ('nc' and 'pc' according to Redfield ratio) pars <- list(nc=16/106, pc=1/106) stoi <- stoiCreate(reactions, eval=TRUE, env=pars, toRight="_f", toLeft="_b") print(stoi)
Declaration of variables of the bacteria growth example model.
A data frame with the following fields:
Name of the variable
Unit
Short description (text)