#  --------------------------------------------------------------

monfn <- function(argvals, Wfdobj, basislist=vector("list",JMAX)) {
#  evaluates a monotone function of the form
#            h(x) = [D^{-1} exp Wfdobj](x)
#  where  D^{-1} means taking the indefinite integral.
#  The interval over which the integration takes places is defined in
#  the basis object in Wfdobj.
#  Arguments:
#  ARGVALS   ... argument values at which function and derivatives are evaluated
#  WFDOBJ    ... a functional data object
#  BASISLIST ... a list containing values of basis functions
#  Returns:
#  HVAL   ... matrix or array containing values of h.

#  Last modified 24 April 2009

  JMAX <- 15
  JMIN <- 11
  EPS  <- 1E-5

  coef  <- Wfdobj$coefs
  coefd <- dim(coef)
  ndim  <- length(coefd)
  if (ndim > 1 && coefd[2] != 1) stop("WFDOBJ is not a single function")

  basisobj <- Wfdobj$basis
  rangeval <- basisobj$rangeval

  #  set up first iteration

  width <- rangeval[2] - rangeval[1]
  JMAXP <- JMAX + 1
  h <- rep(1,JMAXP)
  h[2] <- 0.25
  #  matrix SMAT contains the history of discrete approximations to the
  #    integral
  smat <- matrix(0,JMAXP)
  #  array TVAL contains the argument values used in the approximation
  #  array FVAL contains the integral values at these argument values,
  #     rows corresponding to argument values
  #  the first iteration uses just the endpoints
  tval <- rangeval
  j   <- 1
  if (is.null(basislist[[j]])) {
      bmat <- getbasismatrix(tval, basisobj)
      basislist[[j]] <- bmat
  } else {
      bmat <- basislist[[j]]
  }
  fx   <- exp(bmat %*% coef)
  fval <- fx
  smat[1,]  <- width*apply(fx,2,sum)/2
  tnm <- 0.5
  for (j in 2:JMAX) {
    tnm  <- tnm*2
    del  <- width/tnm
    flag <- ifelse(rangeval[1]+del/2 >= rangeval[2]-del/2, -1, 1)
    tj   <- seq(rangeval[1]+del/2, rangeval[2]-del/2, by=flag*abs(del))
    tval <- c(tval, tj)
    if (is.null(basislist[[j]])) {
        bmat <- getbasismatrix(tj, basisobj)
        basislist[[j]] <- bmat
    } else {
        bmat <- basislist[[j]]
    }
    fx   <- exp(bmat %*% coef)
    fval <- c(fval,fx)
    smat[j] <- (smat[j-1] + width*apply(fx,2,sum)/tnm)/2
    if (j >= JMIN) {
      ind <- (j-4):j
      result <- polintmat(h[ind],smat[ind],0)
      ss  <- result[[1]]
      dss <- result[[2]]
      if (all(abs(dss) < EPS*max(abs(ss)))) {
        # successful convergence
        # sort argument values and corresponding function values
        ordind <- order(tval)
        tval   <- tval[ordind]
        fval   <- fval[ordind]
        nx     <- length(tval)
        del    <- tval[2] - tval[1]
        fval   <- del*(cumsum(fval) - 0.5*(fval[1] + fval))
        hval   <- approx(tval, fval, argvals)$y
        return(hval)
      }
    }
    smat[j+1] <- smat[j]
    h[j+1]    <- 0.25*h[j]
  }
  stop(paste("No convergence after",JMAX," steps in MONFN"))
}