The package paropt is build in order to optimize parameters of odesystems. Thus, the aim is that the output of the integration matches previously measured states. The user has to supply an odesystem in the form of a Rcppfunction. The information about states and parameters are passed via textfiles. Additional information such as e.g. the relative tolerance are passed in R. This vignette uses a predatorprey system as example to show how the package can be used.
The ode system should be a Rcppfunction with a specific signature. The function is called very often, thus using a R function is not advisable. The name of the function is free to choose. The following parameters have to be passed: a double t, a std::vector
#include <Rcpp.h>
// [[Rcpp::export]]
Rcpp::NumericVector ode_system(double t, std::vector<double> params,
Rcpp::NumericVector states) {
// set dimension of vector which should be returned
Rcpp::NumericVector states_deriv(2); //two states in the example above (prey and predator)
// store parameters in variables
double a = params[0];
double b = params[1];
double c = params[2];
double d = params[3];
// store states in variables
double predator = states[0];
double prey = states[1];
// the actual odesystem
double ddtpredator = states_deriv[0] = predator*c*prey  predator*d;
double ddtprey = states_deriv[1] = prey*a  prey*b*predator;
return states_deriv;
}
The information of the states, have to be supplied as textfiles tab seperated. The delimiter has to be a '.'. It is possible to use scientific notations ('.', '+', '', 'e', 'E' and naturally all digits are allowed). It is compulsory that the first column has the name time. This column contains the independent variable across the solver integrate (it has not to be the time. In this document it is always called time). It is followed by the information of the states at the specific time points. It is pivotal that the order of the states is the same as in the odesystem in order to correctly calculate the error. If a state is not availabe at all timepoints use NA in order to ignore this state at the specified timepoint for calculating the error (see Table1). However, the first line below the header cannot contain NA due to the fact that in this line the start values for the solver are defined.
time  predator  prey 

0  10.00  10.000 
2  2.13  0.022 
4  0.23  3.023 
6  4.30  0.028 
8  0.37  5.000 
10  8.60  0.330 
There exists two different types of parameters: constant and variable parameters. For instance (see Table2) parameters p1, p2 and p3 are variable whereas p4 is constant during the entire time from t = 0 to t = 10.



In the predatorprey example followed here all parameters which should be optimized are constant (see tables below).



The startvalues, lowerbounds and upperbounds for the parameters have to be supplied as textfiles which is tab seperated. The delimiter has to be a '.'. It is possible to use scientific notations ('.', '+', '', 'e', 'E' and naturally all digits are allowed). As for the states it is mandatory that the first column has the name time. This column contains the independent variable across the solver integrate (it has not to be the time) The startvalues have to be between the lower and upper boundary. Naturally, the lowerbounds have to be smaller than the upperbounds.
During the optimization the optimizer creates a bunch of possible solutions which consists of the four different values for a, b, c and d. Each solution is passed to the odesolver which integrate along the time and return the states at the timepoints specified in the textfile containing the stateinformation. The in silico solution is compared to the measured states in order to evaluate the parameterset. The error used is the sum of the absolute differences, between in silico and measured states, divided by the number of timepoints.
During the integration the odesolver can reach a timepoint t where no information for the parameter is measured. For instance in table 2 the timepoint 1.5 is not defined for the parameters p1, p2 and p3. In this case the parameters have to be interpolated. This is conducted using a CatmullRomSpline. The parametervector passed to the odesystem contains already the splined parameters at timepoint t. Thus, using the table above (see Table2) the parametervector could contain for instance at timepoint t = 1.5 (p1 = 1.11, p2 = 1.99, p3 = 0.66, p4 = 0.38) (for the startvalues). Due to that, the parameters passed to the odesystem are in the same order as defined in the textfiles but already interpolated.
When you have the four textfiles defining the states, parameter startvalues, lower and upperbounds and the odesystem you can start running the optimization.
path < system.file("examples", package = "paropt")
library(paropt)
#Rcpp::sourceCpp(paste(path,"/ode.cpp", sep = ""))
#if you want compile odesystem on your system (already precompiled in package)
df < read.table(paste(path,"/states_LV.txt", sep = ""), header = TRUE)
time < df$time
param_start < paste(path, "/start.txt", sep = "")
param_lb < paste(path, "/lb.txt", sep = "")
param_ub < paste(path, "/ub.txt", sep = "")
states < paste(path, "/states_LV.txt", sep = "")
state_output < paste(tempdir(), "/final_stateoutput.txt", sep = "")
par_output < paste(tempdir(), "/optimized_params.txt", sep = "")
set.seed(1)
optimizer(time, paropt:::ode_example, 1e6, c(1e8, 1e8),
param_start, param_lb, param_ub,
states, npop = 40, ngen = 200, error = 1,
state_output, par_output, "bdf")
df_in_silico < read.table(paste(tempdir(), "/final_stateoutput.txt", sep = ""), header = TRUE)
optimizer(time, ode, 1e6, c(1e8, 1e8),
param_start, param_lb, param_ub,
states, npop = 40, ngen = 200, error = 1,
state_output, par_output, "bdf")
The first argument is the timevector containing either the same information as in the timecolumn defined in the textfile containing the states (see Table1) or only a subset (it can only be shorted at the end).
It is mandatory to start with the first entry of the timecolumn, however it is possible to stop at a certain timepoint before the last one. Thus, it is possible to optimize only a part of the problem. Next argument is the compiled odesystem (e.g. see code above which is compiled using the Rcpp::sourceCpp('odesystem.cpp') code in R). The relative tolerance and the absolute tolerances are used by the odesolver. Next the textfiles for the startvalues for the parameters (used for testintegration), the lower and upperbounds and the states are defined. In order to optimize the parameters a particle swarm optimizer (PSO) is used. Therefore, the number of particles (npop = 40 is a good starting point for many problems) and the number of generations (ngen) have to be passed to the function. The number of generations defines the maximum number of generations the PSO should run. However, if the PSO finds a suitable parameterset which posess an error below or equal a threshold defined by the user it stops. This threshold is defined by the errorargument. Next two names of (probably best nonexisiting) textfiles are passed. Using this names the function creates two textfiles. The first one containing the stateoutput of the integration using the optimized parameters and the second one contains the optimized parameters itself. The last argument defines the typ of solver to be used. There exists four different solver: 'bdf', 'ADAMS', 'ERK' and 'ARK. All solvers are part of the SUNDIALS projexct. For details see 'https://computing.llnl.gov/projects/sundials' and Hindmarsh et al. (2005). The bdf solver is the backward differential formular solver of CVODE and is best used for stiff problems. It uses a dense matrix module (SUNDense_Matrix) and the corresponding nonlinear solver (SUNLinsol_Dense). The ADAMS solver is the ADAMSMOULTON solver of CVODE and is most suitable for nonstiff problems. The 'ERK' solver is an explicite RungeKutta solver which is also as the ADAMSMOULTON solver used for nonstiffproblems. The 'ARK' solver is a fully implicte RungeKuttasolver which uses the same matrix and nonlinear solver module as 'bdf'. The integration itself occures for all four solvers in a forloop using the 'CV_NORMAL' stepfunction. If integration for a specific parameterset is not possible the error is set to 1.79769e+308 (which is the maximum of a double).
If you want to test a specific parameterset just call the function ode_solving. The function requires the same parameter as optimizer. Naturally, the arguments for the PSO, the errorthreshold as well as the parameter lower and upperbounds are not needed.
This PSOimplementation is a modified version of 'https://github.com/kthohr/optim' (for a general overview see Sengupta, Basak, and Peters (2018)). The PSO has several key features. First of all a bunch (number of particles defined by user) of possible parametersets is created within the boundaries defined by the user. Each parameterset is called a particle. All particles together are called the swarm. Each possible parameterset is passed to the solver which integrates the system. The result is used to calculate the error. Thus, each particle has a current solution and a current personal best value. Furthermore, each particle possess a neighberhood which consists of 03 other particles (for details see D. Akman, Akman, and Schaefer (2018)). The PSO uses a combination of its history (personal best value) and the information of the best particle of the neighberhood to change its current value. Notably, in this package a randomly adaptive topology is used. This means, that the neighberhood is recalculated each time when the global best solution (best solution of the entire swarm) has not improved within one generation.
Akman, Devin, Olcay Akman, and Elsa Schaefer. 2018. “Parameter Estimation in Ordinary Differential Equations Modeling via Particle Swarm Optimization.” Journal of Applied Mathematics 2018. doi:10.1155/2018/9160793.
Hindmarsh, Alan C., Peter N. Brown, Keith E. Grant, Steven L. Lee, Radu Serban, Dan E. Shumaker, and Carol S. Woodward. 2005. “SUNDIALS: Suite of nonlinear and differential/algebraic equation solvers.” ACM Transactions on Mathematical Software 31 (3): 363–96. doi:10.1145/1089014.1089020.
Sengupta, Saptarshi, Sanchita Basak, and Richard Peters. 2018. “Particle Swarm Optimization: A Survey of Historical and Recent Developments with Hybridization Perspectives.” Machine Learning and Knowledge Extraction 1 (1): 157–91. doi:10.3390/make1010010.