integrate.polySpline <- function(spl=NA, lower=NA, upper=NA,
                                  oob.action=stop) {
  ## Verify that all input arguments are given
  if (any(is.na(spl))) {
    stop("The spline to integrate must be given")
  } else if (any(is.na(lower))) {
    stop("The lower bound(s) to integrate must be given")
  } else if (any(is.na(upper))) {
    stop("The upper bound(s) to integrate must be given")
  }
  ## Verify that the inputs are within the interpolation limits
  if (any(min(lower) < min(spl$knots) |
          max(lower) > max(spl$knots) |
          min(upper) < min(spl$knots) |
          max(upper) > max(spl$knots))) {
    oob.action("The bounds are not within the bounds of the spline")
  }

  ## We won't require that the lower bound is smaller than the upper
  ## bound because there are times that it may be desired to have them
  ## inverted, but we will internally swap the arguments to simplify
  ## the integration.
  negate <- lower > upper
  tmp <- lower
  tmp[negate] <- upper[negate]
  upper[negate] <- lower[negate]
  lower[negate] <- tmp[negate]
  negate <- negate*-1+(!negate)*1
  
  coef <- spl$coefficients
  ## divide the coefficients by their exponents to prevent the need at
  ## every iteration.
  for (i in 2:length(coef[1,])) {
    coef[,i] <- coef[,i]/i
  }
  r <- vector(mode="numeric", length=length(lower))
  for (i in 2:length(spl$knots)) {
    ## Only work on parts of the integral that are requied
    keep <- lower <= spl$knots[i] & upper >= spl$knots[i-1]
    ## Only do the integration if necessary
    if (any(keep)) {
      ## Lower and upper bounds for this piecewise integration
      lb <- pmax(lower[keep], spl$knots[i-1]) - spl$knots[i-1]
      ub <- pmin(upper[keep], spl$knots[i]) - spl$knots[i-1]
      ## This is pulled out of the loop for efficiency (removing the ^1)
      r[keep] <- r[keep] + coef[i-1,1]*(ub-lb)
      for (j in 2:length(coef[i,])) {
        r[keep] <- r[keep] + coef[i-1,j]*(ub^j-lb^j)
      }
    }
  }

  return(negate*r)
}