#' @importFrom parallel mclapply #' step_optim <- function(model, data, kseq = NA, prm=list(NA), direction = c("backward","forward","backwardboth","forwardboth"), optimfun = rls_optim, scorefun = rmse, ...){ # Do: # - change all lapply # - Maybe have "cloneit" argument in optimfun, then don't clone inside optim. # - Add argument controlling how much is kept in each iteration (e.g all fitted models) # # - Help: prm <- list(I__degree = c(min=1, max=7), mu_tday__nharmonics = c(min=1, max=7)) # - help: It's not checked that it's the score is calculated on the same values! WARNING should be printed if some models don't forecast same points # # First direction is default if(length(direction) > 1){ direction <- direction[1] } # Don't change the given mfull <- model$clone_deep() # Insert the starting prm values if(!is.na(prm[1])){ if( direction == "backward" ){ mfull$insert_prm(unlist(getse(prm, "max"))) }else if( direction == "forward" ){ mfull$insert_prm(unlist(getse(prm, "min"))) }else{ # Both directions, then start at init, or halfway between min and max i <- which(sapply(prm, function(x){ "init" %in% names(x) })) if(length(i)){ mfull$insert_prm(unlist(getse(prm[i], "init"))) } i <- which(sapply(prm, function(x){ !"init" %in% names(x) })) if(length(i)){ mfull$insert_prm(round(unlist(getse(prm[i], "max")) - unlist(getse(prm[i], "min")) / 2)) } } } # For keeping all the results L <- list() # m <- mfull$clone_deep() if(length(grep("backward",direction))){ # Optimize from the full model res <- optimfun(m, data, kseq, printout=TRUE, ...) # Keep it istep <- 1 L[[istep]] <- list(model = m$clone_deep(), result = res) }else{ # Optimize from the null model m$inputs <- list() # Must be set istep <- 0 res <- list(value=Inf, par=m$get_prmbounds("init")) } # Helper c_mStep <- function(l, nms){ names(l) <- nms l <- l[!is.na(l)] c(mStep, l) } # Go done <- FALSE while(!done){ # Next step istep <- istep + 1 # Insert the optimized parameters from last step m$prmbounds[names(res$par),"init"] <- res$par # message("------------------------------------") message("Reference score value: ",res$value) # -------- mStep <- list() # Generate the input modified models if(length(grep("backward|both", direction))){ # Remove input from the current model one by one if(length(m$inputs) > 1){ tmp <- mclapply(1:length(m$inputs), function(i){ ms <- m$clone_deep() # Remove one input ms$inputs[[i]] <- NULL return(ms) }) mStep <- c_mStep(tmp, pst("-",names(m$inputs))) } } if(length(grep("forward|both", direction))){ # Add input one by one iin <- which(!names(mfull$inputs) %in% names(m$inputs)) if(length(iin)){ tmp <- mclapply(iin, function(i){ ms <- m$clone_deep() # Add one input ms$inputs[[length(ms$inputs) + 1]] <- mfull$inputs[[i]] names(ms$inputs)[length(ms$inputs)] <- names(mfull$inputs)[i] return(ms) }) mStep <- c_mStep(tmp, pst("+",names(mfull$inputs)[iin])) } } # Step parameters if(!is.na(prm[1])){ if(length(grep("backward|both", direction))){ # Count down the parameters one by one tmp <- mclapply(1:length(prm), function(i){ p <- m$get_prmvalues(names(prm[i])) # If the input is not in the current model, then p is NA, so don't include it for fitting if(!is.na(p)){ # Only the ones with prms above minimum if(prm[[i]]["min"] < p){ ms <- m$clone_deep() p <- p - 1 ms$insert_prm(p) # Return both the prm value and the name of the model to be printed return(list(ms, pst(names(p),"=",p))) } } return(list(NA,NA)) }) mStep <- c_mStep(getse(tmp,1), getse(tmp,2)) } if(length(grep("forward|both", direction))){ # Count up the parameters one by one tmp <- mclapply(1:length(prm), function(i){ p <- m$get_prmvalues(names(prm[i])) # If the input is not in the current model, then p is NA, so don't include it for fitting if(!is.na(p)){ # Only the ones with prms above minimum if(p < prm[[i]]["max"]){ ms <- m$clone_deep() p <- p + 1 ms$insert_prm(p) # Return both the prm value and the name of the model to be printed return(list(ms, pst(names(p),"=",p))) } } return(list(NA,NA)) }) mStep <- c_mStep(getse(tmp,1), getse(tmp,2)) } } # Optimize all the step models resStep <- mclapply(1:length(mStep), function(i, ...){ res <- try(optimfun(mStep[[i]], data, kseq, printout=FALSE, ...)) if(class(res) == "try-error"){ browser() } message(names(mStep)[[i]], ": ", res$value) return(res) }, ...) names(resStep) <- names(mStep) # Is one the step models score smaller than the current ref? valStep <- unlist(getse(resStep, "value")) imin <- which.min(valStep) if( valStep[imin] < res$value ){ # Keep the best model m <- mStep[[imin]] res <- resStep[[imin]] # Keep for the result L[[istep]] <- list(model = m$clone_deep(), result = resStep[[imin]]) }else{ # No improvement obtained from reduction, so return the current model (last in the list) message("------------------------------------\n\nDone") message(print(m)) done <- TRUE } } invisible(L) }