library("basefun")

### check first and second derivatives
if (require("numDeriv")) {

OR <- 5
cf <- cumsum(runif(OR + 1))
### exclude support b/c of lack of differentiability
(y <- c(1:9, 11:19, 21:29) / 10)

test <- expression({
  yv <- numeric_var("y", support = c(1, 2), bounds = c(.1, 2.9))
  bs1 <- Bernstein_basis(yv, log_first = LF, order = OR, extrapolate = EX)
  h1 <- model.matrix(bs1, data = data.frame(y = y)) %*% cf
  h1p <- model.matrix(bs1, data = data.frame(y = y), deriv = c("y" = 1)) %*% cf
  h1pp <- model.matrix(bs1, data = data.frame(y = y), deriv = c("y" = 2)) %*% cf

  fun <- function(x)
    model.matrix(bs1, data = data.frame(y = x), deriv = c("y" = 0)) %*% cf

  stopifnot(all.equal(sapply(y, function(x) grad(fun, x)), c(h1p),
                      tol = 1e-5))

  fun <- function(x)
    model.matrix(bs1, data = data.frame(y = x), deriv = c("y" = 1)) %*% cf

  stopifnot(all.equal(sapply(y, function(x) grad(fun, x)), c(h1pp),
                      tol = 1e-5))
})

LF <- FALSE
EX <- FALSE

eval(test)

LF <- TRUE
EX <- FALSE

eval(test)

LF <- FALSE
EX <- TRUE

eval(test)

LF <- TRUE
EX <- TRUE

eval(test)

}

N <- 100L
d <- data.frame(x = runif(N), y = gl(5, N / 5))
d$x[1:10] <- NA
d$y[11:20] <- NA

xv <- numeric_var("x", support = c(0.3, .7), bounds = c(0, 1))
OK <- nrow(model.matrix(Bernstein_basis(xv), data = d)) == N & 
      nrow(model.matrix(Legendre_basis(xv), data = d)) == N &
      nrow(model.matrix(log_basis(xv), data = d)) == N &
      nrow(model.matrix(as.basis(~ x + y, data = d), data = d)) == N
stopifnot(OK)