# R/ode.r In KGode: Kernel Based Gradient Matching for Parameter Inference in Ordinary Differential Equations

#' The 'ode' class object
#'
#' This class provide all information about odes and methods for numerically solving odes.
#' @docType class
#' @importFrom R6 R6Class
#' @import pracma
#' @import mvtnorm
#' @import pspline
#' @export
#' @keywords data
#' @return an  \code{\link{R6Class}} object which can be used for gradient matching.
#' @field ode_par vector(of length n_p) containing ode parameters. n_p is the number of ode parameters.
#' @field ode_fun function containing the ode function.
#' @field t vector(of length n_o) containing time points of observations. n_o is the length of time points.
#' @section Methods:
#' \describe{
#'   \item{\code{solve_ode(par_ode,xinit,tinterv)}}{This method is used to solve ode numerically.}
#'   \item{\code{optim_par(par,y_p,z_p)}}{This method is used to estimate ode parameters by standard gradient matching.}
#'   \item{\code{lossNODE(par,y_p,z_p)}}{This method is used to calculate the mismatching between gradient of interpolation and gradient from ode.}
#'  }
#' @export
#' @examples
#' noise = 0.1  ## set the variance of noise
#' SEED = 19537
#' set.seed(SEED)
#' ## Define ode function, we use lotka-volterra model in this example.
#' ## we have two ode states x[1], x[2] and four ode parameters alpha, beta, gamma and delta.
#' LV_fun = function(t,x,par_ode){
#'   alpha=par_ode[1]
#'   beta=par_ode[2]
#'   gamma=par_ode[3]
#'   delta=par_ode[4]
#'   as.matrix( c( alpha*x[1]-beta*x[2]*x[1] , -gamma*x[2]+delta*x[1]*x[2] ) )
#' }
#' ## Define the gradient of ode function against ode parameters
#' ## df/dalpha,  df/dbeta, df/dgamma, df/ddelta where f is the differential equation.
#' alpha = par[1]; beta= par[2]; gamma = par[3]; delta = par[4]
#' dres= c(0)
#' dres[1] = sum( -2*( z_p[1,]-grad_ode[1,])*y_p[1,]*alpha )
#' dres[2] = sum( 2*( z_p[1,]-grad_ode[1,])*y_p[2,]*y_p[1,]*beta)
#' dres[3] = sum( 2*( z_p[2,]-grad_ode[2,])*gamma*y_p[2,] )
#' dres[4] = sum( -2*( z_p[2,]-grad_ode[2,])*y_p[2,]*y_p[1,]*delta)
#' dres
#' }
#'
#' ## create a ode class object
#' kkk0 = ode$new(2,fun=LV_fun,grfun=LV_grlNODE) #' ## set the initial values for each state at time zero. #' xinit = as.matrix(c(0.5,1)) #' ## set the time interval for the ode numerical solver. #' tinterv = c(0,6) #' ## solve the ode numerically using predefined ode parameters. alpha=1, beta=1, gamma=4, delta=1. #' kkk0$solve_ode(c(1,1,4,1),xinit,tinterv)
#'
#' ## Create another ode class object by using the simulation data from the ode numerical solver.
#' ## If users have experiment data, they can replace the simulation data with the experiment data.
#' ## set initial values for ode parameters.
#' init_par = rep(c(0.1),4)
#' init_yode = kkk0$y_ode #' init_t = kkk0$t
#' kkk = ode$new(1,fun=LV_fun,grfun=LV_grlNODE,t=init_t,ode_par= init_par, y_ode=init_yode ) #' #' @author Mu Niu, \email{ [email protected]} ode<-R6Class("ode", public = list( sample = NULL, t = NULL, y_ode = NULL, ode_par=NULL, ode_fun=NULL, gr_lNODE=NULL, initialize = function(sample = NULL,fun = NULL,grfun=NULL,t=NULL,ode_par = NULL,y_ode=NULL) { self$sample <- sample
self$ode_par <- ode_par self$t<- t
self$y_ode<- y_ode self$ode_fun = fun
self$gr_lNODE = grfun self$greet()
},

greet = function() {
cat(paste0("ode is sample ", self$sample, ".\n")) }, ## ordinary differential equation in matrix form ## ode solver solve_ode = function(par_ode,xinit,tinterv){ x0 = xinit#as.matrix(c(0.5,1)) t0 <- tinterv[1]#0; tf <- tinterv[2]#6 solt <- ode23s(self$ode_fun, t0, tf, x0, par_ode=par_ode,hmax=(tf-t0) )
pick = seq(1,length(solt$t),self$sample)
self$t = solt$t[pick]
rownames(solt$y) <- NULL self$y_ode= t(solt$y[pick,]) }, rmsfun = function(par_ode,state,M1,true_par){ out <- ode(y =state,times =self$t, func = M1,parms=par_ode,method="ode23")
funlos = sqrt( sum( ( t(out[,-1]) - self$y_ode)^2 ) / length(c(self$y_ode)) )
parlos= sqrt( sum( (par_ode - true_par)^2 ) / length(true_par)  )
relalos= sqrt( sum( ( (par_ode - true_par)/true_par )^2 ) / length(true_par)  )
c(funlos,parlos,relalos)
},
## y_p need to be 50*2  apply to each row
ydif = apply(y_p,2,self$ode_fun,t=self$t,par_ode)
ydif
},
lossNODE= function( par,y_p,z_p) {
par_ode= exp(par)
grad_ode = self$gradient(y_p,par_ode) dim(grad_ode) = dim(z_p) lm = max(dim(y_p)) res=0 start = 2 res = sum( (z_p[,start:lm]-grad_ode[,start:lm])^2 ) #res = sum( (z_p-grad_ode)^2 ) res }, grlNODE= function( par,y_p,z_p) { par_ode= exp(par) ## ode par lm = max(dim(y_p)) grad_ode = self$gradient(y_p,par_ode)
start = 2
dres=self$gr_lNODE( par_ode,grad_ode[,start:lm],y_p[,start:lm],z_p[,start:lm] ) #dres=self$gr_lNODE( par_ode,grad_ode,y_p,z_p )
dres
},

loss32NODE= function( par,y_p,z_p) {
par_ode= exp(par)
grad_ode = self$gradient(y_p,par_ode) lm = max(dim(y_p)) res=0 res = sum( (z_p[,]-grad_ode[,])^2 ) res }, grl32NODE= function( par,y_p,z_p) { par_ode= exp(par) ## ode par lm = max(dim(y_p)) grad_ode = self$gradient(y_p,par_ode)
#dres= c(0)
dres=self$gr_lNODE( par_ode,grad_ode,y_p,z_p ) dres }, optim_par = function(par,y_p,z_p){ grlNODE = if (is.null(self$gr_lNODE)) NULL else self$grlNODE op = optim(log(par),self$lossNODE,grlNODE,y_p,z_p,method="L-BFGS-B",control=list(trace=3, REPORT=1))
self$ode_par = exp(op$par)
exp(op\$par)
}

)

)

## Try the KGode package in your browser

Any scripts or data that you put into this service are public.

KGode documentation built on Nov. 17, 2017, 4:25 a.m.