update_rq.R

#' @title update_rq
#'
#' @description This function plots the quantiles
#' @param Lfor List of the data needed to plot
#' @export
#' @examples
#' plotprob()


update_rq <- function(Xny, tau, k,  model = NULL, debug = TRUE){
    if(is.null(model)) stop("Model needs to be initialized before update")

    N <- nrow(Xny)

    res <- model$info[[paste0("q", tau)]][[paste0("k", k)]]
    ## Maybe this is no correct - we should put weights on this!
    Xold <- res$X

    
    # Prepare for keeping for the parameter estimates
    BETA <- matrix(as.numeric(NA), nrow = N, ncol = model$K)
    i <- 1
    counter <- 1

    Ypred <- matrix(as.numeric(NA), nrow = N, ncol = 1)

    if(debug) {
    	print("yeah")
    	list_idx <- matrix(as.numeric(NA), nrow = 1, ncol = length(res$Ih)+1)
    	list_idx[1,] <- c(0,res$Ih)
    }

    while(counter <= N){
        #cat("counter", counter, "\n")
        

        beta <- res$xB[1:model$K]
        BETA[counter,] <- t(beta)

        j <- 0
        res <- rq_update_cpp(newX = matrix(Xny[counter,], nrow = 1),
                             X = Xold,
                             IX = model$IX,
                             Iy = model$Iy,
                             Ih =  res$Ih,
                             Ihc = res$Ihc,
                             beta = beta,
                             Rny = res$R,
                             K = model$K, 
                             n = nrow(Xold),
                             xB = res$xB,
                             P = res$P,
                             tau = tau,
                             k = 0,
                             kk = 0,
                             i = i,
                             n_in_bin = model$n_in_bin,
                             W = model$regprm$lambda)
        Ypred[counter] <- res$Ypred
        i = res$i
        
        Xold <- res$X
        
        res$IH <- res$Ih

        #cat("Here", "\n")
        #if(counter == 119) browser()


    	if(debug) {
    		list_idx <- rbind(list_idx, c(counter, res$Ih))
    	}


        resAlg <- rq_simplex_cpp(X = Xold[,1:model$K],
                                 Ih =  res$Ih,
                                 Ihc = res$Ihc,
                                 IH = as.matrix(res$Ih),
                                 K = model$K,
                                 n = nrow(Xold),
                                 xB = res$xB,
                                 P = res$P,
                                 tau = tau)

        res$CON <- resAlg$CON
        res$s <- resAlg$s
        res$g <- resAlg$g
        res$gain <- resAlg$gain
        res$md <- resAlg$md
        res$alpha <- resAlg$alpha
        res$h <- resAlg$h
        res$IH <- resAlg$IH
        res$cq <- resAlg$cq
        res$q <- resAlg$q

        if(length(res$md) == 0) res$md <- 1
        
        while(res$gain <= 0 & res$md < 0 & j < 24 & res$CON < 10^6){
            
            j <- j + 1
            
            z = res$xB - as.numeric(res$alpha) * res$h
            
            IhM <- res$Ih[res$s + 1]
            IhcM <- res$Ihc[res$q + 1]
            res$Ih[res$s + 1] <- IhcM
            
            res$Ihc[res$q+1] <- IhM
            res$P[res$q+1] <- res$cq
            
            res$xB <- z
            res$xB[res$q + 1 + model$K] <- res$alpha
            
            
            dummy <- sortIdx(res$Ih)
            res$Ih <- dummy$newVal; IndexIh <- dummy$idxVal
            dummy <- sortIdx(res$Ihc)
            res$Ihc <- dummy$newVal; IndexIhc <- dummy$idxVal
            
            if(debug) {
    			list_idx <- rbind(list_idx, c(counter, res$Ih))
    		}
            
            res$P <- res$P[IndexIhc]
            xBm <- res$xB[((model$K+1):length(res$xB))]
            xBm <- xBm[IndexIhc]
            res$xB[((model$K+1):length(res$xB))] <- xBm

            #cat("Here2", "\n")            
            resAlg <- rq_simplex_cpp(X = Xold[,1:model$K],
                                     Ih =  res$Ih,
                                     Ihc = res$Ihc,
                                     IH = as.matrix(res$Ih),
                                     K = model$K,
                                     n = nrow(Xold),
                                     xB = res$xB,
                                     P = res$P,
                                     tau = tau)

            res$CON <- resAlg$CON
            res$s <- resAlg$s
            res$g <- resAlg$g
            res$gain <- resAlg$gain
            res$md <- resAlg$md
            res$alpha <- resAlg$alpha
            res$h <- resAlg$h
            res$IH <- resAlg$IH
            res$cq <- resAlg$cq
            res$q <- resAlg$q
            
            if(length(res$md) == 0) res$md <- 1
            
            
        }
        counter = counter + 1
    }

    
    model$info[[paste0("q", tau)]][[paste0("k", k)]] <- res
    model$beta[[paste0("q", tau)]][[paste0("k", k)]] <- rbind(model$beta[[paste0("q", tau)]][[paste0("k", k)]], BETA)
    if(debug) model$listIH[[paste0("q", tau)]][[paste0("k", k)]] <- list_idx


    return(Ypred)
}