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.

> # verify wiggle returns the integeral of the squared second derivative
> 
> library(cpr)
> 
> ################################################################################
> f <- function(x) {
+   #(x + 2) * (x - 1) * (x - 3)
+   x^3 - 2 * x^2 - 5 * x + 6
+ }
> 
> fprime <- function(x) { # first derivatives of f(x)
+   3 * x^2 - 4 * x - 5
+ }
> 
> fdoubleprime <- function(x) { # second derivatives of f(x)
+   6 * x - 4
+ }
> 
> fdoubleprime_squared <- function(x) {
+   fdoubleprime(x)^2
+ }
> 
> bknots = c(-3, 5)
> x     <- (runif(n = 100, min = -3, max = 5))
> bmat  <- bsplines(x, bknots = bknots)
> theta <- matrix(coef(lm(f(x) ~ bmat + 0)), ncol = 1)
> 
> bmatD1 <- bsplineD(x, bknots = bknots, derivative = 1L)
> bmatD2 <- bsplineD(x, bknots = bknots, derivative = 2L)
> 
> stopifnot(isTRUE(all.equal(
+   integrate(fdoubleprime_squared, lower = -3, upper = 5)$value
+   ,
+   wiggle(cp(bmat, theta), lower = -3, upper = 5)$value
+ )))
> 
> ################################################################################
> # verify the result is the same if the bounds are not specified
> stopifnot(isTRUE(all.equal(
+   integrate(fdoubleprime_squared, lower = -3, upper = 5)$value
+   ,
+   wiggle(cp(bmat, theta))$value
+ )))
> 
> ################################################################################
> # verify the wiggle result is the same if the bounds are set too wide
> stopifnot(isTRUE(all.equal(
+   integrate(fdoubleprime_squared, lower = -3, upper = 5)$value
+   ,
+   wiggle(cp(bmat, theta, lower = -10, upper = 20))$value
+ )))
> 
> ################################################################################
> #                              Count Sign Changes                              #
> 
> # sign changes in fprime
> # the data need to be sorted for fprime call
> stopifnot(identical(
+   sum(abs(diff(sign(fprime(sort(x)))))) / 2L
+   ,
+   sign_changes( object = cp(bmat, theta), lower = -3, upper = 5, derivative = 1)
+   )
+ )
> 
> # sign_changes in fdoubleprime
> stopifnot(identical(
+   sum(abs(diff(sign(fdoubleprime(sort(x)))))) / 2L
+   ,
+   sign_changes( object = cp(bmat, theta), lower = -3, upper = 5, derivative = 2)
+   )
+ )
> 
> ################################################################################
> # verify the sign_changes result is the same if the bounds are set too wide
> 
> warn <- tryCatch(sign_changes(cp(bmat, theta), lower = -10, upper = 20, derivative = 1)
+                  , warning = function(w) w)
> stopifnot(inherits(warn, "warning"))
> stopifnot(warn$message == "setting lower = min(object$bknots)")
> 
> 
> stopifnot(isTRUE(all.equal(
+   sign_changes( object = cp(bmat, theta), lower = -3, upper = 5, derivative = 1)
+   ,
+   suppressWarnings(sign_changes(cp(bmat, theta), lower = -10, upper = 20, derivative = 1))
+ )))
> 
> 
> ################################################################################
> #                                 End of File                                  #
> ################################################################################
> 
> proc.time()
   user  system elapsed 
   0.28    0.09    0.37