R/Utilities.R

In bvpSolve: Solvers for Boundary Value Problems of Differential Equations

### ============================================================================
### S3 methods
### ============================================================================

## An approximation function for bvpSolve objects that were solved with \
## bvpcol
approx <- function (x, ...) UseMethod("approx")

approx.default <- function (x, ...) {
if (inherits (x, "bvpSolve" ))
  approx.bvpSolve(x,...)
else  
  stats::approx(x,...)
#nextmethod()  
}

### ============================================================================

approx.bvpSolve <- function(x, xout=NULL, ...){
 
  Attr <- attributes(x)
  if (Attr$name != "bvpcol")
    stop("can only use 'approx.bvpSolve' if problem was solved with 'bvpcol'")
  istate <- Attr$istate[-(1:6)]   # first 6 elements have nothing to do with continuation
  rstate <- Attr$rstate
  il <- length(xout)
  if (il <= 0)
    stop ("'approx' requires at least one value to approximate")   

  z <- rep(1, istate[4])    #### WAS 4
  if (Attr$colmod)
    appone <- function(x)
      .Fortran("mappsln", as.double(x),
            result = as.double(z), as.double(rstate), as.integer(istate))$result
  else if (! Attr$bspline)
    appone <- function(x)
      .Fortran("appsln", as.double(x),
            result = as.double(z), as.double(rstate), as.integer(istate))$result
  else 
    appone <- function(x)
      .Fortran("sysappsln", as.double(x), 
            result = as.double(z), as.double(rstate), as.integer(istate))$result
  
  Out <- NULL
  for (i in 1:il)
    Out <- rbind(Out,c(xout[i],appone(xout[i])))
  colnames(Out) <- colnames(x)
  class (Out) <- c( "bvpSolve", "matrix" ) 
  attr(Out, "name") <-  "approx"
  dimnames(Out) <- dimnames(x)
  Out
}

### ============================================================================

print.bvpSolve <- function(x, ...)
   print(as.data.frame(x), ... )

### ============================================================================
### Merge two observed data files; assumed that first column = 'x' and ignored
### ============================================================================

# from 3-columned format (what, where, value) to wide format...
convert2wide <- function(Data) {
    cnames   <- as.character(unique(Data[,1]))

    MAT      <- Data[Data[,1] == cnames[1], 2:3]
    colnames.MAT <- c("x", cnames[1])

    for ( ivar in cnames[-1]) {
      sel <- Data[Data[,1] == ivar, 2:3]
      nt  <- cbind(sel[,1],matrix(nrow = nrow(sel), ncol = ncol(MAT)-1, data = NA),sel[,2])
      MAT <- cbind(MAT, NA)
      colnames(nt) <- colnames(MAT)
      MAT <- rbind(MAT, nt)
      colnames.MAT <- c(colnames.MAT, ivar)
    }
  colnames(MAT) <- colnames.MAT
  return(MAT)
}


mergeObs <- function(obs, Newobs) {

  if (! inherits(Newobs, c("data.frame","matrix")))
    stop ("the elements in 'obs' should be either a 'data.frame' or a 'matrix'")

  if (is.character(Newobs[,1]) | is.factor(Newobs[,1]))
    Newobs <- convert2wide(Newobs)

  obsname <- colnames(obs)

## check if some observed variables in NewObs are already in obs
  newname <- colnames(Newobs)[-1]    # 1st column = x-var and ignored
  ii <- which (newname %in% obsname)
  if (length(ii) > 0)
    obsname <- c(obsname, newname[-ii] )
  else
    obsname <- c(obsname, newname)

## padding with NA of the two datasets
  O1 <- matrix(nrow = nrow(Newobs), ncol = ncol(obs), data = NA)
  O1[ ,1] <- Newobs[,1]
  for (j in ii) {   # obseerved data in common are put in correct position
    jj <- which (obsname == newname[j])
    O1[,jj] <- Newobs[,j+1]
  }
  O1 <- cbind(O1, Newobs[,-c(1,ii+1)] )
  colnames(O1) <- obsname

  nnewcol <- ncol(Newobs)-1 - length (ii)  # number of new columns
  if (nnewcol > 0) {
     O2 <- matrix(nrow = nrow(obs), ncol = nnewcol, data = NA)
     O2 <- cbind(obs, O2)
     colnames(O2) <- obsname
  } else O2 <- obs

  obs <- rbind(O2, O1)
  return(obs)
}

## =============================================================================
## Set the mfrow parameters and whether to "ask" for opening a new device
## =============================================================================

setplotpar <- function(nmdots, dots, nv, ask) {
    if (!any(match(nmdots, c("mfrow", "mfcol"), nomatch = 0))) {
      nc <- min(ceiling(sqrt(nv)),3)
      nr <- min(ceiling(nv/nc),3)
      mfrow <- c(nr, nc)
    }
    else if ("mfcol" %in% nmdots)
        mfrow <- rev(dots$mfcol)
    else mfrow <- dots$mfrow

    if (! is.null(mfrow)) {
      mf <- par(mfrow=mfrow)
    }

   ## interactively wait if there are remaining figures
    if (is.null(ask))
      ask <- prod(par("mfrow")) < nv && dev.interactive()

    return(ask)
}
## =============================================================================
## find a variable
## =============================================================================

selectvar <- function (which, var, NAallowed = FALSE) {

    if (!is.numeric(which)) {
        ln <- length(which)
        # keep ordering...
        Select <- NULL
        for ( i in 1:ln) {
          ss <- which(which[i]==var)
          if (length(ss) ==0 & ! NAallowed) 
            stop("variable ", which[i], " not in var")
          else if (length(ss) == 0)
            Select <- c(Select,NA)
          else
            Select <- c(Select,ss)
        }
    }
    else {
        Select <- which + 1  # "Select now refers to the column number
        if (max(Select) > length(var))
            stop("index in 'which' too large")
        if (min(Select) < 1)
            stop("index in 'which' should be > 0")
    }
  return(Select)
}

# Update range, taking into account neg values for log transformed values
Range <- function(Range, x, log) {
   if (log) 
      x[x<=0] <- min(x[x>0])  # remove zeros
   return( range(Range, x, na.rm = TRUE) )
}

## =============================================================================
## function for checking and expanding arguments in dots (...) with default
## =============================================================================

expanddots <- function (dots, default, n) {
  dots <- if (is.null(dots)) default else dots
  rep(dots, length.out = n)
}

# ks->Th: for xlim and ylim....
expanddotslist <- function (dots, n) {
  if (is.null(dots)) return(dots)
  dd <- if (!is.list(dots )) list(dots) else dots
  rep(dd, length.out = n)
}

## =============================================================================
## functions for expanding arguments in dots  (...)
## =============================================================================

repdots <- function(dots, n) 
  if (is.function(dots)) dots else rep(dots, length.out = n)

setdots <- function(dots, n) lapply(dots, repdots, n)

## =============================================================================
## function for extracting element 'index' from dots  (...)
## =============================================================================

extractdots <- function(dots, index) {
  ret <- lapply(dots, "[", index)
  ret <- lapply(ret, unlist) ## thpe: flatten list (experimental)
  return(ret)
}

### ============================================================================
### Plotting bvpSolve objects
### ============================================================================

plot.bvpSolve <- function (x, ..., which = NULL, ask = NULL, obs = NULL, 
    obspar = list()) {

## check observed data - can be a list
    nobs <- 0

    if (! is.null(obs)) {

      if (!is.data.frame(obs) & is.list(obs)) { # a list with different data sets
       Obs <- obs
       obs <- Obs[[1]]  
       obs.pos <- matrix(nrow = 1, c(1, nrow(obs)))
       if (! inherits(obs, c("data.frame", "matrix")))
         stop ("'obs' should be either a 'data.frame' or a 'matrix'")
       if (length(Obs) > 1)
         for ( i in 2 : length(Obs)) {
           obs <- mergeObs(obs, Obs[[i]])
           obs.pos <- rbind(obs.pos, c(obs.pos[nrow(obs.pos),2] +1, nrow(obs)))
         }
       obsname <- colnames(obs) 
      } else {
       if (is.character(obs[,1]) | is.factor(obs[,1]))   # long format - convert
          obs <- convert2wide(obs)
       obsname <- colnames(obs) 
       if (! inherits(obs, c("data.frame", "matrix")))
         stop ("'obs' should be either a 'data.frame' or a 'matrix'")
       obs.pos <- matrix(nrow = 1, c(1, nrow(obs)))
      }                       
    DD <- duplicated(obsname)
    if (sum(DD) > 0)  
      obs <- mergeObs(obs[,!DD], cbind(obs[,1],obs[,DD]))
    nobs <- nrow(obs.pos)   
    }

## variables to be plotted
    varnames <- colnames(x)
    Which   <- which 
    
    if (is.null(Which) & is.null(obs))  # All variables plotted
      Which <- 1 : (ncol(x)-1)
      
    else if (is.null(Which)) {          # All common variables in x and obs plotted
     Which <- which(varnames %in% obsname)
     Which <- Which[Which != 1]         # remove first element (x-value)
     Which <- varnames[Which]           # names rather than numbers
    } 

## Position of variables to be plotted in "x" 
    t       <- 1     # column with "times" 
    xWhich  <- selectvar(Which, varnames)
    np      <- length(xWhich)

    ldots   <- list(...)
    ndots   <- names(ldots)

## number of figures in a row and interactively wait if remaining figures
    ask <- setplotpar(ndots, ldots, np, ask)
    if (ask) {
      oask <- devAskNewPage(TRUE)
      on.exit(devAskNewPage(oask))
    }

## create two lists: x2:   other bvpSolve objects, 
##                   dots: remaining (plotting) parameters
    x2     <- list()
    dots   <- list()
    nd     <- 0
    nother <- 0            

    if (length(ldots) > 0) 
     for ( i in 1:length(ldots))
      if (inherits(ldots[[i]], "bvpSolve")) {
        x2[[nother <- nother + 1]] <- ldots[[i]]  
        names(x2)[nother] <- ndots[i]
      } else if (! is.null(ldots[[i]])) {
        dots[[nd <- nd+1]] <- ldots[[i]]
        names(dots)[nd] <- ndots[i]
      }
    nmdots <- names(dots)

## check compatibility of all bvpSolve objects    
    if (nother > 0) {
      for ( i in 1:nother) {            
        if (min(colnames(x2[[i]]) == varnames) == 0)
          stop("'x' is not compatible with other bvpSolve objects - colnames not the same")
      }
    } 

    nx <- nother + 1 # total number of bvpSolve objects to be plotted

## Position of variables in "obs" (NA = not observed)
    if (nobs > 0) {
      ObsWhich <- selectvar(varnames[xWhich], obsname, NAallowed = TRUE)
      ObsWhich [ ObsWhich > ncol(obs)] <- NA  
      Obspar <- setdots(obspar, nobs)
    } else 
      ObsWhich <- rep(NA, np)

## plotting parameters : split in plot parameters and point parameters
    plotnames <- c("xlab","ylab","xlim","ylim","main","sub","log","asp",
                   "ann","axes","frame.plot","panel.first","panel.last",
                   "cex.lab","cex.axis","cex.main")    
    
    # plot.default parameters
    ii <- names(dots) %in% plotnames
    dotmain <- dots[ii]
    dotmain <- setdots(dotmain, np)  # expand to np for each plot

    # these are different from the default
    dotmain$xlab <- expanddots(dots$xlab, varnames[t]     , np)
    dotmain$ylab <- expanddots(dots$ylab, ""              , np)
    dotmain$main <- expanddots(dots$main, varnames[xWhich], np)

    # ylim and xlim can be lists and are at least two values
    yylim  <- expanddotslist(dots$ylim, np)
    xxlim  <- expanddotslist(dots$xlim, np)

    # point parameters
    ip <- !names(dots) %in% plotnames
    dotpoints <- dots[ip]
    dotpoints <- setdots(dotpoints, nx)   # expand all dots to nx values

    # these are different from default
    dotpoints$type <- expanddots(dots$type, "l", nx)
    dotpoints$lty  <- expanddots(dots$lty, 1:nx, nx)
    dotpoints$pch  <- expanddots(dots$pch, 1:nx, nx)
    dotpoints$col  <- expanddots(dots$col, 1:nx, nx)
    dotpoints$bg   <- expanddots(dots$bg,  1:nx, nx)
 
## for each output variable (plot)
    iobs <- 0
    for (i in 1 : np) {
      ii <- xWhich[i]     # position of variable in 'x'
      io <- ObsWhich[i]   # position of variable in 'obs'
      
      # plotting parameters for bvpSolve output 1 (opens a plot)
      Dotmain   <- extractdots(dotmain, i)
      Dotpoints <- extractdots(dotpoints, 1)
      
      Xlog <- Ylog <- FALSE
      if (! is.null(Dotmain$log)) { 
        Ylog  <- length(grep("y",Dotmain$log))
        Xlog  <- length(grep("x",Dotmain$log))
      }       
      
      # ranges
      if ( is.null (yylim[[i]])) {
        yrange <- Range(NULL, x[, ii], Ylog)
        if (nother>0) 
         for (j in 1:nother) 
           yrange <- Range(yrange, x2[[j]][,ii], Ylog)
        if (! is.na(io)) yrange <- Range(yrange, obs[,io], Ylog)
          Dotmain$ylim <- yrange
      } else  
        Dotmain$ylim  <- yylim[[i]]
       

      if ( is.null (xxlim[[i]])) {
        xrange <- Range(NULL, x[, t], Xlog)
        if (nother>0) 
         for (j in 1:nother) 
           xrange <- Range(xrange, x2[[j]][,t], Xlog)
        if (! is.na(io)) xrange <- Range(xrange, obs[,1], Xlog)
          Dotmain$xlim <- xrange
      } else  
        Dotmain$xlim  <- xxlim[[i]]
      
      # first bvpSolve object plotted (new plot created)
      do.call("plot", c(alist(x[, t], x[, ii]), Dotmain, Dotpoints))
      
      # if other bvpSolve outputs
      if (nother > 0) 
        for (j in 2:nx)   
          do.call("lines", c(alist(x2[[j-1]][, t], x2[[j-1]][, ii]), 
                  extractdots(dotpoints, j)) )
         
      # if observed variables: select correct pars
      if (! is.na(io))   
           for (j in 1: nobs) 
              if (length (i.obs <- obs.pos[j, 1]:obs.pos[j, 2]) > 0) 
                do.call("points", c(alist(obs[i.obs, 1], obs[i.obs, io]), 
                         extractdots(Obspar, j) ))        
            
    }
}


### ============================================================================

diagnostics.bvpSolve<- function(obj, ...) {
  if (!inherits(obj, "bvpSolve" )) return(NULL)
  Attr <- attributes(obj)
  istate <- Attr$istate
  rstate <- Attr$rstate
  idid <- istate[1]

  if (is.null(istate) || is.null (rstate)) return(NULL)
  cat("\n--------------------\n")
  cat(paste( "solved with ",Attr$name))
  cat("\n--------------------\n")

  if (Attr$name == "bvpshoot") {
    cat("\n---------------------------------------------\n")
    cat("diagnostics of BVP solver ")
    cat("\n---------------------------------------------\n")
    df <- c( "The number of function evaluations              :", 
             "The number of jacobian evaluations +LU decomp   :",	
             "The number of steps                             :",
             "The number of calls to the ivp solver           :")
    printmessage(df, Attr$istate2)

    cat("\n---------------------------------------------\n")
    cat("diagnostics of the last run of the IVP solver ")
    cat("\n---------------------------------------------\n")
      diagnostics.deSolve(obj)
    
  } else if (Attr$name == "bvptwp") {
  
    if (idid ==0)  cat("  Integration was successful.\n") else
			{if (idid < 0 &  (Attr$acdc == FALSE))	cat("  Integration was successful but conditioning parameters NOT stabilized.\n")  else
       cat("  Integration was NOT successful\n")}
    df <- c( "The return code                               :",   #1
             "The number of function evaluations            :",
             "The number of jacobian evaluations            :",
             "The number of boundary evaluations            :",
             "The number of boundary jacobian evaluations   :",
             "The number of steps                           :",
             "The number of mesh resets                     :",
             "The maximal number of mesh points             :",   #2
             "The actual number of mesh points              :",
             "The size of the real work array               :",
             "The size of the integer work array            :"
             )

    printmessage(df, istate[c(1:7,10:13)])

    cat("\n--------------------\n")
    cat(paste( "conditioning pars"))
    cat("\n--------------------\n")

    df <- c( "kappa1  :",
             "gamma1  :",
             "sigma   :",
             "kappa   :",
             "kappa2  :")
    printmessage(df, rstate)
    
    if (Attr$acdc == TRUE)
      cat(paste("    The problem was solved for final eps equal to :",
               Attr$eps[1]," \n"))
  } else if (Attr$name == "bvpcol") {
    if (idid ==1)  cat("  Integration was successful.\n") else
       cat("  Integration was NOT successful\n")
    if (! Attr$colmod)  {
     df <- c( "The return code                                   :",   #1
              "The number of function evaluations                :",
              "The number of jacobian evaluations                :",
              "The number of boundary evaluations                :",
              "The number of boundary jacobian evaluations       :",
              "The number of steps                               :",
              "The actual number of mesh points                  :",
              "The number of collocation points per subinterval  :",
              "The number of equations                           :",
              "The number of components (variables)              :")
       printmessage(df, istate[1:10])
    } else {
     df <- c( "The return code                                   :",   #1
              "The number of function evaluations                :",
              "The number of jacobian evaluations                :",
              "The number of boundary evaluations                :",
              "The number of boundary jacobian evaluations       :",
              "The number of continuation steps                  :",
              "The number of succesfull continuation steps       :",
              "The actual number of mesh points                  :",
              "The number of collocation points per subinterval  :",
              "The number of equations                           :",
              "The number of components (variables)              :")
       printmessage(df, c(idid,Attr$icount[1:6],istate[7:10]))
    cat(paste("    The problem was solved for final eps equal to     :",Attr$eps[1]," \n"))
    }
  }

}

### ============================================================================
## internal helper functions for printing solver return code messages
## these functions are not exported

printmessage <-function(message1, state, message2 = NULL, Nr = 1:length(message1)) {
  if (is.null(message2)) {
    cat("\n", paste(formatC(Nr, "##", width = 2), message1,
              signif(state, digits = getOption("digits")), "\n"), "\n")
  } else {
    cat("\n", paste(formatC(Nr, "##", width = 2), message1,
              signif(state, digits = getOption("digits")), message2, "\n"), "\n")

  }
}