R/forcings.R

In deSolve: Solvers for Initial Value Problems of Differential Equations ('ODE', 'DAE', 'DDE')

### ============================================================================
### Check forcing function data set, event inputs and time-lag input
### ============================================================================


checkforcings <- function (forcings, times, dllname,
                           initforc, verbose, fcontrol = list()) {


## Check the names of the initialiser function

 if (is.null(initforc))
   stop(paste("initforc should be loaded if there are forcing functions ",initforc))

 if (inherits (initforc, "CFunc")) {
      ModelForc <- body(initforc)[[2]]
 }  else if (is.loaded(initforc, PACKAGE = dllname, type = "") ||
     is.loaded(initforc, PACKAGE = dllname, type = "Fortran")) {
   ModelForc <- getNativeSymbolInfo(initforc, PACKAGE = dllname)$address
 } else
   stop(paste("initforc should be loaded if there are forcing functions ",initforc))

## Check the type of the forcing function data series

  if (is.data.frame(forcings)) forcings <- list(a=forcings)
  if (! is.list(forcings)) forcings <- list(a=forcings)
  nf <- length(forcings)
  #1 check if each forcing function consists of a 2-columned matrix
  for (i in 1:nf) {
    if (ncol(forcings[[i]]) != 2)
      stop("forcing function data sets should consist of two-column matrix")
  }

## Check the control elements (see optim code)

  con <- list(method="linear", rule = 2, f = 0, ties = "ordered")
  nmsC <- names(con)
  con[(namc <- names(fcontrol))] <- fcontrol
  if (length(noNms <- namc[!namc %in% nmsC]) > 0)
     warning("unknown names in fcontrol: ", paste(noNms, collapse = ", "))

  method <- pmatch(con$method, c("linear", "constant"))
    if (is.na(method))
        stop("invalid interpolation method for forcing functions")
  # 1 if linear, 2 if constant...

## Check the timespan of the forcing function data series

  # time span of forcing function data sets should embrace simulation time...
  # although extrapolation is allowed if con$rule = 2 (the default)
  r_t <- range(times)

  for (i in 1:nf) {
    r_f <- range(forcings[[i]][,1])   # time range of this forcing function

    if (r_f[1] > r_t[1]) {
      if (con$rule == 2) {
        mint <- c(r_t[1],forcings[[i]][1,2] )
        forcings[[i]] <- rbind(mint,forcings[[i]])
        if(verbose)
          warning(paste("extrapolating forcing function data sets to first timepoint",i))
       } else
         stop(paste("extrapolating forcing function data sets to first timepoint",i))
    }

    nr   <- nrow(forcings[[i]])
    if (r_f[2] < r_t[2]) {
      if (con$rule == 2) {
        maxt <- c(r_t[2],forcings[[i]][nr,2] )
        forcings[[i]] <- rbind(forcings[[i]],maxt)
        if(verbose)
          warning(paste("extrapolating forcing function data sets to last timepoint",i))
      } else
        stop(paste("extrapolating forcing function data sets to last timepoint",i))
    }

  }

## Check what needs to be done in case the time series is not "ordered"

   if (!identical(con$ties, "ordered")) { # see approx code

     for (i in 1:nf) {

       x <- forcings[[i]][,1]
       nx <- length(x)
       if (length(ux <- unique(x)) < nx) {  # there are non-unique values
         y <- forcings[[i]][,2]
         ties <- con$tiesn
         if (missing(ties))
           warning("collapsing to unique 'x' values")
          y <- as.vector(tapply(y, x, ties))
          x <- sort(ux)
          forcings[[i]] <- cbind(x, y)

       } else {                             # values are unique, but need sorting
          y <- forcings[[i]][,2]
          o <- order(x)
          x <- x[o]
          y <- y[o]
          forcings[[i]] <-  cbind(x,y)
       }
    } # i
  }

## In case the interpolation is of type "constant" and f not equal to 0
## convert y-series, so that always the left value is taken
  if (method == 2 & con$f != 0) {
     for (i in 1:nf) {
       y <- forcings[[i]][,2]
       YY <- c(y,y[length(y)])[-1]
       forcings[[i]][,2] <- (1-con$f)*y + con$f*YY
     }
  }
## all forcings in one vector; adding index to start/end

  fmat <- tmat <- NULL
  imat <- rep(1,nf+1)

  for (i in 1:nf) {
    # Karline: check for NA in forcing series and remove those
    ii <- apply(forcings[[i]],1,function(x)any(is.na(x)))
    if (sum(ii) > 0) forcings[[i]] <- forcings[[i]][!ii,]
    tmat <- c(tmat, forcings[[i]][,1])
    fmat <- c(fmat, forcings[[i]][,2])
    imat[i+1]<-imat[i]+nrow(forcings[[i]])
  }

  storage.mode(tmat) <- storage.mode(fmat) <- "double"
  storage.mode(imat) <- "integer"

  # DIRTY trick not to inflate the number of arguments:
  # add method (linear/constant) to imat
  return(list(tmat = tmat, fmat = fmat, imat = c(imat, method),
              ModelForc = ModelForc))
}

### ============================================================================
### Check timelags data set - also passes "dllname" now  (not yet used)
### ============================================================================

checklags <- function (lags, dllname) {
  if (!is.null(lags)) {
    lags$islag = 1L
    if (is.null(lags$mxhist))
       lags$mxhist <- 1e4
    if (lags$mxhist <1)
      lags$mxhist <- 1e4
    lags$mxhist<-as.integer(lags$mxhist)
    if (is.null(lags$interpol))   # 1= hermitian, 2 = higher order interpolation
       lags$interpol <- 1
    lags$interpol<-as.integer(lags$interpol)
    lags$isfun <- 0L
  } else
    lags$islag <- 0L
  return(lags)
}