R Under development (unstable) (2024-02-14 r85901 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(cpr)
> require(splines)
Loading required package: splines
> 
> ################################################################################
> e <- new.env()
> with(e, {
+   base_bs <- unclass(splines::bs(0:10, knots = c(2, 2.6, 7.8), Boundary.knots = c(0, 10.2), intercept = TRUE))
+   cpr_bspline <-  unclass(bsplines(0:10, iknots = c(2, 2.6, 7.8), bknots = c(0, 10.2), order = 4))
+ 
+   base_bs |> tail()
+   cpr_bspline |> tail()
+ 
+   stopifnot("Equivalent Basis Matrix between cpr and splines::bs when all(x < max(bknots)" =
+     isTRUE(
+       all.equal(current = cpr_bspline, target = base_bs, check.attributes = FALSE)
+     )
+   )
+ })
> 
> ################################################################################
> # print.cpr_bs
> 
> e <- new.env()
> with(e, {
+   bmat       <- bsplines(runif(n = 100, 1, 100), bknots = c(1, 100), df = 52)
+ 
+   stopifnot(identical(bmat, print(bmat)))
+ 
+   print_bmat <- capture.output(bmat)
+ 
+   stopifnot(any(grepl("First\\s6\\srows:", print_bmat))) # default
+   stopifnot(identical(print_bmat[1], "Basis matrix dims: [100 x 52]"))
+   stopifnot(identical(print_bmat[2], "Order: 4"))
+   stopifnot(identical(print_bmat[3], "Number of internal knots: 48"))
+   stopifnot(identical(print_bmat[4], ""))
+   stopifnot(identical(print_bmat[5], "First 6 rows:"))
+   stopifnot(identical(print_bmat[6], ""))
+ 
+   print_bmat <- capture.output(print(bmat, n = 1000))
+   stopifnot(!any(grepl("First\\s\\d+\\srows:", print_bmat)))
+ 
+   print_bmat <- capture.output(print(bmat, n = 1001))
+   stopifnot(!any(grepl("First\\s\\d+\\srows:", print_bmat)))
+ 
+   print_bmat <- capture.output(print(bmat, n = 10))
+   stopifnot(any(grepl("First\\s10\\srows:", print_bmat)))
+ })
Basis matrix dims: [100 x 52]
Order: 4
Number of internal knots: 48

First 6 rows:

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]      [,11]        [,12]
[1,]    0    0    0    0    0    0    0    0    0     0 0.00000000 0.000000e+00
[2,]    0    0    0    0    0    0    0    0    0     0 0.00000000 0.000000e+00
[3,]    0    0    0    0    0    0    0    0    0     0 0.00000000 0.000000e+00
[4,]    0    0    0    0    0    0    0    0    0     0 0.00000000 1.227005e-05
[5,]    0    0    0    0    0    0    0    0    0     0 0.02241223 4.238390e-01
[6,]    0    0    0    0    0    0    0    0    0     0 0.00000000 0.000000e+00
         [,13]        [,14]     [,15] [,16] [,17] [,18]     [,19]     [,20]
[1,] 0.0000000 0.0000000000 0.0000000     0     0     0 0.0000000 0.0000000
[2,] 0.0000000 0.0000000000 0.0000000     0     0     0 0.0000000 0.0000000
[3,] 0.0000000 0.0000000000 0.0000000     0     0     0 0.2597555 0.5725182
[4,] 0.2906961 0.6546879055 0.0546037     0     0     0 0.0000000 0.0000000
[5,] 0.5531100 0.0006388071 0.0000000     0     0     0 0.0000000 0.0000000
[6,] 0.0000000 0.0000000000 0.0000000     0     0     0 0.0000000 0.0000000
         [,21]       [,22] [,23] [,24]      [,25]     [,26]   [,27]       [,28]
[1,] 0.0000000 0.000000000     0     0 0.06833075 0.6123196 0.31714 0.002209655
[2,] 0.0000000 0.000000000     0     0 0.00000000 0.0000000 0.00000 0.000000000
[3,] 0.1666976 0.001028736     0     0 0.00000000 0.0000000 0.00000 0.000000000
[4,] 0.0000000 0.000000000     0     0 0.00000000 0.0000000 0.00000 0.000000000
[5,] 0.0000000 0.000000000     0     0 0.00000000 0.0000000 0.00000 0.000000000
[6,] 0.0000000 0.000000000     0     0 0.00000000 0.0000000 0.00000 0.000000000
     [,29] [,30] [,31] [,32] [,33] [,34] [,35]        [,36]     [,37]     [,38]
[1,]     0     0     0     0     0     0     0 0.000000e+00 0.0000000 0.0000000
[2,]     0     0     0     0     0     0     0 1.589744e-05 0.2541797 0.6829743
[3,]     0     0     0     0     0     0     0 0.000000e+00 0.0000000 0.0000000
[4,]     0     0     0     0     0     0     0 0.000000e+00 0.0000000 0.0000000
[5,]     0     0     0     0     0     0     0 0.000000e+00 0.0000000 0.0000000
[6,]     0     0     0     0     0     0     0 0.000000e+00 0.0000000 0.0000000
         [,39] [,40] [,41] [,42] [,43] [,44] [,45] [,46] [,47] [,48]
[1,] 0.0000000     0     0     0     0     0     0     0     0     0
[2,] 0.0628301     0     0     0     0     0     0     0     0     0
[3,] 0.0000000     0     0     0     0     0     0     0     0     0
[4,] 0.0000000     0     0     0     0     0     0     0     0     0
[5,] 0.0000000     0     0     0     0     0     0     0     0     0
[6,] 0.0000000     0     0     0     0     0     0     0     0     0
            [,49]      [,50]     [,51]    [,52]
[1,] 0.0000000000 0.00000000 0.0000000 0.000000
[2,] 0.0000000000 0.00000000 0.0000000 0.000000
[3,] 0.0000000000 0.00000000 0.0000000 0.000000
[4,] 0.0000000000 0.00000000 0.0000000 0.000000
[5,] 0.0000000000 0.00000000 0.0000000 0.000000
[6,] 0.0002474507 0.02100171 0.3377149 0.641036
> 
> e <- new.env()
> with(e, {
+   bmat       <- bsplines(runif(n = 2, 1, 100), bknots = c(1, 100), df = 52)
+   print_bmat <- capture.output(bmat)
+ 
+   stopifnot(!any(grepl("First\\s\\d+\\srows:", print_bmat)))
+   stopifnot(identical(print_bmat[1], "Basis matrix dims: [2 x 52]"))
+   stopifnot(identical(print_bmat[2], "Order: 4"))
+   stopifnot(identical(print_bmat[3], "Number of internal knots: 48"))
+   stopifnot(identical(print_bmat[4], ""))
+ })
> 
> ################################################################################
> # Unsorted knots will result in an error
> 
> e <- new.env()
> with(e, {
+   xvec <- runif(n = 25, min = -1, max = 1)
+   iknots <- c(0.34, -0.23)
+   bknots = c(-1, 1)
+ 
+   x <- tryCatch(bsplines(xvec, iknots = iknots, bknots = bknots), error = function(e) {e})
+   stopifnot(inherits(x, "simpleError"))
+   stopifnot(identical(x$message, "Knots are not sorted."))
+ 
+   x <- tryCatch(bsplines(xvec, iknots = sort(iknots), bknots = bknots), error = function(e) {e})
+   stopifnot(inherits(x, "cpr_bs"))
+ 
+   x <- tryCatch(bsplines(xvec, iknots = sort(iknots), bknots = bknots), error = function(e) {e})
+   stopifnot(inherits(x, "cpr_bs"))
+ 
+   x <- tryCatch(bsplines(xvec, iknots = sort(iknots), bknots = c(-0.21, 1)), warning = function(w) {w})
+   stopifnot(inherits(x, "warning"))
+   stopifnot(identical(x$message, "At least one x value < min(bknots)"))
+ 
+   x <- tryCatch(bsplines(xvec, iknots = sort(iknots), bknots = c(-1, 0)), warning = function(w) {w})
+   stopifnot(inherits(x, "warning"))
+   stopifnot(identical(x$message, "At least one x value >= max(bknots)"))
+ 
+ })
> 
> ################################################################################
> # Error if a list is passed to bsplines
> e <- new.env()
> with(e, {
+   x <- tryCatch(bsplines(list(1:10)), error = function(e) { e })
+   stopifnot(inherits(x, "simpleError"))
+   stopifnot("error if list passed to bsplines" = identical(x$message, "x is a list. Use btensor instead of bsplines."))
+ })
> 
> ################################################################################
> # combinations of iknots and df are handled well
> e <- new.env()
> with(e, {
+   xvec <- seq(-1, 1-1e-5, length = 25)
+   bknots <- c(-1, 1)
+ 
+   x <- tryCatch(bsplines(xvec, df = 2, bknots = bknots), warning = function(w) {w})
+   stopifnot(identical(class(x), c("simpleWarning", "warning", "condition")))
+   stopifnot(identical(x$message, "df being set to order"))
+ 
+   x <- suppressWarnings(bsplines(xvec, df = 2, bknots = bknots))
+   y <- bsplines(xvec, df = 4, bknots = bknots)
+   stopifnot(isTRUE(all.equal(x, y, check.attributes = FALSE)))
+ 
+   x <- suppressWarnings(bsplines(xvec, df = 4, order = 5, bknots = bknots))
+   y <- bsplines(xvec, df = 5, order = 5, bknots = bknots)
+   z <- bsplines(xvec, order = 5, bknots = bknots)
+   stopifnot(isTRUE(all.equal(x, y, check.attributes = FALSE)))
+   stopifnot(isTRUE(all.equal(x, z, check.attributes = FALSE)))
+ 
+   x <- tryCatch(bsplines(xvec, iknots = 0.2, df = 6, bknots = bknots), warning = function(w) {w})
+   stopifnot(identical(class(x), c("simpleWarning", "warning", "condition")))
+   stopifnot(identical(x$message, "Both iknots and df defined, using iknots"))
+ 
+   x <- suppressWarnings(bsplines(xvec, iknots = 0.2, df = 6, bknots = bknots))
+   y <- bsplines(xvec, iknots = 0.2, bknots = bknots)
+   stopifnot(isTRUE(all.equal(x, y, check.attributes = FALSE)))
+ 
+   x <- bsplines(xvec, df = 7, bknots = bknots)
+   y <- bsplines(xvec, iknots = trimmed_quantile(xvec, probs = 1:3/4), bknots = bknots)
+   stopifnot(isTRUE(all.equal(x, y, check.attributes = FALSE)))
+ })
> 
> ################################################################################
> # derivatives are as expected
> 
> e <- new.env()
> with(e, {
+   f0 <- function(x) {
+     #(x + 2) * (x - 1) * (x - 3)
+     x^3 - 2 * x^2 - 5 * x + 6
+   }
+   f1 <- function(x) {
+     3 * x^2 - 4 * x - 5
+   }
+   f2 <- function(x) {
+     6 * x - 4
+   }
+ 
+ 
+   #x <- seq(-3, 5, length = 100)
+   x <- sort(runif(n = 100, min = -3, max = 5))
+   bknots = c(-3, 5)
+   bmat <- bsplines(x, bknots = bknots)
+   theta <- coef(lm(f0(x) ~ bsplines(x, bknots = bknots) + 0) )
+ 
+   baseD1 <- splines::spline.des(c(-3, -3, -3, -3, 5, 5, 5, 5), x, derivs = rep(1, length(x)))$design
+   cprD1 <- bsplineD(x, bknots = bknots, derivative = 1L)
+ 
+   head(baseD1); head(cprD1)
+   tail(baseD1); tail(cprD1)
+ 
+   baseD2 <- splines::spline.des(c(-3, -3, -3, -3, 5, 5, 5, 5), x, derivs = rep(2, length(x)))$design
+   cprD2 <- bsplineD(x, bknots = bknots, derivative = 2L)
+ 
+   # verify that I can get cubic correct
+   stopifnot(isTRUE(all.equal(f0(x), as.numeric(bmat %*% theta))))
+ 
+   # verify that I can get the first derivative correct -- this is different from
+   # base R in that, by my construction, the derivative at the boundardy knots
+   # should be NA, as the k-fold knots should result in no differentiability at
+   # these points.
+   stopifnot(isTRUE(all.equal(f1(x), as.numeric(cprD1 %*% theta))))
+ 
+   # verify that I can get the second derivative correct
+   stopifnot(isTRUE(all.equal(f2(x), as.numeric(cprD2 %*% theta))))
+ 
+   if (interactive()) {
+     par(mfrow = c(1, 3))
+     plot(x, f0(x), type = "l")
+     points(x, bmat %*% theta, pch = 2, col = 'blue')
+ 
+     plot(x, f1(x), type = "l")
+     points(x, baseD1 %*% theta, pch = 16, col = 'red', cex = 0.7)
+     points(x, cprD1 %*% theta, pch = 2, col = 'blue')
+ 
+     plot(x, f2(x), type = "l")
+     points(x, baseD2 %*% theta, pch = 16, col = 'red', cex = 0.7)
+     points(x, cprD2 %*% theta, pch = 2, col = 'blue')
+   }
+ 
+   # another set of tests
+   xvec <- seq(-1, 5, length = 100)
+   iknots <- c(-0.8, 0.2, 2.3, 2.5)
+ 
+   xi <- sort(c(rep(range(xvec), 4), iknots))
+ 
+   baseD1 <- spline.des(xi, xvec, derivs = rep(1, length(xvec)))$design
+   baseD2 <- spline.des(xi, xvec, derivs = rep(2, length(xvec)))$design
+ 
+   cprD1 <- bsplineD(xvec, iknots = iknots, derivative = 1L)
+   cprD2 <- bsplineD(xvec, iknots = iknots, derivative = 2L)
+ 
+   stopifnot(isTRUE(
+                    all.equal(
+                              cprD1[-100, ]
+                              ,
+                              baseD1[-100, ]
+                    )
+                    ))
+   stopifnot(isTRUE( all.equal(cprD2[-100, ], baseD2[-100, ])))
+ 
+   x <- tryCatch(bsplineD(xvec, derivative = 1.5), error = function(e) {e})
+   stopifnot(identical(class(x), c("simpleError", "error", "condition")))
+   stopifnot(identical(x$message, "Only first and second derivatives are supported"))
+   rm(x)
+ 
+   x <- tryCatch(bsplineD(xvec, derivative = 3), error = function(e) {e})
+   stopifnot(identical(class(x), c("simpleError", "error", "condition")))
+   stopifnot(identical(x$message, "Only first and second derivatives are supported"))
+ 
+   x <- tryCatch(bsplineD(xvec, order = 2, derivative = 2), error = function(e) {e})
+   stopifnot(identical(class(x), c("simpleError", "error", "condition")))
+   stopifnot(identical(x$message, "(order - 2) <= 0"))
+ })
> 
> ################################################################################
> # verify that order is an integer of at least 2
> e <- new.env()
> with(e, {
+   xvec    <- runif(n = 100, min = -1, max = 5)
+   bknots  <- c(-1, 5)
+   bmat    <- tryCatch(bsplines(xvec, bknots = bknots, order = 1), error = function(e) e)
+   bmatD1  <- tryCatch(bsplineD(xvec, bknots = bknots, order = 1, derivative = 1L), error = function(e) e)
+   bmatD2  <- tryCatch(bsplineD(xvec, bknots = bknots, order = 1, derivative = 2L), error = function(e) e)
+ 
+   stopifnot(inherits(bmat, "simpleError"))
+   stopifnot(inherits(bmatD1, "simpleError"))
+   stopifnot(inherits(bmatD2, "simpleError"))
+ 
+   stopifnot(identical(bmat$message, "order needs to be an integer value >= 2."))
+   stopifnot(identical(bmatD1$message, "order needs to be an integer value >= 2."))
+   stopifnot(identical(bmatD2$message, "order needs to be an integer value >= 2."))
+ 
+   bmat <- tryCatch(bsplines(xvec, bknots = bknots, order = 1.9), error = function(e) e)
+   stopifnot(inherits(bmat, "simpleError"))
+   stopifnot(identical(bmat$message, "order needs to be an integer value >= 2."))
+ 
+   bmat <- tryCatch(bsplines(xvec, bknots = bknots, order = 2.9), error = function(e) e)
+   stopifnot(inherits(bmat, "cpr_bs"))
+   stopifnot(identical(attr(bmat, "order"), 2))
+ 
+ })
> 
> ################################################################################
> # attributes of bsplines
> e <- new.env()
> with(e, {
+   xvec <- runif(n = 100, min = -1, max = 5)
+   iknots <- c(-0.8, 0.2, 2.3, 2.5)
+   bknots <- c(-1.2, 5.2)
+ 
+   bmat   <- tryCatch(bsplines(xvec, iknots = iknots, bknots = bknots), error = function(e) e)
+ 
+   stopifnot(inherits(bmat, "cpr_bs"))
+ 
+   stopifnot(identical(
+       names(attributes(bmat))
+       ,
+       c("dim", "order", "df", "iknots", "bknots", "xi", "xi_star", "class", "call", "environment")
+     )
+   )
+ 
+   stopifnot(identical(attr(bmat, "dim"), c(100L, 8L)))
+ 
+   stopifnot(identical(attr(bmat, "order"), 4))
+ 
+   stopifnot(identical(attr(bmat, "df"), 4 + length(iknots)))
+ 
+   stopifnot(identical(attr(bmat, "iknots"),   iknots))
+ 
+   stopifnot(identical(attr(bmat, "bknots"),   bknots))
+ 
+   stopifnot(identical(attr(bmat, "xi"),   sort(c(rep(bknots, 4), iknots))))
+ 
+   stopifnot(identical(attr(bmat, "derivative"),   NULL))
+ 
+ })
> 
> ################################################################################
> ##                                End of File                                 ##
> ################################################################################
> 
> proc.time()
   user  system elapsed 
   0.42    0.09    0.50