Nothing
### ============================================================================
### mebdfi -- solves differential algebraic and ordinary differential equation
### systems defined in res (DAE) or func (ODE)
### and outputs values for the times in `times'
### on input, y and dy contains the initial values of the state
### variables and rates of changes for times[1]
### parms is a vector of parameters for func. They should not
### change during the integration.
### ============================================================================
mebdfi <- function(y, times, func=NULL, parms, dy=NULL, res=NULL,
nind=c(length(y),0,0), rtol=1e-6, atol=1e-6,
jacfunc=NULL, jacres=NULL,
jactype = "fullint", mass = NULL, verbose=FALSE, tcrit = NULL,
hini=0, ynames=TRUE, maxord =7, bandup=NULL,
banddown=NULL, maxsteps=5000, dllname=NULL, initfunc=dllname,
initpar=parms, rpar=NULL, ipar=NULL,nout=0, outnames=NULL,
forcings=NULL, initforc = NULL, fcontrol=NULL, ...) {
### check input
if (is.null(res) && is.null(func))
stop("either `func' or 'res' must be specified")
if (!is.null(res) && !is.null(func))
stop("either `func' OR 'res' must be specified, not both")
if (is.list(func)) { # a list of compiled codes
if (!is.null(jacfunc) & "jacfunc" %in% names(func))
stop("If 'func' is a list that contains jacfunc, argument 'jacfunc' should be NULL")
if (!is.null(initfunc) & "initfunc" %in% names(func))
stop("If 'func' is a list that contains initfunc, argument 'initfunc' should be NULL")
if (!is.null(dllname) & "dllname" %in% names(func))
stop("If 'func' is a list that contains dllname, argument 'dllname' should be NULL")
if (!is.null(initforc) & "initforc" %in% names(func))
stop("If 'func' is a list that contains initforc, argument 'initforc' should be NULL")
if (!is.null(func$jacfunc)) jacfunc <- func$jacfunc
if (!is.null(func$initfunc)) initfunc <- func$initfunc
if (!is.null(func$dllname)) dllname <- func$dllname
if (!is.null(func$initforc)) initforc <- func$initforc
func <- func$func
}
if (is.list(res)) { #
if (!is.null(jacres) & "jacres" %in% names(res))
stop("If 'res' is a list that contains jacres, argument 'jacres' should be NULL")
if (!is.null(initfunc) & "initfunc" %in% names(res))
stop("If 'res' is a list that contains initfunc, argument 'initfunc' should be NULL")
if (!is.null(dllname) & "dllname" %in% names(res))
stop("If 'res' is a list that contains dllname, argument 'dllname' should be NULL")
if (!is.null(initforc) & "initforc" %in% names(res))
stop("If 'res' is a list that contains initforc, argument 'initforc' should be NULL")
if (!is.null(res$jacres)) jacres <- res$jacres
if (!is.null(res$initfunc)) initfunc <- res$initfunc
if (!is.null(res$dllname)) dllname <- res$dllname
if (!is.null(res$initforc)) initforc <- res$initforc
res <- res$res
}
if (!is.numeric(y))
stop("`y' must be numeric")
n <- length(y)
if (! is.null(times)&&!is.numeric(times))
stop("`times' must be NULL or numeric")
if (!is.null(jacres) && !is.null(jacfunc))
stop("either `jacfunc' OR 'jacres' must be specified, not both")
if (!is.null(func) && !is.function(func) && !is.character(func) && ! inherits(func, "CFunc"))
stop("`func' must be a function, a character vector, of class 'CFunc' or NULL")
if (!is.null(res) && !is.function(res) && !is.character(res) && ! inherits(res, "CFunc"))
stop("`res' must be NULL, a function, or character vector, or of class CFunc")
if ((is.character(res) || is.character(func))&& (is.null(dllname) || !is.character(dllname)))
stop("You need to specify the name of the dll or shared library where res or func can be found (without extension)")
if (!is.null(mass) && !(is.null(res) ))
stop("cannot combine `res' with 'mass' - use 'func' instead, or set 'mass' = NULL")
if (!is.numeric(rtol))
stop("`rtol' must be numeric")
if (!is.numeric(atol))
stop("`atol' must be numeric")
if (!is.null(tcrit) & !is.numeric(tcrit))
stop("`tcrit' must be numeric")
if (!is.null(jacfunc) && !(is.function(jacfunc) ))
stop("`jacfunc' must be a function or NULL")
if (!is.null(jacres) && !(is.function(jacres) || is.character(jacres)))
stop("`jacres' must be a function or character vector")
if (length(atol) > 1 && length(atol) != n)
stop("`atol' must either be a scalar, or as long as `y'")
if (length(rtol) > 1 && length(rtol) != n)
stop("`rtol' must either be a scalar, or as long as `y'")
if (hini < 0)
stop("`hini' must be a non-negative value")
if (hini == 0)
hini = 1e-5 # arbitrary
if (!is.numeric(maxord))
stop("`maxord' must be numeric")
if(maxord < 1 || maxord > 8)
stop("`maxord' must be >1 and <=8")
if (maxsteps<1)
stop("`maxsteps' must be >1 ")
### Jacobian, method flag - note: different from livermore solvers
if (jactype == "fullint" )
imp <- 22 # full, calculated internally
else if (jactype == "fullusr" )
imp <- 21 # full, specified by user function
else if (jactype == "bandusr" )
imp <- 23 # banded, specified by user function
else if (jactype == "bandint" )
imp <- 24 # banded, calculated internally
else stop("'jactype' must be one of 'fullint', 'fullusr', 'bandusr' or 'bandint'")
if (imp %in% c(23,24) && is.null(bandup))
stop("'bandup' must be specified if banded Jacobian")
if (imp %in% c(23,24) && is.null(banddown))
stop("'banddown' must be specified if banded Jacobian")
# if (imp == 23) Jacobian should have empty banddown empty rows
# similar as in vode+daspk only! CHECK IT !
if (imp == 23)
erow<-matrix(ncol = n, nrow = banddown, 0)
else erow<-NULL
if (is.null(banddown))
banddown <-1
if (is.null(bandup ))
bandup <-1
if (is.null(dy))
dy <- rep(0,n)
if (!is.numeric(dy))
stop("`dy' must be numeric")
if (length(nind) != 3)
stop("length of `nind' must be =3")
if (sum(nind) != n)
stop("sum of of `nind' must equal n, the number of equations")
if(length(dy) != length(y))
stop("dy and y should be equally long, and equal n, the number of equations")
### model and Jacobian function
Ynames <- attr(y,"names")
dYnames <- attr(dy,"names")
Res <- NULL
JacRes <- NULL
funtype <- 1
ModelInit <- NULL
flist<-list(fmat=0,tmat=0,imat=0,ModelForc=NULL)
if (!is.null(dllname))
# Karline.... to avoid wrong address to initfunc ... added 24/7/2014
if (sum(duplicated (c(func, initfunc, jacfunc, res, jacres))) > 0)
stop("func, initfunc, jacfunc, res, jacres cannot share the same name")
if (!is.null(dllname)| inherits(func, "CFunc") | inherits(res, "CFunc")) {
if (inherits(initfunc, "CFunc"))
ModelInit <- body(initfunc)[[2]]
else if (is.character(initfunc)) # KS: ADDED THAT to allow absence of initfunc
if (is.loaded(initfunc, PACKAGE = dllname, type = "") ||
is.loaded(initfunc, PACKAGE = dllname, type = "Fortran")) {
ModelInit <- getNativeSymbolInfo(initfunc, PACKAGE = dllname)$address
} else if (initfunc != dllname)
stop(paste("cannot integrate: initfunc not loaded ",initfunc))
# Easier to deal with NA in C-code
if (is.null(initfunc)) initfunc <- NA
if (! is.null(forcings))
flist <- checkforcings(forcings,times,dllname,initforc,verbose,fcontrol)
}
## If res or func is a character vector, make sure it describes
## a function in a loaded dll
if (is.character(res) || is.character(func) || inherits(res, "CFunc") || inherits(func, "CFunc")) {
if (is.character(res)){
resname <- res
if (is.loaded(resname, PACKAGE = dllname)) {
Res <- getNativeSymbolInfo(resname, PACKAGE = dllname)$address
} else stop(paste("cannot integrate: res function not loaded",resname))
} else if (inherits(res, "CFunc")) {
Res <- body(res)[[2]]
} else if (is.character(func)) {
funtype <- 2
if (!is.null(mass)) funtype <- 3
resname <- func
if (is.loaded(resname, PACKAGE = dllname)) {
Res <- getNativeSymbolInfo(resname, PACKAGE = dllname)$address
} else stop(paste("cannot integrate: derivs function not loaded",resname))
} else if (inherits(func, "CFunc")) {
funtype <- 2
Res <- body(func)[[2]]
}
if (!is.null(jacres)) {
if (!is.character(jacres) & !inherits(jacres, "CFunc" ))
stop("If 'res' is dynloaded, so must 'jacres' be")
jacname <- jacres
if (inherits(jacres, "CFunc"))
JacRes <- body(jacres)[[2]]
else if (is.loaded(jacname, PACKAGE = dllname)) {
JacRes <- getNativeSymbolInfo(jacname, PACKAGE = dllname)$address
} else
stop(paste("cannot integrate: Jacobian function jacres not loaded ",jacres))
}
## If we go this route, the number of "global" results is in nout
## and output variable names are in outnames
## THIS DOES NOT WORK IN mebdfi...
Nglobal <- nout
if (nout > 0)
stop(" cannot use mebdfi with a model specified in DLL WITH output variables")
rho <- NULL
if (is.null(outnames))
{ Nmtot <- NULL} else
if (length(outnames) == nout)
{ Nmtot <- outnames} else
if (length(outnames) > nout)
Nmtot <- outnames[1:nout] else
Nmtot <- c(outnames,(length(outnames)+1):nout)
if (is.null(ipar))
ipar<-0
if (is.null(rpar))
rpar<-0
} else {
if (is.null(initfunc))
initpar <- NULL # parameter initialisation not needed if function is not a DLL
## func or res and jac are overruled, either including ynames, or not
## This allows to pass the "..." arguments and the parameters
if (is.null(res) && is.null(mass)) { # res is NOT specified, func is
rho <- environment(func)
Res <- function(time,y,dy) {
if (ynames) attr(y,"names") <- Ynames
FF <-func (time,y,parms,...)
c(dy-unlist(FF[1]), unlist(FF[-1]))
}
Res2 <- function(time,y,dy) {
if (ynames) attr(y,"names") <- Ynames
func (time,y,parms,...)
}
} else if (is.null(res)) { # func with mass
rho <- environment(func)
Res <- function(time,y,dy) {
if (ynames) attr(y,"names") <- Ynames
FF <-func (time,y,parms,...)
c(mass %*% dy-unlist(FF[1]), unlist(FF[-1]))
}
Res2 <- function(time,y,dy) { # just for testing
if (ynames) attr(y,"names") <- Ynames
func (time,y,parms,...)
}
} else { # res is specified
rho <- environment(res)
Res <- function(time,y,dy){
if (ynames) {
attr(y,"names") <- Ynames
attr(dy,"names") <- dYnames
}
unlist(res (time,y,dy,parms,...))
}
Res2 <- function(time,y,dy) {
if(ynames) {
attr(y,"names") <- Ynames
attr(dy,"names") <- dYnames
}
res (time,y,dy,parms,...)
}
}
## the Jacobian
if (! is.null(jacfunc)) { # Jacobian associated with func
tmp <- eval(jacfunc(times[1], y, parms, ...), rho)
if (! is.matrix(tmp))
stop("jacfunc must return a matrix\n")
if (is.null(mass))
JacRes <- function(Rin,y,dy) {
if(ynames) {
attr(y,"names") <- Ynames
attr(dy,"names") <- dYnames
}
JF <- -1* jacfunc(Rin[1],y,parms,...)
if (imp %in% c(23,24)) {
JF[bandup+1,]<-JF[bandup+1,]+1/Rin[2]
JF <- rbind(JF,erow)
} else
JF <-JF + diag(ncol = n, nrow = n, x = 1/Rin[2])
return(JF)
}
else
{
if (imp %in% c(23,24))
stop("cannot combine banded jacobian with mass")
JacRes <- function(Rin,y,dy) {
if(ynames) {
attr(y,"names") <- Ynames
attr(dy,"names") <- dYnames
}
JF <- -1* jacfunc(Rin[1],y,parms,...)
JF <- JF + 1/Rin[2]*mass
return(JF)
}
}
} else if (! is.null(jacres)) { # Jacobian associated with res
tmp <- eval(jacres(times[1], y, dy, parms, 1, ...), rho)
if (! is.matrix(tmp))
stop("jacres must return a matrix\n")
dd <- dim(tmp)
if ((imp ==23 && any(dd != c(bandup+banddown+1,n))) ||
(imp ==21 && any(dd != c(n,n))))
stop("Jacobian dimension not ok")
JacRes <- function(Rin,y,dy) {
if (ynames) {
attr(y,"names") <- Ynames
attr(dy,"names") <- dYnames
}
rbind(erow,jacres(Rin[1],y,dy,parms,Rin[2],...))
}
} else JacRes <- NULL
## Call res once to figure out whether and how many "global"
## results it wants to return and some other safety checks
tmp <- eval(Res2(times[1], y, dy), rho)
if (!is.list(tmp))
stop("Model function must return a list\n")
if (length(tmp[[1]]) != length(y))
stop(paste("The number of derivatives returned by func() (",
length(tmp[[1]]), ") must equal the length of the initial conditions vector (",
length(y), ")", sep = ""))
Nglobal <- if (length(tmp) > 1)
length(unlist(tmp[-1])) else 0
## check for NULL? stop("Problem interpreting model output - check for NULL values")
Nmtot <- attr(unlist(tmp[-1]),"names")
} # is.character(res)
### the mbnd vector
mbnd <- vector("integer",4)
mbnd[] <- 0
if (imp %in% c(23,24)) { # sparse Jacobian
mbnd[1] <- banddown
mbnd[2] <- bandup
mbnd[3] <- bandup + banddown +1
mbnd[4] <- 2*banddown+bandup+1
}
if (imp<23) mbnd[]<-n
### tolerances
if (length(atol)==n) {
if (length(rtol) == n)
itol <- 5
else
itol <- 3
} else {
if (length(rtol)==n)
itol <- 4
else
itol <-2
}
if (imp %in% c(23,21) && is.null(jacfunc) && is.null(jacres))
stop ("mebdfi: *jacfunc* or *jacres* NOT specified; either specify *jacfunc* or *jacres* or change *jactype*")
if (is.null(tcrit)) tcrit = times[length(times)]
### length of rwork and iwork
lrw <- 4 + n*32 + 2*(mbnd[4]*n)
liw <- 14 +n
### calling solver
storage.mode(y) <- storage.mode(dy) <- storage.mode(times) <- "double"
storage.mode(rtol) <- storage.mode(atol) <- "double"
out <- .Call("call_mebdfi", y, dy, times, Res, unlist(initpar),
rtol, atol, as.integer(itol), rho, as.double(tcrit),
as.double(hini), as.integer(maxord), as.integer(maxsteps),
as.integer(nind), JacRes, ModelInit, as.integer(verbose),
as.integer(imp), as.integer(mbnd),
as.integer(liw),as.integer(lrw), as.integer(Nglobal),
as.double (rpar), as.integer(ipar), flist, as.integer(funtype),
as.double(mass), PACKAGE = "deTestSet")
### saving results
istate <- attr(out, "istate")
rstate <- attr(out, "rstate")
out <- saveOutDAE(out, y, dy, n, Nglobal, Nmtot, res, Res2,
iin = c(1,5:14), iout = c(1,5,2,13,3,4,10,19,20,21,18), nr = 1)
## ordinary output variables already estimated
nm <- c("time", if (!is.null(attr(y, "names"))) names(y) else as.character(1:n))
if (Nglobal > 0) nm <- c(nm, if (!is.null(Nmtot)) Nmtot else as.character((n +
1):(n + Nglobal)))
attr(out, "type") <- "mebdfi"
if (verbose) diagnostics(out)
t(out)
}
saveOutDAE <- function (out, y, dy, n, Nglobal, Nmtot, res, Res2,
iin, iout, nr = 4) {
istate <- attr(out,"istate")
istate <- setIstate(istate,iin,iout)
Rstate <- attr(out, "rstate")
rstate <- rep(NA,5)
rstate[1:nr] <- Rstate[1:nr]
dYnames <- attr(dy,"names")
nm <- c("time",
if (!is.null(attr(y,"names"))) names(y) else as.character(1:n))
if (Nglobal > 0) {
dout <- attr(out,"dy")
if (!is.character(res)) { # if a DLL: already done...
out2 <- matrix( nrow=Nglobal, ncol=ncol(out))
for (i in 1:ncol(out2)) {
y <- out[-1,i]
dy <- dout[1:length(y),i]
names(y) <- nm[-1]
names(dy) <- dYnames
out2[, i] <- unlist(Res2(out[1, i], y, dy)[-1]) # KS: Res rather than func
}
out <- rbind(out,out2)
} # end !is.character func
nm <- c(nm,
if (!is.null(Nmtot)) Nmtot else
as.character((n+1) : (n + Nglobal)))
}
attr(out,"istate") <- istate
attr(out, "rstate") <- rstate
class(out) <- c("deSolve","matrix") # a differential equation
dimnames(out) <- list(nm,NULL)
return ( out)
}
## =============================================================================
## Make Istate vector similar for all solvers.
## =============================================================================
setIstate <- function(istate, iin, iout)
{
IstateOut <- rep(NA,21)
IstateOut[iout]<- istate[iin]
IstateOut
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.