Dev - parallelize influence_of_iknots

.Rbuildignore

@@ -17,3 +17,5 @@
 ^_pkgdown\.yml$
 ^docs$
 ^pkgdown$
+^lib
+^docs

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: cpr
 Title: Control Polygon Reduction
-Version: 0.3.0.9000
+Version: 0.3.0.9005
 Authors@R:
     c(person("Peter", "DeWitt", email = "dewittpe@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-6391-0795")),
       person("Samantha", "MaWhinney", email = "sam.mawhinney@ucdenver.edu", role = c("ths")),

NAMESPACE

@@ -28,6 +28,7 @@ S3method(plot,cpr_cnr)
 S3method(plot,cpr_cp)
 S3method(plot,cpr_cpr)
 S3method(plot,cpr_influence_of_iknots)
+S3method(plot,cpr_summary_cpr_cpr)
 S3method(predict,cpr_cn)
 S3method(predict,cpr_cp)
 S3method(print,cpr_bs)

NEWS.md

@@ -2,6 +2,12 @@
 
 ## New Features
 
+* `cpr`'s `progress` argument has been extended to control if a progress bar is
+  used for just the cpr steps, or if a more detailed progress for of the
+  influence weight calculations is reported.
+
+* `influence_of_iknots` gains parallel execution via `pbapply` (#17)
+
 * `plot.cpr_cp` gains the argument `comparitive` which, when set to `FALSE` and
   only one `cpr_cp` is passed in for plotting, the graphic will appear more like
   the `plot.cpr_bs` results.  When `comparitive = TRUE` or more than one

R/RcppExports.R

@@ -25,8 +25,12 @@ coarsen_theta <- function(j, xi, k, theta) {
     .Call('_cpr_coarsen_theta', PACKAGE = 'cpr', j, xi, k, theta)
 }
 
-hat_theta <- function(j, xi, k, theta, calculate_F, Sigma) {
-    .Call('_cpr_hat_theta', PACKAGE = 'cpr', j, xi, k, theta, calculate_F, Sigma)
+hat_theta <- function(j, xi, k, theta) {
+    .Call('_cpr_hat_theta', PACKAGE = 'cpr', j, xi, k, theta)
+}
+
+test_statistic <- function(j, xi, k, theta, Sigma) {
+    .Call('_cpr_test_statistic', PACKAGE = 'cpr', j, xi, k, theta, Sigma)
 }
 
 #' Rank of a Matrix

R/bsplineD.R

@@ -58,11 +58,25 @@
 #' points(x, bmat %*% theta, col = 'blue')
 #' grid()
 #'
-#' plot(x, fprime(x), type = "l", main = bquote(frac(d,dx)~f(x)), ylab = "", xlab = "")
+#' plot(
+#'     x
+#'   , fprime(x)
+#'   , type = "l"
+#'   , main = bquote(frac(d,dx)~f(x))
+#'   , ylab = ""
+#'   , xlab = ""
+#' )
 #' points(x, bmatD1 %*% theta, col = 'blue')
 #' grid()
 #'
-#' plot(x, fdoubleprime(x), type = "l", main = bquote(frac(d^2,dx^2)~f(x)), ylab = "", xlab = "")
+#' plot(
+#'     x
+#'   , fdoubleprime(x)
+#'   , type = "l"
+#'   , main = bquote(frac(d^2,dx^2)~f(x))
+#'   , ylab = ""
+#'   , xlab = ""
+#' )
 #' points(x, bmatD2 %*% theta, col = 'blue')
 #' grid()
 #'
@@ -76,8 +90,8 @@
 #' iknots <- sort(runif(rpois(1, 3), 1, 9))
 #' bknots <- c(0, 10)
 #'
-#' # basis matrix and the first and second derivatives thereof, for cubic (order =
-#' # 4) b-splines
+#' # basis matrix and the first and second derivatives thereof, for cubic
+#' # (order = 4) b-splines
 #' bmat  <- bsplines(xvec, iknots, bknots = bknots)
 #' bmat1 <- bsplineD(xvec, iknots, bknots = bknots, derivative = 1)
 #' bmat2 <- bsplineD(xvec, iknots, bknots = bknots, derivative = 2)
@@ -87,9 +101,11 @@
 #'
 #' # plot data
 #' plot_data <-
-#'   data.frame(Spline            = as.numeric(bmat %*% theta),
-#'              First_Derivative  = as.numeric(bmat1 %*% theta),
-#'              Second_Derivative = as.numeric(bmat2 %*% theta))
+#'   data.frame(
+#'       Spline            = as.numeric(bmat %*% theta)
+#'     , First_Derivative  = as.numeric(bmat1 %*% theta)
+#'     , Second_Derivative = as.numeric(bmat2 %*% theta)
+#'   )
 #' plot_data <- stack(plot_data)
 #' plot_data <- cbind(plot_data, data.frame(x = xvec))
 #'

R/bsplines.R

@@ -3,7 +3,8 @@
 #' An implementation of Carl de Boor's recursive algorithm for building
 #' B-splines.
 #'
-#' There are several differences between this function and \code{\link[splines]{bs}}.
+#' There are several differences between this function and
+#' \code{\link[splines]{bs}}.
 #'
 #' The most important difference is how the two methods treat the right-hand end
 #' of the support.  \code{\link[splines]{bs}} uses a pivot method to allow for
@@ -22,7 +23,8 @@
 #' @references
 #' C. de Boor, "A practical guide to splines. Revised Edition," Springer, 2001.
 #'
-#' H. Prautzsch, W. Boehm, M. Paluszny, "Bezier and B-spline Techniques," Springer, 2002.
+#' H. Prautzsch, W. Boehm, M. Paluszny, "Bezier and B-spline Techniques,"
+#' Springer, 2002.
 #'
 #' @param x a numeric vector
 #' @param iknots internal knots
@@ -59,8 +61,12 @@
 #' # The x-axis, by default, show the knot locations.  Other options are numeric
 #' # values, and/or to use a second x-axis
 #'
-#' plot(bmat, show_xi = TRUE,  show_x = FALSE) # default, knot, symbols, on lower axis
-#' plot(bmat, show_xi = FALSE, show_x = TRUE)  # Numeric value for the knot locations
+#' plot(bmat, show_xi = TRUE,  show_x = FALSE) # default, knot, symbols, on lower
+#'                                             # axis
+#'
+#' plot(bmat, show_xi = FALSE, show_x = TRUE)  # Numeric value for the knot
+#'                                             # locations
+#'
 #' plot(bmat, show_xi = TRUE,  show_x = TRUE)  # symbols on bottom, numbers on top
 #'
 #' # quadratic splines
@@ -92,6 +98,7 @@ bsplines <- function(x, iknots = NULL, df = NULL, bknots = range(x), order = 4L)
   iknots <- iknots_or_df(x, iknots, df, order)
 
   rtn <- cpp_bsplines(x = x, iknots = iknots, bknots = bknots, order = order)
+  #colnames(rtn) <- paste0("bspline", seq(1, ncol(rtn)))
   attr(rtn, "call") <- match.call()
   attr(rtn, "environment") <- parent.frame()
   class(rtn) <- c("cpr_bs", "matrix")

R/btensor.R

@@ -111,6 +111,10 @@ btensor <- function(x, df = NULL, iknots = NULL, bknots, order) {
 #' @method print cpr_bt
 #' @export
 print.cpr_bt <- function(x, ...) {
-  cat("Tensor Product Matrix dims: [", paste(format(dim(x), big.mark = ",", trim = TRUE), collapse = " x "), "]\n\n", sep = "")
+  cat("Tensor Product Matrix dims: ["
+      , paste(format(dim(x), big.mark = ",", trim = TRUE), collapse = " x ")
+      , "]\n\n"
+      , sep = ""
+  )
   invisible(x)
 }

R/cn.R

@@ -25,18 +25,18 @@
 #'  \item{rse}{the residual standard error for the regression models}
 #'  }
 #'
-#' @seealso \code{\link{summary.cpr_cn}}, \code{\link{cnr}}
+#' @seealso \code{\link{summary.cpr_cn}}, \code{\link{cnr}},
+#' \code{\link{plot.cpr_cn}} for plotting control nets
 #'
 #' @examples
 #'
-#' acn <- cn(log10(pdg) ~ btensor(  x = list(day, age)
-#'                                 , df = list(30, 4)
-#'                                 , bknots = list(c(-1, 1), c(44, 53)))
+#' acn <- cn(log10(pdg) ~
+#'               btensor(   x     = list(day, age)
+#'                       , df     = list(30, 4)
+#'                       , bknots = list(c(-1, 1), c(44, 53))
+#'               )
 #'            , data = spdg)
-#'
-#' # plot3D
-#' plot(acn, rgl = FALSE)
-#'
+#' str(acn, max.level = 1)
 #'
 #' @export
 cn <- function(x, ...) {
@@ -45,7 +45,8 @@ cn <- function(x, ...) {
 
 #' @export
 #' @rdname cn
-#' @param theta a vector of (regression) coefficients, the ordinates of the control net.
+#' @param theta a vector of (regression) coefficients, the ordinates of the
+#' control net.
 cn.cpr_bt <- function(x, theta, ...) {
   xi_stars <- lapply(attr(x, "bspline_list"), attr, which = "xi_star")
 
@@ -70,16 +71,22 @@ cn.cpr_bt <- function(x, theta, ...) {
 #' @rdname cn
 #' @param formula a formula that is appropriate for regression method being used.
 #' @param data a required \code{data.frame}
-#' @param method the regression method such as \code{\link[stats]{lm}}, \code{\link[stats]{glm}}, \code{\link[lme4]{lmer}}, etc.
-#' @param method.args a list of additional arguments to pass to the regression method.
-#' @param keep_fit (logical, defaults to \code{FALSE}).  If \code{TRUE} the regression model fit is retained and returned in the the \code{fit} element. If \code{FALSE} the regression model is not saved and the \code{fit} element will be \code{NA}.
-#' @param check_rank (logical, defaults to \code{TRUE}) if TRUE check that the design matrix is full rank.
-cn.formula <- function(formula, data, method = stats::lm, method.args = list(),  keep_fit = TRUE, check_rank = TRUE, ...) {
+#' @param method the regression method such as \code{\link[stats]{lm}},
+#' \code{\link[stats]{glm}}, \code{\link[lme4]{lmer}}, etc.
+#' @param method.args a list of additional arguments to pass to the regression
+#' method.
+#' @param keep_fit (logical, defaults to \code{FALSE}).  If \code{TRUE} the
+#' regression model fit is retained and returned in the the \code{fit} element.
+#' If \code{FALSE} the regression model is not saved and the \code{fit} element
+#' will be \code{NA}.
+#' @param check_rank (logical, defaults to \code{TRUE}) if TRUE check that the
+#' design matrix is full rank.
+cn.formula <- function(formula, data, method = stats::lm, method.args = list(), keep_fit = TRUE, check_rank = TRUE, ...) {
+
   # check for some formula specification issues
-  fterms <- stats::terms(formula)
-  fterms
-  if (sum(grepl("btensor", attr(fterms, "term.labels"))) != 1) {
-    stop("btensor() must appear once, with no effect modifiers, on the right hand side of the formula.")
+  rhs_check <- grepl("btensor", attr(stats::terms(formula), "term.labels"))
+  if ( !rhs_check[1] | any(rhs_check[-1]) ) {
+    stop("btensor() must appear first, once, and with no effect modifiers, on the right hand side of the formula.")
   }
 
   # this function will add f_for_use and data_for_use into this environment
@@ -91,7 +98,8 @@ cn.formula <- function(formula, data, method = stats::lm, method.args = list(),
   cl <- c(cl, method.args)
 
   fit <- do.call(regression, cl)
-  COEF_VCOV <- coef_vcov(fit)
+  Bmat <- stats::model.frame(fit)[[2]]
+  COEF_VCOV <- coef_vcov(fit, theta_idx = seq(1, ncol(Bmat), by = 1))
 
   if (check_rank) {
     m <- stats::model.matrix(lme4::nobars(f_for_use), data_for_use)
@@ -105,8 +113,6 @@ cn.formula <- function(formula, data, method = stats::lm, method.args = list(),
   cl[[1]] <- as.name("cn")
   cl <- as.call(cl)
 
-  Bmat <- eval(extract_cpr_bsplines(formula), data, environment(formula))
-
   out <- cn.cpr_bt(Bmat, COEF_VCOV$theta)
 
   # update elements of the cpr_bs object
@@ -135,20 +141,3 @@ print.cpr_cn <- function(x, ...) {
   print(x$cn, ...)
   invisible(x)
 }
-
-
-extract_cpr_bsplines <- function(form) {
-  B <- NULL
-  rr <- function(x) {
-    if (is.call(x) && grepl("bsplines|btensor", deparse(x[[1]]))) {
-      B <<- x
-    } else if (is.recursive(x)) {
-      as.call(lapply(as.list(x), rr))
-    } else {
-      x
-    }
-  }
-
-  z <- lapply(as.list(form), rr)
-  B
-}

R/cnr.R

@@ -16,7 +16,7 @@
 #' \code{\link{influence_weights}}.
 #' @param n_polycoef the number of polynomial coefficients to use when assessing
 #' the influence of each internal knot.
-#' @param progress show a progress bar.
+#' @param progress controls the level of progress messaging.
 #' @param ... not currently used
 #'
 #' @return A \code{cpr_cnr} object.  This is a list of \code{cpr_cn} objects.
@@ -40,24 +40,26 @@
 #' plot(cnr0)
 #'
 #' @export
-cnr <- function(x, margin, n_polycoef = 20L, progress = interactive(), ...) {
+cnr <- function(x, margin, n_polycoef = 20L, progress = c('cnr', 'influence', 'none'), ...) {
   UseMethod("cnr")
 }
 
 #' @export
-cnr.cpr_cn <- function(x, margin = seq_along(x$bspline_list), n_polycoef = 20L, progress = interactive(), ...) {
+cnr.cpr_cn <- function(x, margin = seq_along(x$bspline_list), n_polycoef = 20L, progress = c('cnr', 'influence', 'none'), ...) {
+
+  progress <- match.arg(progress, several.ok = FALSE)
 
   out <- vector("list", length = sum(sapply(lapply(x$bspline_list[margin], attr, which = "iknots"), length)) + 1L)
 
-  if (progress) {
+  if (progress == 'cnr') {
     pb <- utils::txtProgressBar(max = length(out), style = 3) # nocov
     prg <- 0 # nocov
     utils::setTxtProgressBar(pb, prg) # nocov
   }
 
   for(i in rev(seq_along(out)[-1])) {
     out[[i]] <- x
-    w <- summary(influence_of_iknots(out[[i]], margin, n_polycoef))
+    w <- summary(influence_of_iknots(out[[i]], margin = margin, n_polycoef = n_polycoef, verbose = (progress == 'influence'), ...))
     w <- w[w$influence_rank > 1, ]
     nkts <- lapply(split(w, f = w$margin), getElement, "iknot")
 
@@ -66,16 +68,25 @@ cnr.cpr_cn <- function(x, margin = seq_along(x$bspline_list), n_polycoef = 20L,
       nkts <- nkts[sort(names(nkts))]
     }
 
-    x <- eval(stats::update(x, formula = newknots(x$call$formula, nkts), keep_fit = TRUE, check_rank = FALSE, evaluate = FALSE), parent.frame())
+    x <-
+      eval(
+        stats::update(
+             x
+           , formula = newknots(x$call$formula, nkts)
+           , keep_fit = TRUE
+           , check_rank = FALSE
+           , evaluate = FALSE)
+        , parent.frame()
+      )
 
-    if (progress) {
+    if (progress == 'cnr') {
       utils::setTxtProgressBar(pb, prg <- prg + 1) # nocov
     }
   }
 
   out[[1]] <- x
 
-  if (progress) {
+  if (progress == 'cnr') {
     utils::setTxtProgressBar(pb, prg <- prg + 1) # nocov
     close(pb) # nocov
   }