| solve {paropt} | R Documentation |
Solves an ode-system
Description
Solves an ode equation and calculate an error based on the difference on a user-defined state-data.frame.
Usage
solve(
ode,
parameter,
reltol,
abstol,
states,
solvertype,
own_error_fct,
own_spline_fct,
own_jac_fct,
verbose
)
Arguments
ode |
the ode-system for which the parameter should be optimized. |
parameter |
a data.frame containing the parameters. |
reltol |
a number defining the relative tolerance used by the ode-solver. The default value is 1e-06 |
abstol |
a vector containing the absolute tolerance(s) for each state used by the ode-solver. The default value is 1e-08 |
states |
a data.frame containing the predetermined course of the states. |
solvertype |
a string defining the type of solver which should be used "bdf" or "adams" are the possible values. The default value is "bdf". |
own_error_fct |
An optional function to calculate the error between in silico value and the specified value in the data.frame states. The default error calculation is specified in the section notes. |
own_spline_fct |
An optional function to interpolate the values for variable parameters. The default function is a CatmullRome spline interpolation. |
own_jac_fct |
An optional function which returns the jacobian function. Furthermore it is possible to calculate the jacobian using the R package dfdr. If this is desired "dfdr" has to be passed as argument. If nothing is passed the jacobian matrix is numerically calculated. |
verbose |
A logical value defining whether the output during compilation should be shown or not. The default value is FALSE |
Details
The ode system:
___________________________________________________________
The ode system is an R function which accepts four arguments and returns one.
the first argument is t which defines the (time-) point of then independent variable at which the ode-system is evaluated.
the second argument is a vector called y which defines the current states at timepoint t
the third argument is a vector called ydot which should be filled with the derivative (left hand side) of the ode-system. It has already the correct length! This vector has to be returned.
the last argument is called parameter and is a vector containing the current parameter-set which is tested by the optimization algorithm.
If the parameters can change over time. The already interpolated value is passed to the ode-system.
# theoretical Example: The parameter 'a' can change over time whereas 'b' is constant over time.
parameter_set <- data.frame(
time = c(0, 10, 20, 30, 40),
a = c(1, 2, 3, 4, 5),
b = c(1, NA, NA, NA, NA, NA))
t <- 5
# Interpolation would result in 1.5 for parameter 'a'
parameter <- c(1.5, 1) # 'a', 'b'
y <- 1
ydot <- vector(length(1))
ode(t, y, ydot, parameter)
The function returns ydot. It is only necessary to fill the vector ydot. Check the package ast2ast for more details how this works.
The R function is translated to a C++ function using the package ast2ast, see also ast2ast on CRAN and
ast2ast::translate(). Therefore, if you are calling 'solve' for the first time the function has to be compiled. This can require a bit of time.
The parameters:
___________________________________________________________
The lower and upper boundaries are defined as a data.frame that contains 'time' as the first column.
The subsequent columns contain the information of the parameter.
# Here some examples
# all parameters are constant over the entire integration_time
example1 <- data.frame(
time = 0,
a = 0.4,
b = 1.1,
c = 0.1,
d = 0.4)
# The parameter a, b, and c are constant whereas the parameter d can change over time
example2 <- data.frame(
time = c(0, 5, 10, 15),
a = c(0.4, NA, NA, NA),
b = c(1.1, NA, NA, NA),
c = c(0.1, NA, NA, NA),
d = c(0.4, 0.5, 0.3, 0.4))
# The parameter a, b are constant
# whereas parameter c and d can change over time.
# However, d is not known for all points of c
example3 <- data.frame(
time = c(0, 5, 10, 15, 20, 25),
a = c(1.1, NA, NA, NA, NA, NA),
b = c(0.1, NA, NA, NA, NA, NA),
c = c(0.2, 0.2, 0, 0, 0, 0),
d = c(0, 0.1, 0.2, 1, NA, NA))
The states data.frame:
___________________________________________________________
The states are defined as a data.frame that contains the 'time' as the first column.
The subsequent columns are the individual states.
# Here some examples
# Only the initial values are defined.
example1 <- data.frame(
time = seq(0, 100, 0.5),
prey = c(10, rep(NA, 200)),
predator = c(10, rep(NA, 200)))
# All values are defined at each timepoint
example2 <- data.frame(
time = seq(0, 100, 0.5),
prey = c(10, runif(200)),
predator = c(10, runif(200)) )
# Only the values for prey are known and are used during optimization
example3 <- data.frame(
time = seq(0, 100, 0.5),
prey = c(10, runif(200)),
predator = c(10, rep(NA, 200)) )
solvertype:
___________________________________________________________
For solving the ode system the SUNDIALS Software is used check the Sundials homepage for more informations.
The solver-type which is used during optimization: “bdf“, “adams“. bdf is an abbreviation for Backward Differentiation Formulas and adams means Adams-Moulton.
All solvers are used in the NORMAL-Step method in a for-loop using the time-points defined in the first column of the 'states' data.frame.
The bdf solver use the SUNLinSol_Dense as linear solver.
own_error_fct:
___________________________________________________________
The error function calculates the error at one of the possible time-points. Moreover, the function expects three numerical scalars as arguments. The first one is the number of data-points at which the error is calculated. The second argument describes the in silico value at one specific time-point. The third argument is the input of the user at the specific time-point which should be matched.
Here is one example shown using the sum of squares as an alternative error function.
error_fct <- function(num_points, insilico, measured) {
ret = (insilico - measured)^2
return(ret/num_points)
}
own_spline_fct:
___________________________________________________________
The spline function is called, directly before the ode-system is evaluated. However, the function is only called for non-constant parameters. See example Nr.2 and Nr.3 parameter d as described above. The results of the spline function is then stored in the vector parameter which is passed to the ode-system. The function expects three arguments:
The time-point at which the function is evaluated
a vector containing the time-points for which parameters are defined
a vector containing the parameters at the respective time-points
The function has to return a scalar value.See the example above for a linear interpolation:
linear_interpolation <- function(t, time_vec, par_vec) {
left = 0
left_time = 0
right = 0
right_time = 0
for(i in 1:length(time_vec)) {
if(t == time_vec[i]) {
return(par_vec[i])
}
if(t < time_vec[i]) {
left = par_vec[i - 1]
right = par_vec[i]
left_time = time_vec[i - 1]
right_time = time_vec[i]
break
}
}
timespan = right_time - left_time
m = (right - left) / timespan
ret = left + m*(t - left_time)
return(ret)
}
Mentionable, is that it hasn't to be a interpolation function. See the example above:
fct <- function(t, time_vec, par_vec) {
ret = 0
for(i in par_vec) {
ret = ret + i
}
return(ret)
}
own_jac_fct:
___________________________________________________________
The jacobian function expects 5 arguments.
the first argument is t which defines the (time-) point of then independent variable at which the ode-system is evaluated.
the second argument is a vector called y which defines the current states at timepoint t
the third argument is a vector called ydot which should be filled with the derivative (left hand side) of the ode-system. It has already the correct length! Please do not return the vector.
the foruth argument is a matrix called J which should be filled with the respective derivatives of ydot. The matrix has already the correct dimensions. This matrix has to be returned.
the last argument is called parameter and is a vector containing the current parameter-set which is tested by the optimization algorithm.
If the parameters can change over time. The already interpolated value is passed to the ode-system.
Value
A list is returned which contains two elements. The first one is the error of the best particle. The other element is a data.frame containing the in silico states returned from the ode-solver using the parameter-set passed by the user..
Note
The error between the defined states and the in silico states is the absolute difference normalised using the true state.
See Also
optimize(),
ast2ast::translate()
Examples
# Solve an ode-system
ode <- function(t, y, ydot, parameter) {
a_db = at(parameter, 1)
b_db = at(parameter, 2)
c_db = at(parameter, 3)
d_db = at(parameter, 4)
predator_db = at(y,1)
prey_db = at(y, 2)
ydot[1] = predator_db*prey_db*c_db - predator_db*d_db
ydot[2] = prey_db*a_db - prey_db*predator_db*b_db
return(ydot)
}
path <- system.file("examples", package = "paropt")
states <- read.table(paste(path,"/states_LV.txt", sep = ""), header = TRUE)
parameter <- data.frame(time = 0, a = 1.1, b = 0.4, c = 0.1, d = 0.4)
res <- paropt::solve(ode,
parameter = parameter,
reltol = 1e-06, abstol = c(1e-08, 1e-08),
states = states, verbose = FALSE)
# solving with own error, spline and jacobian function
jac <- function(t, y, ydot, J, parameter) {
a_db = at(parameter, 1)
b_db = at(parameter, 2)
c_db = at(parameter, 3)
d_db = at(parameter, 4)
predator_db = at(y,1)
prey_db = at(y, 2)
J[1, 1] = prey_db*c_db - d_db
J[2, 1] = - prey_db*b_db
J[1, 2] = predator_db*c_db
J[2, 2] = a_db - predator_db*b_db
return(J)
}
error_fct <- function(c, a, b) {
ret = (a - b)^2
return(ret)
}
spline_fct <- function(t, time_vec, par_vec) {
ret = 0
for(i in par_vec) {
ret = ret + i
}
return(ret)
}
path <- system.file("examples", package = "paropt")
states <- read.table(paste(path,"/states_LV.txt", sep = ""), header = TRUE)
parameter <- data.frame(time = 0, a = 1.1, b = 0.4, c = 0.1, d = 0.4)
res <- paropt::solve(ode,
parameter = parameter,
reltol = 1e-06, abstol = c(1e-08, 1e-08),
states = states, verbose = FALSE)
res <- paropt::solve(ode,
parameter = parameter,
reltol = 1e-06, abstol = c(1e-08, 1e-08),
states = states, verbose = FALSE,
own_error_fct = error_fct,
own_spline_fct = spline_fct,
own_jac_fct = jac)