R/predict.DDPstar.regfun.R

In DDPstar: Density Regression via Dirichlet Process Mixtures of Normal Structured Additive Regression Models

predict_DDPstar.regfun <-
function(object, newdata, select = NULL, parallel = c("no", "multicore", "snow"), ncpus = 1, cl = NULL) {
	doMCMC <- function(k, object, Xp, ind.coeff) {
		ncomp <- ncol(object$fit$probs)
		if(ncomp == 1) {
			Beta <- matrix(object$fit$beta[k,,ind.coeff], nrow = 1)
		} else {
			Beta <- object$fit$beta[k,,ind.coeff]
		}
		regfun <- colSums(object$fit$probs[k,]*t(Xp%*%t(Beta)))
		res <- list()
		res$regfun <- regfun
		res$foo <- 1
		res	
	}

	# Select does not include the intercept!
	if(is.null(select)) {
		select.fun <- 1:object$fit$mm$iformula$npartial
	} else {
		if(!all(select %in% 1:object$fit$mm$iformula$npartial)) {
			stop("Some of the selected components are not in the formula")
		}
		select.fun <- select
	}
	
	# Dimension of each component
	# Number of subcomponents per component
	np.e.sc <- cumsum(object$fit$mm$K.length)
	np.s.sc <- np.e.sc - object$fit$mm$K.length + 1 

	dim.terms <- unlist(lapply(object$fit$mm$K, nrow))
	np.e <- cumsum(dim.terms)
	np.s <- np.e - dim.terms + 1 

	Xp <- predict_design.matrix.DDPstar.aux(object = object$fit$mm, newdata = newdata, select = select.fun)
	
	# Which coeff are needed?
	ind.coeff <- NULL
	# Add the intecept (if needed)
	if(is.null(select)) {
		Xp <- cbind(1, Xp)
		ind.coeff <- c(ind.coeff, 1)
	}
	for(i in select.fun) {
		ind.coeff <- c(ind.coeff, np.s[np.s.sc[i+1]]:np.e[np.e.sc[i+1]])	
	}
	
	parallel <- match.arg(parallel)	
	nsim <- nrow(object$fit$sd)
	if(nsim > 0) {
		do_mc <- do_snow <- FALSE
	    if (parallel != "no" && ncpus > 1L) {
	        if (parallel == "multicore") {
	            do_mc <- .Platform$OS.type != "windows"
	        } else if (parallel == "snow") {
	            do_snow <- TRUE
	        }
	        if (!do_mc && !do_snow) {
	            ncpus <- 1L
	        }       
	        loadNamespace("parallel") # get this out of the way before recording seed
	    }

	    resMCMC <- if (ncpus > 1L && (do_mc || do_snow)) {
                if (do_mc) {
                    parallel::mclapply(seq_len(nsim), doMCMC, object = object, Xp = Xp, ind.coeff = ind.coeff, mc.cores = ncpus)
                } else if (do_snow) {                
                    if (is.null(cl)) {
                        cl <- parallel::makePSOCKcluster(rep("localhost", ncpus))
                        if(RNGkind()[1L] == "L'Ecuyer-CMRG") {
                            parallel::clusterSetRNGStream(cl)
                        }
                        res <- parallel::parLapply(cl, seq_len(nsim), doMCMC, object = object, Xp = Xp, ind.coeff = ind.coeff)
                        parallel::stopCluster(cl)
                        res
                    } else {
                        if(!inherits(cl, "cluster")) {
                            stop("Class of object 'cl' is not correct")
                        } else {
                            parallel::parLapply(cl, seq_len(nsim), doMCMC, object = object, Xp = Xp, ind.coeff = ind.coeff)
                        }                        
                    }
                }
            } else {
                lapply(seq_len(nsim), doMCMC, object = object, Xp = Xp, ind.coeff = ind.coeff)
            }

        resMCMC <- simplify2array(resMCMC)
        regfun <- simplify2array(resMCMC["regfun",])

	} else {
		stop("nsave should be larger than zero.")
	}

	res <- list()
	res$Xp <- Xp
	res$regfun <- regfun
	res$select <- select
	class(res) <- "DDPstar.regfun"
	res
}