From 14364805ca81cc17c1c60f043727f2b364931186 Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Fri, 30 Sep 2022 11:54:24 -0600 Subject: [PATCH 01/30] add note about rgl graphics --- R/plot.cpr_cn.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/plot.cpr_cn.R b/R/plot.cpr_cn.R index c3ed8ed..4c03f0e 100644 --- a/R/plot.cpr_cn.R +++ b/R/plot.cpr_cn.R @@ -6,6 +6,8 @@ #' surface, or both, for a \code{cpr_cn} objects. The three-dimensional plots #' are generated by either \code{\link[plot3D]{persp3D}} form the \code{plot3D} #' package or \code{\link[rgl]{persp3d}} from the \code{rgl} package. +#' \code{rgl} graphics may or may not work on your system depending on support +#' for OpenGL. #' #' Building complex and customized graphics might be easier for you if you use #' \code{\link{get_surface}} to generate the needed data for plotting. See From e61c9ebfddd40f446950b493dcb942998b45c72c Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Fri, 30 Sep 2022 12:21:22 -0600 Subject: [PATCH 02/30] move examples out of examples/ and into to R/ files --- Makefile | 3 +- R/build_tensor.R | 17 ++++++++- R/cpr.R | 72 +++++++++++++++++++++++++++++++++++++- R/get_spline.R | 26 +++++++++++++- R/update_bsplines.R | 42 +++++++++++++++++++++- examples/build_tensor.R | 13 ------- examples/cpr.R | 68 ----------------------------------- examples/get_surface.R | 24 ------------- examples/update_bsplines.R | 40 --------------------- man/build_tensor.Rd | 5 ++- man/cpr.Rd | 25 +++++++------ man/get_surface.Rd | 1 + man/plot.cpr_cn.Rd | 2 ++ man/update_bsplines.Rd | 1 + 14 files changed, 176 insertions(+), 163 deletions(-) delete mode 100644 examples/build_tensor.R delete mode 100644 examples/cpr.R delete mode 100644 examples/get_surface.R delete mode 100644 examples/update_bsplines.R diff --git a/Makefile b/Makefile index c91dcd8..22be6ae 100644 --- a/Makefile +++ b/Makefile @@ -8,7 +8,6 @@ CRAN = "https://cran.rstudio.com" SRC = $(wildcard $(PKG_ROOT)/src/*.cpp) RFILES = $(wildcard $(PKG_ROOT)/R/*.R) -EXAMPLES = $(wildcard $(PKG_ROOT)/examples/*.R) TESTS = $(wildcard $(PKG_ROOT)/tests/testthat/*.R) VIGNETTES = $(wildcard $(PKG_ROOT)/vignettes/*) RAWDATAR = $(wildcard $(PKG_ROOT)/data-raw/*.R) @@ -22,7 +21,7 @@ all: $(PKG_NAME)_$(PKG_VERSION).tar.gz -e "options(warn = 2)" @touch $@ -.document.Rout: $(RFILES) $(SRC) $(EXAMPLES) $(RAWDATAR) $(VIGNETTES) $(PKG_ROOT)/DESCRIPTION +.document.Rout: $(RFILES) $(SRC) $(RAWDATAR) $(VIGNETTES) $(PKG_ROOT)/DESCRIPTION if [ -e "$(PKG_ROOT)/data-raw/Makefile" ]; then $(MAKE) -C $(PKG_ROOT)/data-raw/; else echo "Nothing to do"; fi Rscript --vanilla --quiet -e "options(repo = c('$(CRAN)', '$(BIOC)'))" \ -e "options(warn = 2)" \ diff --git a/R/build_tensor.R b/R/build_tensor.R index d106ef1..7139705 100644 --- a/R/build_tensor.R +++ b/R/build_tensor.R @@ -5,10 +5,25 @@ #' @param x a matrix, or list of numeric matrices, build the tensor product #' @param ... additional numeric matrices to build the tensor product #' +#' @seealso +#' \code{vignette("cpr-pkg", package = "cpr")} for details on tensor products. +#' #' @return #' A matrix #' -#' @example examples/build_tensor.R +#' @examples +#' +#' A <- matrix(1:4, nrow = 10, ncol = 20) +#' B <- matrix(1:6, nrow = 10, ncol = 6) +#' +#' # Two ways of building the same tensor product +#' tensor1 <- build_tensor(A, B) +#' tensor2 <- build_tensor(list(A, B)) +#' all.equal(tensor1, tensor2) +#' +#' # a three matrix tensor product +#' tensor3 <- build_tensor(A, B, B) +#' str(tensor3) #' #' @export build_tensor <- function(x, ...) { diff --git a/R/cpr.R b/R/cpr.R index a3d045a..d771eb4 100644 --- a/R/cpr.R +++ b/R/cpr.R @@ -22,7 +22,77 @@ #' @param progress show a progress bar. #' @param ... not currently used #' -#' @example examples/cpr.R +#' @examples +#' ############################################################################# +#' # Example 1: find a model for log10(pdg) = f(day) from the spdg data set +#' \donttest{ +#' # need the lme4 package to fit a mixed effect model +#' require(lme4) +#' +#' # construct the initial control polygon. Forth order spline with fifty +#' # internal knots. Remember degrees of freedom equal the polynomial order +#' # plus number of internal knots. +#' init_cp <- cp(log10(pdg) ~ bsplines(day, df = 54) + (1|id), +#' data = spdg, method = lme4::lmer) +#' cpr_run <- cpr(init_cp) +#' plot(cpr_run, color = TRUE) +#' plot(cpr_run, type = "rmse", to = 10) +#' +#' # preferable model is in index 4 +#' preferable_cp <- cpr_run[[4]] +#' } +#' +#' ############################################################################# +#' # Example 2: logistic regression +#' # simulate a binary response Pr(y = 1 | x) = p(x) +#' p <- function(x) { 0.65 * sin(x * 0.70) + 0.3 * cos(x * 4.2) } +#' +#' set.seed(42) +#' x <- runif(2500, 0.00, 4.5) +#' sim_data <- data.frame(x = x, y = rbinom(2500, 1, p(x))) +#' +#' # Define the initial control polygon +#' init_cp <- cp(formula = y ~ bsplines(x, df = 54, bknots = c(0, 4.5)), +#' family = binomial(), +#' method = glm, +#' data = sim_data) +#' +#' # run CPR, preferable model is in index 7 +#' cpr_run <- cpr(init_cp) +#' plot(cpr_run, color = TRUE, type = "rmse", to = 15) +#' +#' # plot the fitted spline and the true p(x) +#' sim_data$pred_select_p <- +#' plogis(predict(cpr_run[[7]], newdata = sim_data)$pred) +#' +#' ggplot2::ggplot(sim_data) + +#' ggplot2::theme_bw() + +#' ggplot2::aes(x = x) + +#' ggplot2::geom_point(mapping = ggplot2::aes(y = y), alpha = 0.1) + +#' ggplot2::geom_line(mapping = ggplot2::aes(y = pred_select_p, color = "pred_select_p")) + +#' ggplot2::stat_function(fun = p, mapping = ggplot2::aes(color = 'p(x)')) +#' +#' # compare to gam and a binned average +#' sim_data$x2 <- round(sim_data$x, digits = 1) +#' bin_average <- +#' lapply(split(sim_data, sim_data$x2), function(x) { +#' data.frame(x = x$x2[1], y = mean(x$y)) +#' }) +#' bin_average <- do.call(rbind, bin_average) +#' +#' ggplot2::ggplot(sim_data) + +#' ggplot2::theme_bw() + +#' ggplot2::aes(x = x) + +#' ggplot2::stat_function(fun = p, mapping = ggplot2::aes(color = 'p(x)')) + +#' ggplot2::geom_line(mapping = ggplot2::aes(y = pred_select_p, color = "pred_select_p")) + +#' ggplot2::stat_smooth(mapping = ggplot2::aes(y = y, color = "gam"), +#' method = "gam", +#' formula = y ~ s(x, bs = "cs"), +#' se = FALSE, +#' n = 1000) + +#' ggplot2::geom_line(data = bin_average +#' , mapping = ggplot2::aes(y = y, color = "bin_average")) +#' #' #' @export cpr <- function(x, keep = -1, p = 2, progress = interactive(), ...) { diff --git a/R/get_spline.R b/R/get_spline.R index 48c9617..d00f3a6 100644 --- a/R/get_spline.R +++ b/R/get_spline.R @@ -102,7 +102,31 @@ get_spline.cpr_cn <- function(x, margin = 1, at, n = 100) { #' #' @seealso \code{\link{get_spline}} #' -#' @example examples/get_surface.R +#' @examples +#' ## Extract the control net and surface from a cpr_cn object. +#' a_cn <- cn(pdg ~ btensor(list(day, age), df = list(15, 3), order = list(3, 2)), +#' data = spdg) +#' +#' cn_and_surface <- get_surface(a_cn, n = 50) +#' str(cn_and_surface, max.level = 2) +#' +#' par(mfrow = c(1, 2)) +#' with(cn_and_surface$cn, +#' plot3D::persp3D(unique(Var1), +#' unique(Var2), +#' matrix(z, +#' nrow = length(unique(Var1)), +#' ncol = length(unique(Var2))), +#' main = "Control Net") +#' ) +#' with(cn_and_surface$surface, +#' plot3D::persp3D(unique(Var1), +#' unique(Var2), +#' matrix(z, +#' nrow = length(unique(Var1)), +#' ncol = length(unique(Var2))), +#' main = "Surface") +#' ) #' #' @export get_surface <- function(x, margin = 1:2, at, n = 100) { diff --git a/R/update_bsplines.R b/R/update_bsplines.R index 9dcac91..d044121 100644 --- a/R/update_bsplines.R +++ b/R/update_bsplines.R @@ -11,7 +11,47 @@ #' @seealso \code{\link[stats]{update}}, \code{\link{bsplines}}, #' \code{\link{btensor}} #' -#' @example examples/update_bsplines.R +#' @examples +#' ########################### Updating a cpr_bs object ########################### +#' # construct a B-spline basis +#' bmat <- bsplines(seq(1, 10, length = 15), df = 5, order = 3) +#' +#' # look at the structure of the basis +#' str(bmat) +#' +#' # change the order +#' str(update_bsplines(bmat, order = 4)) +#' +#' # change the order and the degrees of freedom +#' str(update_bsplines(bmat, df = 12, order = 4)) +#' +#' ########################### Updating a cpr_bt object ########################### +#' # construct a tensor product +#' tpmat <- btensor(list(x1 = seq(0, 1, length = 10), x2 = seq(0, 1, length = 10)), +#' df = list(4, 5)) +#' +#' tpmat +#' +#' # update the degrees of freedom +#' update_btensor(tpmat, df = list(6, 7)) +#' +#' ####### Updating bsplines or btensor on the right and side of a formula ######## +#' +#' f1 <- y ~ bsplines(x, df = 14) + var1 + var2 +#' f2 <- y ~ btensor(x = list(x1, x2), df = list(50, 31), order = list(3, 5)) + var1 + var2 +#' +#' update_bsplines(f1, df = 13, order = 5) +#' update_btensor(f2, df = list(13, 24), order = list(3, 8)) +#' +#' ########################### Updating a cpr_cp object ########################### +#' data(spdg, package = "cpr") +#' init_cp <- cp(pdg ~ bsplines(day, df = 30) + age + ttm, data = spdg) +#' updt_cp <- update_bsplines(init_cp, df = 5) +#' +#' ########################### Updating a cpr_cn object ########################### +#' init_cn <- cn(pdg ~ btensor(list(day, age), df = list(30, 4)) + ttm, data = spdg) +#' updt_cn <- update_btensor(init_cn, df = list(30, 2), order = list(3, 2)) +#' #' #' @export #' @rdname update_bsplines diff --git a/examples/build_tensor.R b/examples/build_tensor.R deleted file mode 100644 index 964fbf2..0000000 --- a/examples/build_tensor.R +++ /dev/null @@ -1,13 +0,0 @@ -# see vignette("cpr-pkg", package = "cpr") for details on tensor products. - -A <- matrix(1:4, nrow = 10, ncol = 20) -B <- matrix(1:6, nrow = 10, ncol = 6) - -# Two ways of building the same tensor product -tensor1 <- build_tensor(A, B) -tensor2 <- build_tensor(list(A, B)) -all.equal(tensor1, tensor2) - -# a three matrix tensor product -tensor3 <- build_tensor(A, B, B) -str(tensor3) diff --git a/examples/cpr.R b/examples/cpr.R deleted file mode 100644 index d9bcc4d..0000000 --- a/examples/cpr.R +++ /dev/null @@ -1,68 +0,0 @@ -############################################################################### -# Example 1: find a model for log10(pdg) = f(day) from the spdg data set -\dontrun{ -# need the lme4 package to fit a mixed effect model -# library(lme4) - -# construct the initial control polygon. Forth order spline with fifty internal -# knots. Remember degrees of freedom equal the polynomial order plus number of -# internal knots. -init_cp <- cp(log10(pdg) ~ bsplines(day, df = 54) + (1|id), - data = spdg, method = lme4::lmer) -cpr_run <- cpr(init_cp) -plot(cpr_run) -plot(cpr_run, type = "rmse", to = 10) - -# preferable model is in index 4 -preferable_cp <- cpr_run[[4]] - -} - -############################################################################### -# Example 2: logistic regression -# simulate a binary response Pr(y = 1 | x) = p(x) -p <- function(x) { 0.65 * sin(x * 0.70) + 0.3 * cos(x * 4.2) } - -set.seed(42) -x <- runif(2500, 0.00, 4.5) -sim_data <- data.frame(x = x, y = rbinom(2500, 1, p(x))) - -# Define the initial control polygon -init_cp <- cp(formula = y ~ bsplines(x, df = 54, bknots = c(0, 4.5)), - family = binomial(), - method = glm, - data = sim_data) - -# run CPR, preferable model is in index 7 -cpr_run <- cpr(init_cp) -plot(cpr_run, color = TRUE, type = "rmse", to = 15) - -# plot the fitted spline and the true p(x) -sim_data$pred_select_p <- plogis(predict(cpr_run[[7]], newdata = sim_data)$pred) - -ggplot2::ggplot(sim_data) + -ggplot2::theme_bw() + -ggplot2::aes(x = x) + -ggplot2::geom_point(mapping = ggplot2::aes(y = y), alpha = 0.1) + -ggplot2::geom_line(mapping = ggplot2::aes(y = pred_select_p, color = "pred_select_p")) + -ggplot2::stat_function(fun = p, mapping = ggplot2::aes(color = 'p(x)')) - -# compare to gam and a binned average -sim_data$x2 <- round(sim_data$x, digits = 1) -bin_average <- - lapply(split(sim_data, sim_data$x2), function(x) { - data.frame(x = x$x2[1], y = mean(x$y)) - }) -bin_average <- do.call(rbind, bin_average) - -ggplot2::ggplot(sim_data) + -ggplot2::theme_bw() + -ggplot2::aes(x = x) + -ggplot2::geom_line(mapping = ggplot2::aes(y = pred_select_p, color = "pred_select_p")) + -ggplot2::stat_smooth(mapping = ggplot2::aes(y = y, color = "gam"), - method = "gam", - formula = y ~ s(x, bs = "cs"), - se = FALSE, - n = 1000) + -ggplot2::geom_line(data = bin_average, mapping = ggplot2::aes(y = y, color = "bin_average")) - diff --git a/examples/get_surface.R b/examples/get_surface.R deleted file mode 100644 index 4f67c2e..0000000 --- a/examples/get_surface.R +++ /dev/null @@ -1,24 +0,0 @@ -## Extract the control net and surface from a cpr_cn object. -a_cn <- cn(pdg ~ btensor(list(day, age), df = list(15, 3), order = list(3, 2)), - data = spdg) - -cn_and_surface <- get_surface(a_cn, n = 50) -str(cn_and_surface, max.level = 2) - -par(mfrow = c(1, 2)) -with(cn_and_surface$cn, - plot3D::persp3D(unique(Var1), - unique(Var2), - matrix(z, - nrow = length(unique(Var1)), - ncol = length(unique(Var2))), - main = "Control Net") - ) -with(cn_and_surface$surface, - plot3D::persp3D(unique(Var1), - unique(Var2), - matrix(z, - nrow = length(unique(Var1)), - ncol = length(unique(Var2))), - main = "Surface") - ) diff --git a/examples/update_bsplines.R b/examples/update_bsplines.R deleted file mode 100644 index 3bf2665..0000000 --- a/examples/update_bsplines.R +++ /dev/null @@ -1,40 +0,0 @@ -########################### Updating a cpr_bs object ########################### -# construct a B-spline basis -bmat <- bsplines(seq(1, 10, length = 15), df = 5, order = 3) - -# look at the structure of the basis -str(bmat) - -# change the order -str(update_bsplines(bmat, order = 4)) - -# change the order and the degrees of freedom -str(update_bsplines(bmat, df = 12, order = 4)) - -########################### Updating a cpr_bt object ########################### -# construct a tensor product -tpmat <- btensor(list(x1 = seq(0, 1, length = 10), x2 = seq(0, 1, length = 10)), - df = list(4, 5)) - -tpmat - -# update the degrees of freedom -update_btensor(tpmat, df = list(6, 7)) - -####### Updating bsplines or btensor on the right and side of a formula ######## - -f1 <- y ~ bsplines(x, df = 14) + var1 + var2 -f2 <- y ~ btensor(x = list(x1, x2), df = list(50, 31), order = list(3, 5)) + var1 + var2 - -update_bsplines(f1, df = 13, order = 5) -update_btensor(f2, df = list(13, 24), order = list(3, 8)) - -########################### Updating a cpr_cp object ########################### -data(spdg, package = "cpr") -init_cp <- cp(pdg ~ bsplines(day, df = 30) + age + ttm, data = spdg) -updt_cp <- update_bsplines(init_cp, df = 5) - -########################### Updating a cpr_cn object ########################### -init_cn <- cn(pdg ~ btensor(list(day, age), df = list(30, 4)) + ttm, data = spdg) -updt_cn <- update_btensor(init_cn, df = list(30, 2), order = list(3, 2)) - diff --git a/man/build_tensor.Rd b/man/build_tensor.Rd index 6eebda8..2dcd568 100644 --- a/man/build_tensor.Rd +++ b/man/build_tensor.Rd @@ -18,7 +18,6 @@ A matrix Tensor products of Matrices. } \examples{ -# see vignette("cpr-pkg", package = "cpr") for details on tensor products. A <- matrix(1:4, nrow = 10, ncol = 20) B <- matrix(1:6, nrow = 10, ncol = 6) @@ -31,4 +30,8 @@ all.equal(tensor1, tensor2) # a three matrix tensor product tensor3 <- build_tensor(A, B, B) str(tensor3) + +} +\seealso{ +\code{vignette("cpr-pkg", package = "cpr")} for details on tensor products. } diff --git a/man/cpr.Rd b/man/cpr.Rd index bd9454e..c46a5ed 100644 --- a/man/cpr.Rd +++ b/man/cpr.Rd @@ -41,27 +41,26 @@ stored in the first \code{keep + 1} (zero internal knots, one internal knot, limit on the number of stored regression fits is to keep memory usage down. } \examples{ -############################################################################### +############################################################################# # Example 1: find a model for log10(pdg) = f(day) from the spdg data set -\dontrun{ +\donttest{ # need the lme4 package to fit a mixed effect model -# library(lme4) +require(lme4) -# construct the initial control polygon. Forth order spline with fifty internal -# knots. Remember degrees of freedom equal the polynomial order plus number of -# internal knots. +# construct the initial control polygon. Forth order spline with fifty +# internal knots. Remember degrees of freedom equal the polynomial order +# plus number of internal knots. init_cp <- cp(log10(pdg) ~ bsplines(day, df = 54) + (1|id), data = spdg, method = lme4::lmer) cpr_run <- cpr(init_cp) -plot(cpr_run) +plot(cpr_run, color = TRUE) plot(cpr_run, type = "rmse", to = 10) # preferable model is in index 4 preferable_cp <- cpr_run[[4]] - } -############################################################################### +############################################################################# # Example 2: logistic regression # simulate a binary response Pr(y = 1 | x) = p(x) p <- function(x) { 0.65 * sin(x * 0.70) + 0.3 * cos(x * 4.2) } @@ -81,7 +80,8 @@ cpr_run <- cpr(init_cp) plot(cpr_run, color = TRUE, type = "rmse", to = 15) # plot the fitted spline and the true p(x) -sim_data$pred_select_p <- plogis(predict(cpr_run[[7]], newdata = sim_data)$pred) +sim_data$pred_select_p <- + plogis(predict(cpr_run[[7]], newdata = sim_data)$pred) ggplot2::ggplot(sim_data) + ggplot2::theme_bw() + @@ -101,12 +101,15 @@ bin_average <- do.call(rbind, bin_average) ggplot2::ggplot(sim_data) + ggplot2::theme_bw() + ggplot2::aes(x = x) + +ggplot2::stat_function(fun = p, mapping = ggplot2::aes(color = 'p(x)')) + ggplot2::geom_line(mapping = ggplot2::aes(y = pred_select_p, color = "pred_select_p")) + ggplot2::stat_smooth(mapping = ggplot2::aes(y = y, color = "gam"), method = "gam", formula = y ~ s(x, bs = "cs"), se = FALSE, n = 1000) + -ggplot2::geom_line(data = bin_average, mapping = ggplot2::aes(y = y, color = "bin_average")) +ggplot2::geom_line(data = bin_average + , mapping = ggplot2::aes(y = y, color = "bin_average")) + } diff --git a/man/get_surface.Rd b/man/get_surface.Rd index ed789cb..c3318ab 100644 --- a/man/get_surface.Rd +++ b/man/get_surface.Rd @@ -48,6 +48,7 @@ with(cn_and_surface$surface, ncol = length(unique(Var2))), main = "Surface") ) + } \seealso{ \code{\link{get_spline}} diff --git a/man/plot.cpr_cn.Rd b/man/plot.cpr_cn.Rd index 67e45ab..93905c5 100644 --- a/man/plot.cpr_cn.Rd +++ b/man/plot.cpr_cn.Rd @@ -56,6 +56,8 @@ This plotting method generates three-dimensional plots of the control net, surface, or both, for a \code{cpr_cn} objects. The three-dimensional plots are generated by either \code{\link[plot3D]{persp3D}} form the \code{plot3D} package or \code{\link[rgl]{persp3d}} from the \code{rgl} package. +\code{rgl} graphics may or may not work on your system depending on support +for OpenGL. Building complex and customized graphics might be easier for you if you use \code{\link{get_surface}} to generate the needed data for plotting. See diff --git a/man/update_bsplines.Rd b/man/update_bsplines.Rd index 1beb3a3..392a3b9 100644 --- a/man/update_bsplines.Rd +++ b/man/update_bsplines.Rd @@ -62,6 +62,7 @@ updt_cp <- update_bsplines(init_cp, df = 5) init_cn <- cn(pdg ~ btensor(list(day, age), df = list(30, 4)) + ttm, data = spdg) updt_cn <- update_btensor(init_cn, df = list(30, 2), order = list(3, 2)) + } \seealso{ \code{\link[stats]{update}}, \code{\link{bsplines}}, From 128566455a9b6a09f363d792203c671ef18c9799 Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Fri, 30 Sep 2022 12:54:25 -0600 Subject: [PATCH 03/30] remove unneeded self referencing namespace --- R/cp.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/R/cp.R b/R/cp.R index 378d8e1..cf68df6 100644 --- a/R/cp.R +++ b/R/cp.R @@ -195,10 +195,10 @@ plot.cpr_cp <- function(x, ..., show_cp = TRUE, show_spline = FALSE, show_xi = T spline_data <- lapply(list(x, ...), function(xx) { b <- xx$bknots - bmat <- cpr::bsplines(seq(b[1], b[2], length = n), - iknots = xx$iknots, - bknots = b, - order = xx$order) + bmat <- bsplines(seq(b[1], b[2], length = n), + iknots = xx$iknots, + bknots = b, + order = xx$order) data.frame(x = seq(b[1], b[2], length = n), y = as.numeric(bmat %*% xx$cp$theta)) }) From 697687318ecf7e25621ee7fa62497b17dbd9d33f Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Fri, 30 Sep 2022 13:02:39 -0600 Subject: [PATCH 04/30] remove unneeded `+0` from example in cpr.R --- R/cp.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/cp.R b/R/cp.R index cf68df6..f0a3cba 100644 --- a/R/cp.R +++ b/R/cp.R @@ -35,7 +35,7 @@ #' #' # via formula #' dat <- data.frame(x = xvec, y = sin((xvec - 2)/pi) + 1.4 * cos(xvec/pi)) -#' cp3 <- cp(y ~ cpr::bsplines(x) + 0, data = dat) +#' cp3 <- cp(y ~ cpr::bsplines(x), data = dat) #' #' # plot the control polygon, spline and target data. #' plot(cp3, show_spline = TRUE) + From 977bdb9707e6b4fb0e588913ceecacf17b5d6891 Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Mon, 27 Nov 2023 17:39:08 -0700 Subject: [PATCH 05/30] upgrade Roxygen version --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index bc22c66..ff7fb82 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -27,5 +27,5 @@ Suggests: qwraps2, rmarkdown, rgl -RoxygenNote: 7.2.1 +RoxygenNote: 7.2.3 VignetteBuilder: knitr From fd6fb6255eea4d6f9a7a613fb8a8f869c1877b2f Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Tue, 28 Nov 2023 08:26:06 -0700 Subject: [PATCH 06/30] upgrade import versions of ggplot2 and Rcpp --- DESCRIPTION | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index ff7fb82..ac733d3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -14,12 +14,12 @@ Encoding: UTF-8 URL: https://github.com/dewittpe/cpr/ LazyData: true Imports: - ggplot2 (>= 2.2.0), + ggplot2 (>= 3.0.0), lme4, plot3D, - Rcpp (>= 0.12.13), + Rcpp (>= 1.0.11), scales -LinkingTo: Rcpp (>= 0.12.13), +LinkingTo: Rcpp (>= 1.0.11), RcppArmadillo Suggests: covr, From 779d0074875de9c5daa7632f1723f991943ecb6b Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Tue, 28 Nov 2023 09:18:21 -0700 Subject: [PATCH 07/30] remove Makevars and Makevars.win after upgrading Rcpp --- src/Makevars | 3 --- src/Makevars.win | 1 - 2 files changed, 4 deletions(-) delete mode 100644 src/Makevars delete mode 100644 src/Makevars.win diff --git a/src/Makevars b/src/Makevars deleted file mode 100644 index 2a4cd15..0000000 --- a/src/Makevars +++ /dev/null @@ -1,3 +0,0 @@ -PKG_CXXFLAGS = -I../inst/include -DARMA_64BIT_WORD -PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) -CXX_STD = CXX11 diff --git a/src/Makevars.win b/src/Makevars.win deleted file mode 100644 index 22ebc63..0000000 --- a/src/Makevars.win +++ /dev/null @@ -1 +0,0 @@ -PKG_LIBS = $(LAPACK_LIBS) $(BLAS_LIBS) $(FLIBS) From 1451ea27780a0e3bfd58e9ecf911049d2f5a99ce Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Tue, 28 Nov 2023 09:18:45 -0700 Subject: [PATCH 08/30] extend clean recipe --- Makefile | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Makefile b/Makefile index 22be6ae..7f7d500 100644 --- a/Makefile +++ b/Makefile @@ -43,4 +43,6 @@ clean: $(RM) -rf $(PKG_NAME).Rcheck $(RM) -f .document.Rout $(RM) -f .install_dev_deps.Rout + $(RM) -f src/*.o + $(RM) -f src/*.so From 602d9cbf23914bf513d725ba8083e2da8a6a5f47 Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Tue, 28 Nov 2023 09:19:05 -0700 Subject: [PATCH 09/30] wip --- NEWS.md | 8 ++++---- man/cp.Rd | 2 +- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/NEWS.md b/NEWS.md index 5192f16..86687cf 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,7 +10,7 @@ ## Other Changes * Depends on Rcpp >= 0.12.11 to handle registering native routines. -* Moves rgl from `Imports` to `Enhances` (re #36) +* Moves rgl from `Imports` to `Suggests` (re #36) * Refactoring base code to eliminate the use of dplyr, tidyr, tibble, etc. Focus on base R methods to reduce install dependencies and improve long term stability of the package. @@ -48,7 +48,7 @@ Documentation improvements. * `is.` a collection of `is.cpr_cp`, `is.cpr_bs`, ... functions added. * The dataset `spdg` has been added to the package. -## Non-User Visible Changes +## Other Changes * removed a redundant `build_tensor` definition # Version 0.2.0 @@ -80,7 +80,7 @@ Continued development should be focused on bug fixes and minor enhancements. * `keep` is correctly handled in the `cnr` call. * `show_xi` correctly handled in the `plot.cpr_cp` call. -## Non-User visible changes +## Non visible changes * non-exported function `knot_expr` created to help with plotting the knot locations in `cpr:::plot.cpr_bs`. @@ -91,7 +91,7 @@ Continued development should be focused on bug fixes and minor enhancements. When plotting multiple control polygons and splines, this option will make it easier to view the spline functions. -# Non-visible changes: +# Non visible changes: * Extended testing scripts. # Version 0.1.0 diff --git a/man/cp.Rd b/man/cp.Rd index 5925c0b..496e9a4 100644 --- a/man/cp.Rd +++ b/man/cp.Rd @@ -110,7 +110,7 @@ plot(cp1, cp2, show_spline = TRUE, color = TRUE) # via formula dat <- data.frame(x = xvec, y = sin((xvec - 2)/pi) + 1.4 * cos(xvec/pi)) -cp3 <- cp(y ~ cpr::bsplines(x) + 0, data = dat) +cp3 <- cp(y ~ cpr::bsplines(x), data = dat) # plot the control polygon, spline and target data. plot(cp3, show_spline = TRUE) + From 40b7db1ad0a08fb2d76f576776cc60b6b7c530ac Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Tue, 28 Nov 2023 09:19:15 -0700 Subject: [PATCH 10/30] improve testing of bsplines --- tests/test-bsplines.R | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/tests/test-bsplines.R b/tests/test-bsplines.R index 4a0c44a..f83ede1 100644 --- a/tests/test-bsplines.R +++ b/tests/test-bsplines.R @@ -8,12 +8,12 @@ bmat <- bsplines(1:100, df = 52) print_bmat <- capture.output(bmat) stopifnot(any(grepl("First\\s\\d+\\srows:", print_bmat))) -stopifnot(print_bmat[1] == "Basis matrix dims: [100 x 52]") -stopifnot(print_bmat[2] == "Order: 4") -stopifnot(print_bmat[3] == "Number of internal knots: 48") -stopifnot(print_bmat[4] == "") -stopifnot(print_bmat[5] == "First 6 rows:") -stopifnot(print_bmat[6] == "") +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))) @@ -21,10 +21,13 @@ stopifnot(!any(grepl("First\\s\\d+\\srows:", print_bmat))) ################################################################################ # Equivalent Basis Matrix stopifnot( - all.equal( - current = unclass(bsplines(0:10, iknots = c(2, 2.6, 7.8), bknots = c(0, 10), order = 4)) - , target = unclass(splines::bs(0:10, knots = c(2, 2.6, 7.8), Boundary.knots = c(0, 10), intercept = TRUE)) - , check.attributes = FALSE) + isTRUE( + all.equal( + current = unclass(bsplines(0:10, iknots = c(2, 2.6, 7.8), bknots = c(0, 10), order = 4)) + , target = unclass(splines::bs(0:10, knots = c(2, 2.6, 7.8), Boundary.knots = c(0, 10), intercept = TRUE)) + , check.attributes = FALSE + ) + ) ) ################################################################################ From 4cf41a348f9ea5955b6f5b0c37c1723864545604 Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Tue, 28 Nov 2023 09:36:27 -0700 Subject: [PATCH 11/30] setting up for a version 0.2.4 release --- DESCRIPTION | 2 +- NEWS.md | 5 +++-- src/cpr.h | 2 +- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index ac733d3..9a4753d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: cpr Title: Control Polygon Reduction -Version: 0.2.3.9005 +Version: 0.2.4 Authors@R: c(person("Peter", "DeWitt", email = "dewittpe@gmail.com", role = c("aut", "cre")), person("Samantha", "MaWhinney", email = "sam.mawhinney@ucdenver.edu", role = c("ths")), diff --git a/NEWS.md b/NEWS.md index 86687cf..935b052 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,4 @@ -# Version 0.2.3.9005 +# Version 0.2.4 ## New Examples * `cpr` has examples @@ -9,13 +9,14 @@ * A new vignette focused on just B-splines has been started. ## Other Changes -* Depends on Rcpp >= 0.12.11 to handle registering native routines. +* Depends on Rcpp >= 0.12.11 (actually moved to >= 1.0.11) to handle registering native routines. * Moves rgl from `Imports` to `Suggests` (re #36) * Refactoring base code to eliminate the use of dplyr, tidyr, tibble, etc. Focus on base R methods to reduce install dependencies and improve long term stability of the package. * Require R > 3.5.0 * Stop using testthat for testing +* Remove use of the tidyr, dplyr # Version 0.2.3 First public release. diff --git a/src/cpr.h b/src/cpr.h index 960c52d..759426d 100644 --- a/src/cpr.h +++ b/src/cpr.h @@ -10,7 +10,7 @@ struct bspline { arma::vec knots; // full knot sequence unsigned int j; // jth spline unsigned int order; // polynomial order - arma::vec spline; // the jth B-spline + arma::vec spline; // the jth B-spline // constructors bspline(); From 7bb39d9f4dc7f40214cb324638e4fd7f4eca481f Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Tue, 28 Nov 2023 09:59:11 -0700 Subject: [PATCH 12/30] improve testing to be more robust --- tests/test-bsplines.R | 49 ++++++++++++++++++-------------- tests/test-btensor.R | 19 +++++++------ tests/test-influence_of.R | 11 ++++--- tests/test-knot_expr.R | 5 +++- tests/test-matrix_rank.R | 4 +++ tests/test-recover-know-spline.R | 3 ++ tests/test-trimmed_quantile.R | 10 +++++-- 7 files changed, 63 insertions(+), 38 deletions(-) diff --git a/tests/test-bsplines.R b/tests/test-bsplines.R index f83ede1..6251de6 100644 --- a/tests/test-bsplines.R +++ b/tests/test-bsplines.R @@ -36,23 +36,26 @@ xvec <- seq(-1, 1, length = 25) iknots <- c(0.34, -0.23) x <- tryCatch(bsplines(xvec, iknots = iknots), warning = function(w) {w}) + stopifnot("testing if warning is given when sorting knots" = - class(x) == c("simpleWarning", "warning", "condition")) -stopifnot("test warning message" = x$message == "Sorting knots") + identical(class(x), c("simpleWarning", "warning", "condition")) + ) + +stopifnot("test warning message" = identical(x$message, "Sorting knots")) -stopifnot("expected bases generated after sorting knots" = +stopifnot("expected bases generated after sorting knots" = isTRUE( all.equal( current = unclass(suppressWarnings(bsplines(xvec, iknots = c(.34, -.23)))) , target = unclass(splines::bs(xvec, knots = iknots, intercept = TRUE)) , check.attributes = FALSE ) + ) ) ################################################################################ x <- tryCatch(bsplines(list(1:10)), error = function(e) { e }) -stopifnot(class(x) == c("simpleError", "error", "condition")) -stopifnot("error if list passed to bsplines" = - x$message == "x is a list. use cpr::btensor instead of cpr::bsplines.") +stopifnot(identical(class(x), c("simpleError", "error", "condition"))) +stopifnot("error if list passed to bsplines" = identical(x$message, "x is a list. use cpr::btensor instead of cpr::bsplines.")) ################################################################################ @@ -61,31 +64,30 @@ stopifnot("error if list passed to bsplines" = xvec <- seq(-1, 1, length = 25) x <- tryCatch(bsplines(xvec, df = 2), warning = function(w) {w}) -stopifnot(class(x) == c("simpleWarning", "warning", "condition")) -stopifnot(x$message == "df being set to order") - +stopifnot(identical(class(x), c("simpleWarning", "warning", "condition"))) +stopifnot(identical(x$message, "df being set to order")) x <- suppressWarnings(bsplines(xvec, df = 2)) y <- bsplines(xvec, df = 4) -stopifnot(all.equal(x, y, check.attributes = FALSE)) +stopifnot(isTRUE(all.equal(x, y, check.attributes = FALSE))) x <- suppressWarnings(bsplines(xvec, df = 4, order = 5)) y <- bsplines(xvec, df = 5, order = 5) z <- bsplines(xvec, order = 5) -stopifnot(all.equal(x, y, check.attributes = FALSE)) -stopifnot(all.equal(x, z, check.attributes = FALSE)) +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), warning = function(w) {w}) -stopifnot(class(x) == c("simpleWarning", "warning", "condition")) -stopifnot(x$message == "Both iknots and df defined, using iknots") +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)) y <- bsplines(xvec, iknots = 0.2) -stopifnot(all.equal(x, y, check.attributes = FALSE)) +stopifnot(isTRUE(all.equal(x, y, check.attributes = FALSE))) x <- bsplines(xvec, df = 7) y <- bsplines(xvec, iknots = trimmed_quantile(xvec, probs = 1:3/4)) -stopifnot(all.equal(x, y, check.attributes = FALSE)) +stopifnot(isTRUE(all.equal(x, y, check.attributes = FALSE))) ################################################################################ @@ -101,15 +103,18 @@ baseD2 <- spline.des(xi, xvec, derivs = rep(2, length(xvec)))$design cprD1 <- bsplineD(xvec, iknots = iknots) cprD2 <- bsplineD(xvec, iknots = iknots, derivative = 2L) -all.equal(cprD1[-100, ], baseD1[-100, ]) -all.equal(cprD2[-100, ], baseD2[-100, ]) +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(class(x) == c("simpleError", "error", "condition")) -stopifnot(x$message == "Only first and second derivatives are supported") +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(class(x) == c("simpleError", "error", "condition")) -stopifnot(x$message == "Only first and second derivatives are supported") +stopifnot(identical(class(x), c("simpleError", "error", "condition"))) +stopifnot(identical(x$message, "Only first and second derivatives are supported")) +################################################################################ +## End of File ## +################################################################################ diff --git a/tests/test-btensor.R b/tests/test-btensor.R index e61f766..87e7c87 100644 --- a/tests/test-btensor.R +++ b/tests/test-btensor.R @@ -3,16 +3,16 @@ library(cpr) test_warning <- function(x = NULL, msg = NULL) { stopifnot(!is.null(x)) stopifnot(!is.null(msg)) - stopifnot(class(x) == c("simpleWarning", "warning", "condition")) - stopifnot(x$message == msg) + stopifnot(identical(class(x), c("simpleWarning", "warning", "condition"))) + stopifnot(identical(x$message, msg)) invisible(TRUE) } test_error <- function(x, msg = NULL) { stopifnot(!is.null(x)) stopifnot(!is.null(msg)) - stopifnot(class(x) == c("simpleError", "error", "condition")) - stopifnot(x$message == msg) + stopifnot(identical(class(x), c("simpleError", "error", "condition"))) + stopifnot(identical(x$message, msg)) invisible(TRUE) } @@ -32,7 +32,7 @@ mm <- splines::bs(mtcars$hp, knots = c(100, 150), intercept = TRUE) ) -stopifnot(all.equal(mm, unclass(bm), check.attributes = FALSE)) +stopifnot(isTRUE(all.equal(mm, unclass(bm), check.attributes = FALSE))) ################################################################################ @@ -48,7 +48,7 @@ mm <- splines::bs(mtcars$hp, knots = c(100, 150), intercept = TRUE) ) -stopifnot(all.equal(mm, unclass(bm), check.attributes = FALSE)) +stopifnot(isTRUE(all.equal(mm, unclass(bm), check.attributes = FALSE))) ################################################################################ # 3D btensor matrix is constructed as expected @@ -63,7 +63,7 @@ mm <- splines::bs(mtcars$hp, knots = c(100, 150), intercept = TRUE) : splines::bs(mtcars$mpg, knots = c(12.2, 16.3, 21.9), intercept = TRUE)) -stopifnot(all.equal(mm, unclass(bm), check.attributes = FALSE)) +stopifnot(isTRUE(all.equal(mm, unclass(bm), check.attributes = FALSE))) ################################################################################ @@ -73,11 +73,12 @@ bm <- btensor(x = list(mtcars$disp, mtcars$hp, mtcars$mpg), iknots = list(numeric(0), c(100, 150), c(12.2, 16.3, 21.9))) -stopifnot( +stopifnot(isTRUE( all.equal( lapply(attr(bm, "bspline_list"), attr, which = "bknots"), lapply(list(mtcars$disp, mtcars$hp, mtcars$mpg), range) ) + ) ) x <- @@ -95,3 +96,5 @@ x <- tryCatch(btensor(x = list(mtcars$disp, mtcars$hp, mtcars$mpg), bknots = lis test_error(x, "length(bknots) == length(x) is not TRUE") ################################################################################ +## End of File ## +################################################################################ diff --git a/tests/test-influence_of.R b/tests/test-influence_of.R index 6f11b81..885fff8 100644 --- a/tests/test-influence_of.R +++ b/tests/test-influence_of.R @@ -6,8 +6,8 @@ bmat <- bsplines(x = seq(0, 6, length = 500), iknots = c(1.0, 1.5, 2.3, 4.0, 4.5)) theta <- c(1, 0, 3.5, 4.2, 3.7, -0.5, -0.7, 2, 1.5) omit_xi <- influence_of(cp(bmat, theta), c(6, 8)) -stopifnot(all.equal(omit_xi$weight$w[1], 0.5391355923)) #sprintf("%.10f", omit_xi$weight$w[1]) -stopifnot(all.equal(omit_xi$weight$w[2], 0.2775424520)) +stopifnot(isTRUE(all.equal(omit_xi$weight$w[1], 0.5391355923))) #sprintf("%.10f", omit_xi$weight$w[1]) +stopifnot(isTRUE(all.equal(omit_xi$weight$w[2], 0.2775424520))) ################################################################################ # test that influence weights for control polygon is as expected @@ -16,7 +16,7 @@ bmat <- bsplines(x = seq(0, 6, length = 500), theta <- c(1, 0, 3.5, 4.2, 3.7, -0.5, -0.7, 2, 1.5) iw <- influence_weights(cp(bmat, theta)) expected <- structure(list(iknots = c(1, 1.5, 2.3, 4, 4.5), w = c(1.28320371711094, 0.539135592305278, 0.558614618870046, 0.277542452002693, 0.647979474076723)), class = "data.frame", row.names = c(NA, -5L)) -stopifnot(all.equal(iw, expected)) +stopifnot(isTRUE(all.equal(iw, expected))) ################################################################################ # test that influence_weights for control net are as expected @@ -50,6 +50,9 @@ expected <- ) ) -stopifnot(all.equal(iw, expected)) +stopifnot(isTRUE(all.equal(iw, expected))) +################################################################################ +## End of File ## +################################################################################ diff --git a/tests/test-knot_expr.R b/tests/test-knot_expr.R index d3bdaca..3248983 100644 --- a/tests/test-knot_expr.R +++ b/tests/test-knot_expr.R @@ -4,5 +4,8 @@ library(cpr) bmat <- bsplines(mtcars$hp, df = 8) ke <- cpr:::knot_expr.cpr_bs(bmat, digits = 1) -stopifnot(all.equal(ke$xi_expr[[4]], bquote(xi[7]))) +stopifnot(isTRUE(all.equal(ke$xi_expr[[4]], bquote(xi[7])))) +################################################################################ +## End of File ## +################################################################################ diff --git a/tests/test-matrix_rank.R b/tests/test-matrix_rank.R index a0c2044..2e82c23 100644 --- a/tests/test-matrix_rank.R +++ b/tests/test-matrix_rank.R @@ -16,3 +16,7 @@ stopifnot(identical(matrix_rank(bmat), 15)) bmat <- bsplines(seq(0, 1, length = 100), iknots = c(0.001, 0.002)) stopifnot(identical(ncol(bmat), 6L)) stopifnot(identical(matrix_rank(bmat), 5)) + +################################################################################ +## End of File ## +################################################################################ diff --git a/tests/test-recover-know-spline.R b/tests/test-recover-know-spline.R index 3082200..63d084d 100644 --- a/tests/test-recover-know-spline.R +++ b/tests/test-recover-know-spline.R @@ -53,3 +53,6 @@ recover_spline <- function(k = 4L, start_with = 100L, seed, theta_dist_sd = 100, set.seed(42) stopifnot(recover_spline(start_with = 40L, progress = FALSE)$recovered) +################################################################################ +## End of File ## +################################################################################ diff --git a/tests/test-trimmed_quantile.R b/tests/test-trimmed_quantile.R index 1dc8c62..00fda04 100644 --- a/tests/test-trimmed_quantile.R +++ b/tests/test-trimmed_quantile.R @@ -12,11 +12,12 @@ e <- try(trimmed_quantile(1:100, trim = -3.9, prob = 1:23 / 24, name = FALSE), s stopifnot(inherits(e, "try-error")) stopifnot(attr(e, "condition")$message == "(converted from warning) Overruling trim less than 1 with trim = 1L") -stopifnot( +stopifnot(isTRUE( all.equal(suppressWarnings(trimmed_quantile(1:100, trim = 3.9, prob = 1:23 / 24, name = FALSE)) , quantile(4:97, prob = 1:23 / 24) , check.attributes = FALSE ) + ) ) ################################################################################ @@ -27,6 +28,9 @@ ux <- unique(x) xt <- x[!(x %in% range(x))] uxt <- ux[!(ux %in% range(ux))]#-c(which.min(ux), which.max(ux))] -stopifnot(all.equal(trimmed_quantile(x), quantile(uxt))) -stopifnot(all.equal(trimmed_quantile(x, use_unique = FALSE), quantile(xt))) +stopifnot(isTRUE(all.equal(trimmed_quantile(x), quantile(uxt)))) +stopifnot(isTRUE(all.equal(trimmed_quantile(x, use_unique = FALSE), quantile(xt)))) +################################################################################ +## End of File ## +################################################################################ From 1c0427c7288cd1d62e8d6821c99f3d2365cd5b60 Mon Sep 17 00:00:00 2001 From: Peter DeWitt Date: Tue, 28 Nov 2023 10:09:00 -0700 Subject: [PATCH 13/30] use old vignette for minor release --- NEWS.md | 2 +- vignettes/bsplines.Rnw | 417 ------------ vignettes/cpr-pkg.Rmd | 221 +++++++ vignettes/cpr-pkg.Rnw | 1358 ---------------------------------------- 4 files changed, 222 insertions(+), 1776 deletions(-) delete mode 100644 vignettes/bsplines.Rnw create mode 100644 vignettes/cpr-pkg.Rmd delete mode 100644 vignettes/cpr-pkg.Rnw diff --git a/NEWS.md b/NEWS.md index 935b052..b43504a 100644 --- a/NEWS.md +++ b/NEWS.md @@ -92,7 +92,7 @@ Continued development should be focused on bug fixes and minor enhancements. When plotting multiple control polygons and splines, this option will make it easier to view the spline functions. -# Non visible changes: +# Non visible changes * Extended testing scripts. # Version 0.1.0 diff --git a/vignettes/bsplines.Rnw b/vignettes/bsplines.Rnw deleted file mode 100644 index 05f32a7..0000000 --- a/vignettes/bsplines.Rnw +++ /dev/null @@ -1,417 +0,0 @@ -%\VignetteEngine{knitr::knitr} -%\VignetteIndexEntry{B-Splines and Control Polygons} -%\VignetteEncoding{UTF-8} - -% ---------------------------------------------------------------------------- % -\documentclass[10pt]{article} -\usepackage{fullpage} -\usepackage[T1]{fontenc} -\usepackage[utf8]{inputenc} -\usepackage{alltt} - -% ---------------------------------------------------------------------------- % -% Float setup -\usepackage{ctable} -\usepackage{floatrow} -\floatsetup[table]{framestyle=colorbox,framefit=yes,heightadjust=all,framearound=all,capposition=top} -\floatsetup[figure]{framestyle=colorbox,framefit=yes,heightadjust=all,framearound=all,capposition=bottom} -\usepackage{subcaption} - -% ---------------------------------------------------------------------------- % -% JSS style commands: \proglang, \pkg, and \code -\makeatletter -\@ifundefined{proglang}{% - \let\proglang=\textsf - \newcommand{\pkg}[1]{{\fontseries{b}\selectfont #1}} - \newcommand\code{\bgroup\@makeother\_\@makeother\~\@makeother\$\@codex} - \def\@codex#1{{\normalfont\ttfamily\hyphenchar\font=-1 #1}\egroup} -} -\makeatother - -% ---------------------------------------------------------------------------- % -% Math Stuff -\usepackage{amsmath,amssymb,amsfonts,mathrsfs} - -% ---------------------------------------------------------------------------- % -% new commands -\newcommand{\ie}{{\itshape i.e.}} -\newcommand{\eg}{{\itshape e.g.}} -\newcommand{\bs}[1]{\boldsymbol{#1}} -\newcommand{\card}[1]{\left|#1\right|} -\newcommand{\CP}{C\!P} % control polygon, or control point sequence -\newcommand{\CN}{C\!N} % control net -\DeclareMathOperator*{\argmin}{arg\,min} -\DeclareMathOperator*{\argmax}{arg\,max} - -% ---------------------------------------------------------------------------- % -% Bibtex -\usepackage{filecontents} -\usepackage[authoryear]{natbib} - - -% ---------------------------------------------------------------------------- % -% Front Matter -\title{B-Splines and Control Polygons} -\author{Peter DeWitt} - -% ---------------------------------------------------------------------------- % -% Begin the document -\begin{document} -\maketitle -\abstract{% - This vignette covers general theory and mathematics of B-Splines and their - corresponding control polygons. The use of the \pkg{cpr} package to build and - work with B-Splines and control polygons is presented to show examples of the - general concepts. Additional details about the control polygon reduction - approach for finding a parsimonious set of knots for B-spline regression model - are presented in \code{vignette(topic = "cpr-pkg", package = "cpr")}. - -<>= -# knitr::opts_chunk$set() -hook_source = knitr::knit_hooks$get('source') -knitr::knit_hooks$set(source = function(x, options) { - txt = hook_source(x, options) - # extend the default source hook - gsub('\\\\textbackslash\\{\\}ref\\\\\\{(.+)\\\\\\}', '\\\\ref{\\1}', txt) -}) -@ - -<<>>= -set.seed(42) -library(cpr) -@ -} - -\section{B-Splines} -\subsection{General Definition} -Many references are available for B-splines. Two very good references are -\citet{deboor2001} and \citet{prautzsch2002}. The goal of this document is to -provide a detailed exploration of B-splines in the context useful for -statistical models. - -B-splines get there name from `basis'-splines. That this, a set of basis -functions define a B-spline basis and space. A B-spline basis, $\bs{B}_{k, -\bs{\xi}} \left(\bs{x}\right),$ is defined by a polynomial -order\footnote{Recall: the order of a polynomial is the degree plus 1, \eg, -$x^n$ is a $n^{th}$ degree, or $(n+1)^{st}$ order, polynomial.} -$k$ and knot sequence $\bs{\xi}$ over the range of a numeric vector $\bs{x}.$ -A traditional knot sequence construction places $k$ fold boundary knots at -$\min\left(\bs{x}\right)$ and $\max\left(\bs{x}\right),$ with an additional $l -\geq 0$ internal knots located within $\text{range}\left(\bs{x}\right).$ -The basis matrix, $\bs{B}_{k, \bs{\xi}} \left( x\right) \in -\mathbb{R}^{\card{\bs{x}} \times \left(k + l\right)}$ is of the form -\begin{equation} - \label{eq:basis_matrix} - \bs{B}_{k, \bs{\xi}} \left(x \right) = - \begin{pmatrix} - B_{1, k, \bs{\xi}}\left(x_1\right) & B_{2, k, \bs{\xi}}\left(x_1\right) & \cdots & B_{k+l, k, \bs{\xi}}\left(x_1\right) \\ - B_{1, k, \bs{\xi}}\left(x_2\right) & B_{2, k, \bs{\xi}}\left(x_2\right) & \cdots & B_{k+l, k, \bs{\xi}}\left(x_2\right) \\ - \vdots & \vdots & \ddots & \vdots \\ - B_{1, k, \bs{\xi}}\left(x_{\card{\bs{x}}}\right) & B_{2, k, \bs{\xi}}\left(x_{\card{\bs{x}}}\right) & \cdots & B_{k+l, k, \bs{\xi}}\left(x_{\card{\bs{x}}}\right) \\ - \end{pmatrix} -\end{equation} -where $B_{j, k, \bs{\xi}} \left( \bs{x} \right)$ is the $j^{th}$ B-spline defined by de Boor's recursive algorithm -% -\begin{equation} - \label{eq:deboor_recursion} - B_{j, k, \bs{\xi}} \left(x\right) = - \omega_{j, k, \bs{\xi}}\left(x\right) B_{j, k-1, \bs{\xi}}\left(x\right) + \left(1 - - \omega_{j+1, k, \bs{\xi}} \right) B_{j+1, k-1, \bs{\xi}}\left(x\right), -\end{equation} -with -\begin{equation} - \label{eq:recurrence_relation_null} - B_{j, k, \bs{\xi}}\left(x\right) = 0 \quad \text{for } x \notin - \left[\xi_{j}, \xi_{j+k}\right), \quad - B_{j, 1, \bs{\xi}}\left(x\right) = - \begin{cases} - 1 & x \in \left[\xi_{j}, \xi_{j+1} \right)\\ - 0 & \text{otherwise}, - \end{cases} -\end{equation} -and -\begin{equation} - \label{eq:omega} - \omega_{jk}\left(x\right) = - \begin{cases} - \hfil 0 & \hfil x \leq \xi_{j} \\ - \frac{x - \xi_{j}} {\xi_{j + k - 1} - \xi_{j}} & \xi_{j} < x < \xi_{j+ k - 1} \\ - \hfil 1 & \hfil \xi_{j + k - 1} \leq x - \end{cases}. -\end{equation} -The basis is a partition of unity over the support of $\bs{x},$ -\ie, all rows sum to one, and the span of the basis defines a spline -space $\mathbb{S}_{k, \bs{\xi}} = \text{span } \bs{B}_{k, \bs{\xi}}.$ - -\subsection{Building B-spline basis} -Build B-spline basis via \code{bsplines}. -<<>>= -str(bsplines) -@ -<>= -stopifnot(length(formals(bsplines)) == 5L) -# if the above fails verify that the description of the arguments below are -# correct. -@ -\noindent There are \Sexpr{length(formals(bsplines))} arguments for the -\code{bsplines} call. -\begin{itemize} - \item \code{x} is a numeric vector as noted in the above notation. - \item \code{iknots} is a numeric vector of locations for the internal knots. - \item \code{df} an integer value for the number of degrees of freedom. - \code{df} will be equal to the \code{order + length(iknots)}. - \item \code{bknots} boundary knots. This needs only be a numeric vector of - length two. The \code{bsplines} call will generate the full knot sequence - $\bs{\xi}$ with the proper multiplicity of the boundary knots. - \item \code{order} the polynomial order. By default \code{order = 4L} - generates fourth-order B-splines, also called cubic B-splines as polynomial - order the polynomial degree plus one. -\end{itemize} -You only need to specify \code{iknots} or \code{df}. If both are provided a -warning will be given and \code{iknots} will be used in the construction of the -basis matrix. - -<>= -x <- seq(0, 6, length = 500) # will be used again -bmat <- bsplines(x = x, - iknots = c(1.0, 1.5, 2.3, 4.0, 4.5)) -bmat -@ - -\subsection{cpr::bsplines vs splines::bs} -The standard installation of \proglang{R} includes the -\pkg{splines} package and the \code{splines::bs} function for -generating the basis matrix of B-splines, \ie, the matrix shown in -\eqref{eq:basis_matrix}. The provided \pkg{cpr::bsplines} function -for generating B-spline basis matrices with the class -\code{cpr\_bs} has arguments and defaults which are beneficial, and required, -for the use in the control polygon reduction algorithm. -The differences in the functional arguments and the attributes -of the return objects between \code{splines::bs} and \code{cpr::bsplines} are -listed in Table~\ref{tab:bs_vs_bsplines}. -\begin{table} - \caption{Comparison of the arguments, with default values, for \code{splines::bs} and - \code{cpr::bsplines}. The attributes for the resulting \code{splines::bs} and - \code{cpr\_bs} objects are also reported.} - \label{tab:bs_vs_bsplines} - \centering -<>= -stopifnot( - names(attributes(bmat)) == - c("dim", "order", "iknots", "bknots", "xi", "xi_star", "class", "call", "environment")) -@ - -\begin{tabular}{lll} \hline - & \code{splines::bs} & \code{cpr::bsplines} \\ \hline - Arguments & & \\ - & \code{x} & \code{x} \\ - & \code{df} & \code{df} \\ - & \code{knots} & \code{iknots} \\ - & \code{degree = 3} & \code{order = 4L} \\ - & \code{Boundary.knots = range(x)} & \code{bknots = range(x)} \\ - & \code{intercept = FALSE} & -- \\ \hline -Attributes & & \\ - & \code{dim} & \code{dim} \\ - & \code{degree} & \code{order} \\ - & \code{knots} & \code{iknots} \\ - & \code{Boundary.knots} & \code{bknots} \\ - & \code{intercept} & -- \\ - & -- & \code{xi} \\ - & -- & \code{xi\_star} \\ - & \code{class} & \code{class} \\ - & -- & \code{call} \\ - & -- & \code{environment} \\ \hline -\end{tabular} -\end{table} - -A major difference between the two functions is related to the \code{intercept} -argument of \code{splines::bs}. By default, \code{splines::bs} will omit the -first column of the basis whereas \code{cpr::bsplines} will return the whole -basis. The omission of the first column of the basis generated by -\code{splines::bs} allows for additive \code{splines::bs} calls to be used on -the right-hand-side of a regression formula and generate a full rank design -matrix. If additive \code{cpr::bsplines} calls, or additive \code{splines::bs} -with \code{intercept = TRUE}, are on the right-hand-side of the regression -equation the resulting design matrix will be rank deficient. This is a result -of the B-splines being a partition of unity. As the CPR algorithm is based on -having the whole basis, {he \code{cpr::bsplines} function is provided to make it -easy to work with the whole basis without having to remember to use non-default -settings in \code{splines::bs}. The default call \code{splines::bs(x)} is -replicated by \code{cpr::bsplines(x)[, -1]} and the default call -\code{cpr::bsplines(x)} is replicated \code{splines::bs(x, intercept = TRUE)}. - -Specifying the polynomial order and knot sequence between the two functions -differ between \code{splines::bs} and \code{cpr::bsplines}. \code{splines::bs} -uses the polynomial \code{degree} whereas \code{cpr::bsplines} uses the -polynomial \code{order} (order = degree + 1) to define the splines. The default -for both \code{splines::bs} and \code{cpr::bsplines} is to generate cubic -B-splines. - -The call to \code{splines::bs} which will generate the same basis matrix as -\code{bmat} above is: -<<>>= -bmat_bs <- splines::bs(x = seq(0, 6, length = 500), - knots = c(1.0, 1.5, 2.3, 4.0, 4.5), - intercept = TRUE) -all.equal(unclass(bmat), unclass(bmat_bs), check.attributes = FALSE) -@ - -<>= -stopifnot(all.equal(unclass(bmat), unclass(bmat_bs), check.attributes = FALSE)) -@ - -For both \code{splines::bs} and \code{cpr::bsplines} only the degrees of freedom -or the internal knots need to be specified. If the end user specifies both, the -specified knots take precedence. If only \code{df} is specified then \code{df - -order} internal knots will be generated. - -There is a very important difference in way \code{splines::bs} and -\code{cpr::bsplines} build the sequence of internal knots. \code{splines::bs} -uses quantiles of the whole numeric vector $\bs{x}$ whereas \code{cpr::bsplines} -uses quantiles of the unique values of $\bs{x} \backslash \{\min(\bs{x}), -\max(\bs{x}) \}.$ -For example: the internal knots for a B-spline with 9 degrees of freedom on a -forth order spline with an intercept, \code{splines::bs} will generate a -sequence of internal knots via a call equivalent to -<<>>= -y <- sample(runif(100, 0, 6), 500, replace = TRUE) -# k <- degrees of freedom - order + (1L - intercept) -degf <- 9 -ordr <- 4 -k <- degf - ordr + (1 - as.integer(TRUE)) -p <- seq(0, 1, length = k + 2L)[-c(1, k + 2L)] -stats::quantile(y, probs = p) -@ -\noindent whereas \code{cpr::bsplines} will generate a sequence equivalent to -<<>>= -stats::quantile(unique(y)[-c(1, length(unique(y)))], - probs = seq(1, degf - ordr, by = 1) / (degf - ordr + 1)) -@ -\noindent The function \code{cpr::trimmed_quantile} is provided to generate such -sequences. -<<>>= -trimmed_quantile(y, probs = seq(1, degf - ordr, by = 1) / (degf - ordr + 1)) -@ - -The return object from both \code{splines::bs} and \code{cpr::bsplines} is a -matrix. The attributes returned include the argument values used to construct -the basis. The major difference in the attributes between the two objects is -that \code{cpr::bsplines} returns the full knot sequence, $\bs{\xi},$ in the -\code{xi} element and the Greville sites, $\bs{\xi}^{*},$ in the \code{xi\_star} -element which is useful for building control polygons. Lastly, the classes for -the two objects differ: \code{splines::bs} returns a three classes, -\code{c("bs", "basis", "matrix")} and \code{cpr::bsplines} returns two classes, -\code{c("cpr\_bs", "matrix")}. - -\subsection{Plotting B-spline Basis} -The \pkg{cpr} package provides an plotting -method, built on top of \pkg{ggplot2}~\citep{R-ggplot2}, -for \code{cpr\_bs} objects. -<