diff --git a/DESCRIPTION b/DESCRIPTION
index 43032fabd65ae2c21104e21ca802ff8b58106f42..f3ada69bb65572829431a240ec8e60a5986da022 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -15,7 +15,8 @@ Imports:
     R6 (>= 2.2.2),
     splines (>= 3.1.1),
     pbs,
-    digest
+    digest,
+    quantreg
 LinkingTo: Rcpp, RcppArmadillo
 Suggests:
     knitr,
diff --git a/R/RcppExports.R b/R/RcppExports.R
index 9bd609395c449dba4099a290a498b10dd87ab940..abbe84c30e475324368021990eddaf71dcb4345b 100644
--- a/R/RcppExports.R
+++ b/R/RcppExports.R
@@ -36,3 +36,15 @@ rls_update_cpp <- function(y, X, theta, P, lambda, k, n, np, istart, kmax) {
     .Call('_onlineforecast_rls_update_cpp', PACKAGE = 'onlineforecast', y, X, theta, P, lambda, k, n, np, istart, kmax)
 }
 
+rq_init_cpp <- function(X, R, beta, n, weights) {
+    .Call('_onlineforecast_rq_init_cpp', PACKAGE = 'onlineforecast', X, R, beta, n, weights)
+}
+
+rq_simplex_cpp <- function(X, Ih, Ihc, IH, K, n, xB, P, tau) {
+    .Call('_onlineforecast_rq_simplex_cpp', PACKAGE = 'onlineforecast', X, Ih, Ihc, IH, K, n, xB, P, tau)
+}
+
+rq_update_cpp <- function(newX, X, IX, Iy, Ih, Ihc, beta, Rny, K, n, xB, P, tau, k, kk, i, n_in_bin, W = 1) {
+    .Call('_onlineforecast_rq_update_cpp', PACKAGE = 'onlineforecast', newX, X, IX, Iy, Ih, Ihc, beta, Rny, K, n, xB, P, tau, k, kk, i, n_in_bin, W)
+}
+
diff --git a/R/data.R b/R/data.R
index ec202c0a3c7b10c3788bfd4b0580521775943986..2ba8d1ffbde2df700922a9fe604abdba77bd6419 100644
--- a/R/data.R
+++ b/R/data.R
@@ -18,3 +18,26 @@
 #' }
 #' @source See \url{https://onlineforecasting.org/examples/datasets.html}.
 "Dbuilding"
+
+
+#' Observations and weather forecasts from a PV plant, weather station and Danish Meteorological Institute (DMI)
+#'
+#' Data of the period from 2008-12-29 to 2011-01-01.
+#'
+#' Hourly average values. The time point is set in the end of the hour.
+#'
+#' Set in the format of a data.list used as input to forecast models in the onlineforecast package.
+#'
+#' @format A data list with 1854 rows and 7 variables:
+#' \describe{
+#'   \item{t}{Time in GMT as POSIXct}
+#'   \item{Ps}{The generated power of a PV plant in MW}
+#'   \item{Iobs}{Observed global radiation at the weather station in W/m^2}
+#'   \item{I}{Weather forecasts of global radiation up to 36 hours ahead from DMI in W/m^2}
+#'   \item{Ta}{Weather forecasts of ambient temperature up to 36 hours ahead from DMI in Celcius}
+#'   \item{Ws}{Weather forecasts of wind speed up to 36 hours ahead from DMI in m/s}
+#'   \item{Wd}{Weather forecasts of wind direction up to 36 hours ahead from DMI in degrees}
+#'   \item{tday}{The time of the day from the given prediction time at the certain prediction horizon in hours}
+#' }
+#' @source See \url{https://onlineforecasting.org/examples/datasets.html}.
+"Dsolar"
\ No newline at end of file
diff --git a/R/pinball_loss.R b/R/pinball_loss.R
new file mode 100755
index 0000000000000000000000000000000000000000..3773a1b7ba7b322ab747757db5a2b40c719e2202
--- /dev/null
+++ b/R/pinball_loss.R
@@ -0,0 +1,3 @@
+pinball_loss <- function(r, tau){
+    return(sum(ifelse(r < 0, (tau - 1)*r, tau*r)))
+}
diff --git a/R/qmodel.R b/R/qmodel.R
new file mode 100644
index 0000000000000000000000000000000000000000..818ab9b70a26256bb8c87db3923eca39d2f01e3a
--- /dev/null
+++ b/R/qmodel.R
@@ -0,0 +1,246 @@
+qmodel <- R6::R6Class("qmodel", public = list(
+      ####### Store Data and information for algorithm #####
+      info = list(),
+      ####### Store Coefficients Data #######
+      beta = list(),
+      ####### Fixed Parameters ######
+      ## Weighted
+      W = NA,
+      n_in_bin = NA,
+      ## Number of predictiors
+      K = NA,
+      ## Quantiles
+      tau = NA,
+      ## Cold start size
+      N1 = NA,
+      debug = NA,
+
+      ####### Quantile Information ########
+      ## Index of the columns of the X matrix of Xny
+      IX = NA,
+      ## Index of the column of y
+      Iy = NA,
+      ## The design matrix
+      X = list(), 
+
+      ####### debug Information ########
+      listIH = list(),
+
+      ####### Forecast Information ########
+      inputs = list(),
+      kseq = NA,
+      output = NA,
+      prm = NA,
+      regprmexpr = NA,
+      regprm = NA,
+      ## Offline parameters
+      prmbounds = as.matrix(data.frame(lower=NA, init =NA, upper=NA)),
+      datatr = NA,
+      qrFIT = NA,
+      maxlagAR = NA,
+      yAR = NA,
+      Ypred = NA,
+
+      initialize = function(N1 = NULL, debug = FALSE){
+          if(is.null(N1)) stop("The number of data points for the cold start needs to be defined")
+          self$N1 <- N1
+          self$debug <-  debug
+      },
+
+#----------------------------------------------------------------
+    # Get the transformation parameters
+    get_prmbounds = function(nm){
+        if(nm == "init"){
+            if(is.null(dim(self$prmbounds))){
+                val <- self$prmbounds[nm]
+            }else{
+                val <- self$prmbounds[ ,nm]
+                if(is.null(nams(val))){
+                    nams(val) <- rownames(self$prmbounds)
+                }
+            }
+        }
+        if(nm == "lower"){
+            if("lower" %in% nams(self$prmbounds)){
+                val <- self$prmbounds[,"lower"]
+                if(is.null(nams(val))){
+                    nams(val) <- rownames(self$prmbounds)
+                }
+            }else{
+                val <- -Inf
+            }
+        }
+        if(nm == "upper"){
+            if("upper" %in% nams(self$prmbounds)){
+                val <- self$prmbounds[,"upper"]
+                if(is.null(nams(val))){
+                    nams(val) <- rownames(self$prmbounds)
+                }
+            }else{
+                val <- Inf
+            }
+        }
+        names(val) <- row.names(self$prmbounds)
+        return(val)
+    },
+    #----------------------------------------------------------------
+
+        #----------------------------------------------------------------
+    # Add the transformation parameters and bounds for optimization
+    add_prmbounds = function(...) {
+        dots <- list(...)
+        for (i in 1:length(dots)) {
+            nm <- names(dots)[i]
+            if (nm %in% rownames(self$prmbounds)) {
+                self$prmbounds[nm, ] <- dots[[i]]
+            } else {
+                if(nrow(self$prmbounds) == 1 & is.na(self$prmbounds[1,2])){
+                    self$prmbounds[1, ] <- dots[[i]]
+                }else{
+                    self$prmbounds <- rbind(self$prmbounds, dots[[i]])
+                }
+                rownames(self$prmbounds)[nrow(self$prmbounds)] <- nm
+            }
+        }
+    },
+    #----------------------------------------------------------------
+
+
+    #----------------------------------------------------------------
+    # Resets the input states
+    reset_state = function(){
+        # Reset the inputs state
+        lapply(self$inputs, function(input){
+            input$state_reset()
+        })
+        # Reset stored data
+        self$datatr <- NA
+        #self$yAR <- NA
+    },
+    #----------------------------------------------------------------
+
+
+
+
+      #----------------------------------------------------------------
+    # Insert the transformation parameters prm in the input expressions and regression expressions, and keep them (simply string manipulation)
+    insert_prm = function(prm){
+        # If just NA or NULL given, then don't do anything
+        if(is.null(prm) | (is.na(prm)[1] & length(prm) == 1)){
+            return(NULL)
+        }
+        # MUST INCLUDE SOME checks here and print useful messages if something is not right
+        if(any(is.na(prm))){ stop(pst("None of the parameters (in prm) must be NA: prm=",prm)) }
+
+        # Keep the prm
+        self$prm <- prm
+        # Find if any opt parameters, first the one with "__" hence for the inputs
+        pinputs <- prm[grep("__",nams(prm))]
+        # If none found for inputs, then the rest must be for regression
+        if (length(pinputs) == 0 & length(prm) > 0) {
+            preg <- prm
+        } else {
+            preg <- prm[-grep("__",nams(prm))]
+        }
+        # ################
+        # For the inputs, insert from prm if any found
+        if (length(pinputs)) {
+            pnms <- unlist(getse(strsplit(nams(pinputs),"__"), 1))
+            pprm <- unlist(getse(strsplit(nams(pinputs),"__"), 2))
+            #
+            for(i in 1:length(self$inputs)){
+                for(ii in 1:length(pnms)){
+                    # Find if the input i have prefix match with the opt. parameter ii
+                    if(pnms[ii]==nams(self$inputs)[i]){
+                        # if the opt. parameter is in the expr, then replace
+                        self$inputs[[i]]$expr <- private$replace_value(name = pprm[ii],
+                                                                       value = pinputs[ii],
+                                                                       expr = self$inputs[[i]]$expr)
+                    }
+                }
+            }
+        }
+        # ################
+        # For the fit parameters, insert from prm if any found
+        if (length(preg) & any(!is.na(self$regprmexpr))) {
+            nams(preg)
+            for(i in 1:length(preg)){
+                # if the opt. parameter is in the expr, then replace
+                self$regprmexpr <- private$replace_value(name = nams(preg)[i],
+                                                         value = preg[i],
+                                                         expr = self$regprmexpr)
+            }
+        }
+        self$regprm <- eval(parse(text = self$regprmexpr))
+    },
+    #----------------------------------------------------------------
+
+      #----------------------------------------------------------------    
+    # Add inputs to the model
+    add_inputs = function(...){
+        dots <- list(...)
+        for (i in 1:length(dots)){
+            self$inputs[[ nams(dots)[i] ]] <-  onlineforecast:::input_class$new(dots[[i]], model=self)
+        }
+    },
+
+    #----------------------------------------------------------------
+    # Add the expression (as character) which generates the regression parameters
+    add_regprm = function(regprmexpr){
+        self$regprmexpr <- regprmexpr
+        self$regprm <- eval(parse(text = self$regprmexpr))
+    },
+    #----------------------------------------------------------------
+
+
+
+
+      #----------------------------------------------------------------
+    # Function for transforming the input data to the regression data
+    transform_data = function(data){
+        # Evaluate for each input the expresssion to generate the model input data
+        L <- lapply(self$inputs, function(input){
+            # Evaluate the expression (input$expr)
+            L <- input$evaluate(data)
+            # Must return a list
+            if(class(L)[1]=="matrix"){ return(list(as.data.frame(L))) }
+            if(class(L)[1]=="data.frame"){ return(list(L)) }
+            if(class(L)[1]!="list"){ stop(pst("The value returned from evaluating: ",input$expr,", was not a matrix, data.frame or a list of them."))}
+            if(class(L[[1]])[1]=="matrix"){ return(lapply(L, function(mat){ return(as.data.frame(mat)) })) }
+            return(L)
+        })
+        # Put together in one data.list
+        L <- structure(do.call(c, L), class="data.list")
+        #
+        return(L)
+    }
+),
+# Private functions
+    private = list(
+    #----------------------------------------------------------------
+
+
+    #----------------------------------------------------------------
+    # Replace the value in "name=value" in expr
+    replace_value = function(name, value, expr){
+        # First make regex
+        pattern <- gsub("\\.", ".*", name)
+        # Try to find it in the input
+        pos <- regexpr(pattern, expr)
+        # Only replace if prm was found
+        if(pos>0){
+            pos <- c(pos+attr(pos,"match.length"))
+            # Find the substr to replace with the prm value
+            (tmp <- substr(expr, pos, nchar(expr)))
+            pos2 <- regexpr(",|)", tmp)
+            # Insert the prm value and return
+            expr <- pst(substr(expr,1,pos-1), "=", value, substr(expr,pos+pos2-1,nchar(expr)))
+            # Print? Not used now
+            #if(printout){ message(names(value),"=",value,", ",sep="")}
+        }
+        return(expr)
+    }
+    #----------------------------------------------------------------
+
+
+))
diff --git a/R/quantile_fit.R b/R/quantile_fit.R
new file mode 100644
index 0000000000000000000000000000000000000000..1bb543222b5f6ea135834f0a34e9fe11c795cad4
--- /dev/null
+++ b/R/quantile_fit.R
@@ -0,0 +1,115 @@
+#' @title quantile_fit
+#'
+#' @description This function plots the quantiles
+#' @param Lfor List of the data needed to plot
+#' @export
+#' @examples
+#' plotprob()
+
+quantile_fit <- function(prm, model, q, data, printout = TRUE, returnanalysis = TRUE, multistep = NA){
+
+    if(printout){
+        message("----------------")
+        if(is.na(prm[1])){
+            message("prm=NA, so current parameters are used.")
+        }else{
+            print_to_message(prm)
+        }
+    }
+    
+    # First insert the prm into the model input expressions
+    model$insert_prm(prm)
+
+    model$n_in_bin <- round(2*(1/(1-model$regprm$lambda)))
+
+    ## Annoying!
+    if(model$n_in_bin < model$N1){
+        model$n_in_bin = model$N1 + 1
+    }
+
+    ## Transform the data
+    
+    # Since rls_fit is run from scratch, the init the stored inputs data (only needed when running iteratively)
+    #model$datatr <- NA
+    #model$yAR <- NA
+    # Reset the model state (e.g. inputs state, stored iterative data, ...)
+    model$reset_state()
+
+    if(q %in% as.numeric(sub(".*q", "", names(model$info)))){
+        model$info <- list()
+        model$beta <- list()
+    }
+
+
+    # Generate the 2nd stage inputs (i.e. the transformed data)
+    model$datatr <- model$transform_data(data)
+    ## Fit the cold start
+    #X <- data.frame(matrix(unlist(model$datatr), ncol = length(model$datatr), byrow = F))
+    #colnames(X) <- names(model$datatr)
+    Ypred <- NULL
+
+    for(k in model$kseq){
+        cat("horizon :", k, "\n")
+        cat("n_in_bin :", model$n_in_bin, "\n")
+        ## Setup the design matrix for k step
+        X <- as.data.frame(subset(model$datatr, kseq=k))
+        ## Lag them to match to the 
+        X <- onlineforecast:::lagdf.matrix(X, k)
+        X$y <- data[[model$output]]
+        
+        X <- X[-1:-k,]
+        fit <- rq(paste("y ~ -1 + ", paste(names(X[, -which(names(X) %in% c("y"))]),  collapse =  " + "), sep = ""), tau = q, data = X[1:(model$N1-k),])
+            
+            
+
+            beta <- as.matrix(fit$coefficients)
+            
+            r <- as.matrix(fit$residuals, ncol = length(q))
+            for(i in 1:dim(r)[2]){
+                ## So these guys are the zeros, or lie on the constraints?
+                IdxNull <- which(abs(r[,i]) < 10^-2)
+                r[IdxNull, i] <- 0
+            }
+            ### Initialization of algorithm
+            model$info[[paste0("q", q)]][[paste0("k", k)]] <- rq_init_cpp(X = as.matrix(X[1:(model$N1-k),]),
+                                                                            R = r,
+                                                                            beta = beta,
+                                                                            n = model$N1-k,
+                                                                            weights = model$regprm$lambda^rev(seq(1,length(r),1)))
+
+            ## The quantile estimation
+            Ypredres <- update_rq(as.matrix(X[(model$N1+1 - k):nrow(X),]), q, k, model = model)
+            ### Get the results, Ypred, loss
+            Ypred <- cbind(Ypred, Ypredres)
+    }
+    
+    ## Compute the loss function, maybe the residuals are compute in different place
+    loss <- sapply(model$kseq, function(i)  pinball_loss(r = model$info[[paste0("q", q)]][[paste0("k", i)]]$R[data$scoreperiod[(i+1):length(data$scoreperiod)]], q))
+    model$Ypred <- Ypred
+
+    if(!returnanalysis){
+        #Yhat <- quantile_prediction(model = model, multistep = model$kseq)
+        cat("Quantile, tau:", q, "\n")
+        cat("the loss: ", loss, "\n")
+        val <- sum(loss)
+        if(printout){
+            print_to_message(c(loss,sum=val))
+        }        
+        return(val)
+    } else{
+        val <- sum(loss)
+        if(printout){
+            print_to_message(c(loss,sum=val))
+        }     
+        ## Get predictions
+        #Yhat <- model$getPred(multistep = multistep)
+        #Yhat <- quantile_prediction(model = model, multistep = multistep)
+
+        #invisible(structure(list(Yhat = Yhat, model = model, data = data, scoreval = loss), class = c("quantile_fit")))
+    }
+}
+
+
+print_to_message <- function(...) {
+    message(paste(utils::capture.output(print(...)), collapse="\n"))
+}
diff --git a/R/quantile_optim.R b/R/quantile_optim.R
new file mode 100644
index 0000000000000000000000000000000000000000..808de1929967184db9cee67623ed2151f9b66404
--- /dev/null
+++ b/R/quantile_optim.R
@@ -0,0 +1,41 @@
+#' @title quantile_optim
+#'
+#' @description This function plots the quantiles
+#' @param Lfor List of the data needed to plot
+#' @export
+#' @examples
+#' plotprob()
+
+
+quantile_optim <- function(model, data){
+    # Take the parameters bounds from the parameter bounds set in the model
+    init <- model$get_prmbounds("init")
+    lower <- model$get_prmbounds("lower")
+    upper <- model$get_prmbounds("upper")
+    # If bounds are NA, then set
+    if(any(is.na(lower))){ lower[is.na(lower)] <- -Inf}
+    if(any(is.na(upper))){ lower[is.na(upper)] <- Inf}
+
+    resList <- list()
+
+    for(i in (1:length(model$tau))){
+        cat("Quantile, tau: ", model$tau[i], "\n")
+    	resList[[i]] <-  optim(par = init,
+                      fn = quantile_fit,
+                      # Parameters to pass to quantile_fit
+                      model = model,
+                      q = model$tau[i],
+                      data = Dnew,
+                      returnanalysis = FALSE,
+                      # Parameters to pass to optim
+                      lower = lower,
+                      upper = upper,
+                   method = "L-BFGS-B")
+
+    }
+
+    #model$optPar <- sapply(1:length(model$tau), function(q) return(resList[[i]]$par))
+
+
+    return(resList)
+}
diff --git a/R/quantile_predict.R b/R/quantile_predict.R
new file mode 100644
index 0000000000000000000000000000000000000000..48c77ecf7bee81151833c6d3d294480c76286e69
--- /dev/null
+++ b/R/quantile_predict.R
@@ -0,0 +1,36 @@
+#' @title quantile_prediction
+#'
+#' @description This function plots the quantiles
+#' @param Lfor List of the data needed to plot
+#' @export
+#' @examples
+#' plotprob()
+
+quantile_predict <- function(model, datatr){
+    Yhat <- lapply(model$tau, function(q){
+        
+        yHat <- sapply(model$kseq, function(k){
+            ## Setup the design matrix for k step
+            X <- as.data.frame(subset(datatr, kseq=k))
+            ## Lag them to match to the 
+            #X <- onlineforecast:::lagdf.matrix(X, k)
+            yhat <- as.numeric(rep(NA, nrow(X)))
+ 
+            for(i in ((k):nrow(X))){
+                if(i <= (model$N1)) {
+                    j <- 1
+                } else{
+                    j <- i - model$N1
+                }
+                #browser()
+                yhat[i] <- t(as.numeric(X[i,])) %*% model$beta[[paste0("q",q)]][[paste0("k",k)]][j,]
+            }
+
+            return(yhat)
+        })
+        nams(yHat) <- pst("k", model$kseq)
+        return(yHat)
+    })
+    names(Yhat) <- pst("q", model$tau)
+    return(Yhat)
+}
diff --git a/R/quantile_update.R b/R/quantile_update.R
new file mode 100644
index 0000000000000000000000000000000000000000..ea98f5342343583cfa90e89915304ec81db9e42e
--- /dev/null
+++ b/R/quantile_update.R
@@ -0,0 +1,26 @@
+quantile_update <- function(model, datatr = NA, y = NA){
+    # Find the number of parameters for the regression
+    np <- length(datatr)
+
+    # Check if data was kept
+    kept_input_data <- !is.na(model$datatr[1])
+    
+    # Check if data was kept
+    kept_input_data <- !is.na(model$datatr[1])
+
+    # Parameters for quantile
+    lambda <- model$regprm$lambda
+    ## n_in_bin could be change!!
+    n_in_bin <-  round(2*(1/(1-lambda)))
+    
+    model$run(X = X[, -which(names(X) %in% c("y"))],
+              Y = X[, "y"],
+              tau = tau,
+              residuals = as.matrix(r),
+              beta = matrix(beta, nrow = length(datatr)),
+              Time = data[["t"]],
+              n_in_bin = round(2*(1/(1-model$regprm$lambda))),
+              Weights = model$regprm$lambda)
+
+    invisible()
+}
diff --git a/R/sortIdx.R b/R/sortIdx.R
new file mode 100644
index 0000000000000000000000000000000000000000..5ecd608a9601069e3e4f870723025d1230d338ce
--- /dev/null
+++ b/R/sortIdx.R
@@ -0,0 +1,13 @@
+#' @title quantile_fit
+#'
+#' @description This function plots the quantiles
+#' @param Lfor List of the data needed to plot
+#' @export
+#' @examples
+#' plotprob()
+
+sortIdx = function(val){
+    newVal <- sort(val)
+    idxVal <- order(val)
+    return(list(newVal = newVal, idxVal = idxVal))
+}
diff --git a/R/update_rq.R b/R/update_rq.R
new file mode 100644
index 0000000000000000000000000000000000000000..d8e09de824ff2bc39d4bfb74d348b8e40286e221
--- /dev/null
+++ b/R/update_rq.R
@@ -0,0 +1,165 @@
+#' @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)
+}
diff --git a/data/Dsolar.rda b/data/Dsolar.rda
new file mode 100644
index 0000000000000000000000000000000000000000..b9e4e413229537ea4876be882a9988d0762c7807
Binary files /dev/null and b/data/Dsolar.rda differ
diff --git a/misc-R/quantile/quantile_Solar.R b/misc-R/quantile/quantile_Solar.R
new file mode 100644
index 0000000000000000000000000000000000000000..772e427a0791ef014ac42f6797888e740a0dcc67
--- /dev/null
+++ b/misc-R/quantile/quantile_Solar.R
@@ -0,0 +1,166 @@
+### Clean up the workspace
+rm(list = ls())
+
+## additional libraries
+library(devtools)
+library(splines)
+library(quantreg)
+
+## Load the onlineforecast package - this is the development version
+load_all(".") # need to be under the package directory
+
+# Load the data
+D <- Dsolar
+names(D)
+
+# Short Data Anlysis
+tstart <- as.POSIXct("2009-01-01", tz = "GMT")
+tend <- as.POSIXct("2011-01-01", tz = "GMT")
+D <- subset(D, in_range(tstart,D$t,tend))
+
+
+plotmat <- matrix(1:3, nrow = 3, byrow=TRUE)
+layout(plotmat, widths=c(1,1,1), heights=c(1,1,1))
+par(mar=c(0, 3.5, 0.5, 1), lwd = 2, cex.lab = 1.5, cex.axis = 1.5)
+
+plot(D$t, D$Ps, type = "l", bty = "l", xlab = "", ylab = "", axes = F)
+axis(2, seq(0,ceiling(max(D$Ps)+0.2), by = 0.5), las = 1)
+
+plot(D$t, lagvec(D$I$k1,1), type = "l", col = "gold", bty = "l", xlab = "", ylab = "", axes = F)
+lines(D$t, D$Iobs, col = rgb(0, 0, 255, max = 255, alpha = 125))
+axis(2, seq(0,ceiling(max(D$I$k1)+0.2), by = 100), las = 1)
+
+plot(D$t, lagvec(D$Ta$k1,1), type = "l", col = rgb(240, 37, 37, max = 255, alpha = 150), bty = "l", xlab = "", ylab = "", axes = F, ylim = c(min(D$Ta$k1) - 4,max(D$Ta$k1)))
+axis(2, seq(floor(min(D$Ta$k1)),ceiling(max(D$Ta$k1)+0.2), by = 10), las = 1)
+axis.POSIXct(side = 1, x = D$t, xaxt = "s",
+        at = seq(D$t[1], D$t[length(D$t)], by = "1 month"),
+        format = "%Y-%m-%d", mgp = c(4,-1,-2))
+
+
+# Make a training set of 3 months, and a then test set with 1 month
+tstart <- "2009-06-01"
+ttest <- "2009-03-31"
+tend <- "2009-07-01"
+
+# Cut only the necessary period
+D$scoreperiod <- in_range("2009-06-10", D$t)
+Dtrain <- subset(D, in_range(tstart,D$t,tend))
+
+
+model <- qmodel$new(N1 = 50)
+model$output = "Ps"
+
+model$add_inputs(I = "bspline(tday, df=5, Boundary.knots=c(5,18)) %**% I")
+
+model$add_regprm("rls_prm(lambda=0.9)")
+model$add_prmbounds(lambda = c(0.99, 0.999, 0.9999))
+
+
+## First 5 columns are the input Matrix
+model$IX <- 0:4
+
+## The output column is the 6th one
+model$Iy <- 5
+
+model$K <- 5
+
+model$datatr <- model$transform_data(Dtrain)
+
+model$tau <- c(0.05,0.2, 0.5, 0.8,0.95)
+model$kseq <- 1:24
+#opt_model <- quantile_optim(model = model, data = Dnew)
+
+plot(Dtrain$I$k1[1:100], type = "l")
+lines(model$datatr$I.bs1$k1, col = "red")
+lines(model$datatr$I.bs2$k1, col = "red")
+lines(model$datatr$I.bs3$k1, col = "red")
+lines(model$datatr$I.bs4$k1, col = "red")
+lines(model$datatr$I.bs5$k1, col = "red")
+
+
+PAR <- c("lambda" = 0.999)
+quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[1])
+PAR <- c("lambda" = 0.999)
+quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[2])
+PAR <- c("lambda" = 0.999)
+quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[3])
+PAR <- c("lambda" = 0.999)
+quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[4])
+PAR <- c("lambda" = 0.999)
+quantile_fit(prm = PAR, model = model, data = Dtrain, q = model$tau[5])
+
+
+Pred_model <- quantile_predict(model = model, datatr = model$datatr)
+
+
+
+
+
+
+str(Pred_model)
+
+k <- 24
+colSeq <- grey(seq(0.9,0.1,len=3))
+#colSeq <- colorJet(nPolygons)
+
+require( RColorBrewer )
+colSeq <- brewer.pal(11 , "Spectral" )
+
+for(i in k:(length(Dtrain$t)-1)){
+
+    par(mfrow = c(2,1))
+    par(mar = c(2,4,2,2)) # bottom left top right
+    
+    with(Pred_model, {
+
+        plot(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], bty = "l", lwd = 2, col = "black", pch=19, cex=0.5,  axes = FALSE, xaxt = "n",
+             type="n", ylim = range(Dtrain$Ps[(i-10):(i+k)], q0.5[i,],  q0.95[i,], q0.05[i,], na.rm = T), main = "TimeAdaptive with Weights", xlab = "Time", ylab = "Ps")
+        
+        axis(2)
+        
+        
+        
+        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
+                c(Pred_model$q0.05[i,],rev(Pred_model$q0.2[i,])), col= colSeq[3], border = NA)
+        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
+                c(Pred_model$q0.8[i,],rev(Pred_model$q0.95[i,])), col= colSeq[3], border = NA)
+
+        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
+               c(Pred_model$q0.2[i,],rev(Pred_model$q0.5[i,])), col= colSeq[4], border = NA)
+        polygon(c(Dtrain$t[(i+1):(i+k)], rev(Dtrain$t[(i+1):(i+k)])),
+               c(Pred_model$q0.5[i,],rev(Pred_model$q0.8[i,])), col= colSeq[4], border = NA)
+        
+        
+        lines(Dtrain$t[(i-10):(i+k)], Dtrain$Ps[(i-10):(i+k)], col = "black", type = "b", lwd = 2)
+        lines(Dtrain$t[(i+1):(i+k)], q0.5[i,], type = "b", col = "grey", lwd = 2)
+
+        axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)], 
+                     at = seq(from = Dtrain$t[(i-10):(i+k)][1],
+                              to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])],
+                              by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a",
+                     las = 1, cex.axis = 1, srt = 45)
+        
+        #lines(Pred_model$q0.05[i,], type = "b", col = "blue")
+        #lines(Pred_model$q0.95[i,], type = "b", col = "blue")
+
+        plot(Dtrain$t[(i-10):(i+k)], c(rep(NA,11), as.numeric(Dtrain$I[i, 1:k])), type = "b", col = "steelblue", axes = FALSE, xlab = "Time", ylab = "Temp", lwd = 2, ylim = range(Dtrain$I[i, 1:k],  Dtrain$Iobs[(i+1):(i+k)]))
+        #lines(Dtrain$t[i], Dtrain$Iobs[i+1,1], type = "b", col = "red", axes = FALSE)
+        lines(Dtrain$t[(i+1):(i+k)], Dtrain$Iobs[(i+1):(i+k)], type = "b", col = "red", lwd = 2)
+        axis(2)
+        abline(h=0,v=Dtrain$t[i],lty=2, col = "lightgrey", lwd = 2)
+        
+        legend("topleft", legend=c("I Obs", "I Pred"),
+               col=c("red", "steelblue"), lty=1, cex=1, lwd = 2)
+    
+        
+i 
+        axis.POSIXct(side = 1, x = Dtrain$t[(i-10):(i+k)], 
+                     at = seq(from = Dtrain$t[(i-10):(i+k)][1],
+                              to = Dtrain$t[(i-10):(i+k)][length(Dtrain$t[(i-10):(i+k)])],
+                              by = "1 hour"), format = "%Y/%m/%d %H:%M \n %a",
+                     las = 1, cex.axis = 1, srt = 45)
+    })
+
+
+}
+
diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp
index 20e61f5a6db5a3c3525e52ed4f7bbe15ba0aba45..e5f7c71f5b5e4e222fd00ed7effdea7cddd9955d 100644
--- a/src/RcppExports.cpp
+++ b/src/RcppExports.cpp
@@ -43,10 +43,75 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// rq_init_cpp
+Rcpp::List rq_init_cpp(arma::mat X, arma::vec R, arma::vec beta, unsigned int n, arma::vec weights);
+RcppExport SEXP _onlineforecast_rq_init_cpp(SEXP XSEXP, SEXP RSEXP, SEXP betaSEXP, SEXP nSEXP, SEXP weightsSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type R(RSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type beta(betaSEXP);
+    Rcpp::traits::input_parameter< unsigned int >::type n(nSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type weights(weightsSEXP);
+    rcpp_result_gen = Rcpp::wrap(rq_init_cpp(X, R, beta, n, weights));
+    return rcpp_result_gen;
+END_RCPP
+}
+// rq_simplex_cpp
+Rcpp::List rq_simplex_cpp(arma::mat X, arma::uvec Ih, arma::uvec Ihc, arma::mat IH, unsigned int K, unsigned int n, arma::vec xB, arma::vec P, double tau);
+RcppExport SEXP _onlineforecast_rq_simplex_cpp(SEXP XSEXP, SEXP IhSEXP, SEXP IhcSEXP, SEXP IHSEXP, SEXP KSEXP, SEXP nSEXP, SEXP xBSEXP, SEXP PSEXP, SEXP tauSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP);
+    Rcpp::traits::input_parameter< arma::uvec >::type Ih(IhSEXP);
+    Rcpp::traits::input_parameter< arma::uvec >::type Ihc(IhcSEXP);
+    Rcpp::traits::input_parameter< arma::mat >::type IH(IHSEXP);
+    Rcpp::traits::input_parameter< unsigned int >::type K(KSEXP);
+    Rcpp::traits::input_parameter< unsigned int >::type n(nSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type xB(xBSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type P(PSEXP);
+    Rcpp::traits::input_parameter< double >::type tau(tauSEXP);
+    rcpp_result_gen = Rcpp::wrap(rq_simplex_cpp(X, Ih, Ihc, IH, K, n, xB, P, tau));
+    return rcpp_result_gen;
+END_RCPP
+}
+// rq_update_cpp
+Rcpp::List rq_update_cpp(arma::mat newX, arma::mat X, arma::uvec IX, arma::uvec Iy, arma::uvec Ih, arma::uvec Ihc, arma::vec beta, arma::vec Rny, unsigned int K, unsigned int n, arma::vec xB, arma::vec P, double tau, unsigned int k, arma::uvec kk, unsigned int i, unsigned int n_in_bin, double W);
+RcppExport SEXP _onlineforecast_rq_update_cpp(SEXP newXSEXP, SEXP XSEXP, SEXP IXSEXP, SEXP IySEXP, SEXP IhSEXP, SEXP IhcSEXP, SEXP betaSEXP, SEXP RnySEXP, SEXP KSEXP, SEXP nSEXP, SEXP xBSEXP, SEXP PSEXP, SEXP tauSEXP, SEXP kSEXP, SEXP kkSEXP, SEXP iSEXP, SEXP n_in_binSEXP, SEXP WSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat >::type newX(newXSEXP);
+    Rcpp::traits::input_parameter< arma::mat >::type X(XSEXP);
+    Rcpp::traits::input_parameter< arma::uvec >::type IX(IXSEXP);
+    Rcpp::traits::input_parameter< arma::uvec >::type Iy(IySEXP);
+    Rcpp::traits::input_parameter< arma::uvec >::type Ih(IhSEXP);
+    Rcpp::traits::input_parameter< arma::uvec >::type Ihc(IhcSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type beta(betaSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type Rny(RnySEXP);
+    Rcpp::traits::input_parameter< unsigned int >::type K(KSEXP);
+    Rcpp::traits::input_parameter< unsigned int >::type n(nSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type xB(xBSEXP);
+    Rcpp::traits::input_parameter< arma::vec >::type P(PSEXP);
+    Rcpp::traits::input_parameter< double >::type tau(tauSEXP);
+    Rcpp::traits::input_parameter< unsigned int >::type k(kSEXP);
+    Rcpp::traits::input_parameter< arma::uvec >::type kk(kkSEXP);
+    Rcpp::traits::input_parameter< unsigned int >::type i(iSEXP);
+    Rcpp::traits::input_parameter< unsigned int >::type n_in_bin(n_in_binSEXP);
+    Rcpp::traits::input_parameter< double >::type W(WSEXP);
+    rcpp_result_gen = Rcpp::wrap(rq_update_cpp(newX, X, IX, Iy, Ih, Ihc, beta, Rny, K, n, xB, P, tau, k, kk, i, n_in_bin, W));
+    return rcpp_result_gen;
+END_RCPP
+}
 
 static const R_CallMethodDef CallEntries[] = {
     {"_onlineforecast_lp_vector_cpp", (DL_FUNC) &_onlineforecast_lp_vector_cpp, 2},
     {"_onlineforecast_rls_update_cpp", (DL_FUNC) &_onlineforecast_rls_update_cpp, 10},
+    {"_onlineforecast_rq_init_cpp", (DL_FUNC) &_onlineforecast_rq_init_cpp, 5},
+    {"_onlineforecast_rq_simplex_cpp", (DL_FUNC) &_onlineforecast_rq_simplex_cpp, 9},
+    {"_onlineforecast_rq_update_cpp", (DL_FUNC) &_onlineforecast_rq_update_cpp, 18},
     {NULL, NULL, 0}
 };
 
diff --git a/src/rq_init_cpp.cpp b/src/rq_init_cpp.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..4a780b2f0bb2a4663db4393f1dec95ffb0675a49
--- /dev/null
+++ b/src/rq_init_cpp.cpp
@@ -0,0 +1,56 @@
+#include "RcppArmadillo.h"
+#include <Rcpp.h>
+// [[Rcpp::depends(RcppArmadillo)]]
+using namespace Rcpp;
+using namespace std; 
+
+// [[Rcpp::export]]
+Rcpp::List rq_init_cpp(arma::mat X,
+					   arma::vec R,
+					   arma::vec beta, 
+					   unsigned int n,
+					   arma::vec weights){
+
+	//arma::vec idx(n);
+	arma::mat L, U, Pi;
+	unsigned int Lr0 = count(R.begin(), R.end(), 0);
+	//unsigned int idx;
+
+	if(Lr0 > beta.n_elem){
+
+		Rcout << "Worked" << "\n";
+		arma::uvec Irs = sort_index(abs(R));
+		Irs = Irs.subvec(0,Lr0-1);
+		arma::mat Xh = X.rows(Irs);
+
+		lu(L, U, Pi, Xh);
+		arma::vec In = arma::linspace(0,Lr0,Lr0+1);
+		arma::vec rI = arma::linspace(0, (Lr0)-beta.n_elem-1, (Lr0)-beta.n_elem);
+
+
+		for(unsigned int i = (beta.n_elem); i < Lr0; i++){
+			//idx = In(find(Pi.row(i) == 1));
+			R(Irs(find(Pi.row(i) == 1))).fill(1.e-15);
+
+		}
+		//R(Irs(rI)) = 1.e-15;
+	}
+
+	// Find all index where r == 0 - Vertices
+	arma::uvec Ih = find(R == 0);
+	arma::uvec Ihc = find(abs(R) > 0);
+	arma::vec P = sign(R.elem(Ihc));
+	R.elem(find(abs(R) < 1.e-15)).fill(0);
+
+	arma::vec W = weights % R;
+
+
+	arma::vec xB = join_cols(beta, abs(W(Ihc)));
+
+	return List::create(Named("xB") = xB,
+						Named("Ih") = Ih,
+						Named("Ihc") = Ihc,
+						Named("P") = P,
+						Named("R") = R,
+						Named("X") = X);
+}
\ No newline at end of file
diff --git a/src/rq_simplex_cpp.cpp b/src/rq_simplex_cpp.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..5890c20ea03044547e26341fc5c957b5c8a8500c
--- /dev/null
+++ b/src/rq_simplex_cpp.cpp
@@ -0,0 +1,145 @@
+#include "RcppArmadillo.h"
+#include <Rcpp.h>
+// [[Rcpp::depends(RcppArmadillo)]]
+using namespace Rcpp;
+using namespace std; 
+
+// [[Rcpp::export]]
+Rcpp::List rq_simplex_cpp(arma::mat X,
+					   	arma::uvec Ih,
+					   	arma::uvec Ihc,
+					   	arma::mat IH,
+					   	unsigned int K, 
+					   	unsigned int n,
+					   	arma::vec xB,
+					   	arma::vec P,
+					   	double tau){
+
+
+	arma::mat invXh = inv(X.rows(Ih));
+	arma::vec cB = arma::conv_to<arma::vec>::from(P < 0) + P * tau;
+
+	arma::vec cC = arma::join_cols(arma::ones<arma::vec>(K) * tau, arma::ones<arma::vec>(K) * (1-tau));
+
+
+
+	//Rcout << Xny .rows(Ihc) << "\n";
+	arma::mat IB2 = -(P * arma::ones<arma::vec>(K).t() % X.rows(Ihc))*invXh;	
+
+	arma::vec g = IB2.t() * cB;
+	arma::vec d = cC - arma::join_cols(g,-g);
+	d.elem(find(abs(d) < 1.e-15)).fill(0);
+	arma::uvec s = sort_index(d);
+	arma::vec md = sort(d);
+
+
+	s = s.elem(find(md < 0));
+	md = md.elem(find(md < 0));
+
+
+
+	arma::vec c = arma::ones<arma::vec>(s.size());
+
+	for(unsigned int i = 0; i < s.size(); i++){
+		if(s(i)  >= K){
+			s(i) = s(i) - K;
+			c(i) = -1;
+		}
+	}
+
+
+	//c.elem(s>K).fill(-1);
+	
+
+
+
+	arma::mat C(c.size(),c.size(),arma::fill::eye);
+	C.diag() = c;
+
+
+	arma::mat h = arma::join_cols(invXh.cols(s), IB2.cols(s)) * C;
+
+
+	arma::vec xm = xB.elem(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K));
+	xm.elem(find(xm < 0)).fill(0); 
+	arma::mat hm = h.rows(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K));
+	
+
+
+	arma::vec alpha(s.size());
+	arma::uvec q(s.size());
+	arma::vec cq(s.size());
+	arma::uvec idx2 = arma::linspace<arma::uvec>(0,s.size()-1,s.size());
+	
+
+
+	for(unsigned int k = 0; k < s.size(); k++){
+
+		arma::vec sigma = xm;
+		arma::uvec idx = find(hm.col(k) > 1.e-12);
+
+		sigma.elem(idx) = xm.elem(idx) / hm(idx, find(idx2 == k));
+
+		sigma(find(hm.col(k) <= 1.e-12)).fill(1.e30);
+
+
+
+		alpha.elem(find(idx2 == k)).fill(sigma.min());
+		q.elem(find(idx2 == k)).fill(sigma.index_min());
+
+		cq.elem(find(idx2 == k)) = c.elem(find(idx2 == k));
+	}
+
+
+	arma::vec gain = md % alpha;
+	arma::vec Mgain = sort(gain);
+	arma::uvec IMgain = sort_index(gain);
+	arma::uvec IhMid;
+
+	double CON = 1.e30;
+	unsigned int j = 0;
+
+
+
+	if(gain.size() < 1){
+		gain.resize(1);
+		gain(0) = 1;
+	} else {
+		while((CON > 1.e6) & (j < s.size())){
+			IhMid = Ih;
+			IhMid(s(IMgain(j))) = Ihc(q(IMgain(j))); 
+			IhMid = sort(IhMid);
+
+			if(sum(abs(IH.t() - IhMid.t())) == 0){
+				CON = 1.e30;
+			} else{
+				CON = cond(X.rows(IhMid));
+			}
+			s = s(IMgain(j));
+			q = q(IMgain(j));
+			cq = cq(IMgain(j));
+			alpha = alpha(IMgain(j));
+			IH = arma::join_rows(IH, arma::conv_to<arma::vec>::from(IhMid));
+			h = h.col(IMgain(j));
+			gain = gain(IMgain(j));
+			md = md(IMgain(j));
+
+
+			j += 1;
+
+		}
+
+	}
+
+
+	return List::create(Named("CON") = CON,
+						Named("s") = s,
+						Named("g") = g,
+						Named("q") = q,
+						Named("gain") = gain,
+						Named("md") = md,
+						Named("alpha") = alpha,
+						Named("h") = h,
+						Named("IH") = IH,
+						Named("cq") = cq);
+}
\ No newline at end of file
diff --git a/src/rq_update_cpp.cpp b/src/rq_update_cpp.cpp
new file mode 100755
index 0000000000000000000000000000000000000000..17ca40e891255604715d911c0a6e88fa747221eb
--- /dev/null
+++ b/src/rq_update_cpp.cpp
@@ -0,0 +1,179 @@
+#include "RcppArmadillo.h"
+#include <Rcpp.h>
+// [[Rcpp::depends(RcppArmadillo)]]
+using namespace Rcpp;
+using namespace std; 
+
+// [[Rcpp::export]]
+Rcpp::List rq_update_cpp(arma::mat newX,
+					   	arma::mat X,
+					   	arma::uvec IX,
+					   	arma::uvec Iy,
+					   	arma::uvec Ih,
+					   	arma::uvec Ihc,
+					   	arma::vec beta,
+					   	arma::vec Rny,
+					   	unsigned int K, 
+					   	unsigned int n,
+					   	arma::vec xB,
+					   	arma::vec P,
+					   	double tau,
+					   	unsigned int k,
+					   	arma::uvec kk,
+					   	unsigned int i,
+					   	unsigned int n_in_bin,
+					   	double W = 1){
+
+	// Init 
+	arma::vec xBm;
+
+	// If weights are given - then weight the past Xny matrix
+	// Notice, tho the newest observaiton is not weight - maybe need a fix?
+	if(W != 1){
+		X = W * X;
+	}
+
+	X = join_cols(X, newX);
+
+	// Rolling window of leaving the design matrix
+	// i.e. the last observation leaves the matrix 
+	// See Møller 2008 if something else is needed
+	unsigned int Leav = 0;
+	arma::uvec sortIdx = sort_index(Ih);
+
+	//Rcout << "min" << sortIdx << "\n";
+	//Rcout << "minV"  << sortIdx(0) << "\n";
+
+	if(Ih(sortIdx(0)) == Leav){
+		//Rcout << "Here but still hmm" << "\n";
+		//Rcout << "Ih" << Ih << "\n";
+		// Compute the inverse of X(h)
+		arma::mat invXh = inv(X(Ih, IX));
+		//Rcout << invXh.col(sortIdx(0)) << "\n";
+
+		// Computing the c
+		arma::vec cB = arma::conv_to<arma::vec>::from(P < 0) + P * tau;
+		// Computing inverse (B^T) * c_B - P * X(bar(h)) * X(b)^-1 * inverse(X(bar(h)))
+		double g = arma::as_scalar(-(P % (X(Ihc,IX) * invXh.col(sortIdx(0)))).t() * cB);
+		int signG = (g > 0) - (g < 0);
+		//  h is the direction or the d in the paper
+		arma::vec h = join_cols(invXh.col(sortIdx(0)), -(P % (X(Ihc,IX) * invXh.col(sortIdx(0))))) * signG;
+		// Compute the magitude of the step, alpha
+		arma::vec sigma = arma::zeros(n - K);
+		arma::vec hm = h.elem(arma::linspace<arma::uvec>(K,h.size() - 1,h.size() - K));
+
+ 		// The absolute of the residuals which are not equal to zero
+		arma::vec xBm = xB.elem(arma::linspace<arma::uvec>(K,h.size()-1,h.size()-K));
+		xBm.elem(find(xBm < 0)).fill(0);
+		//Rcout << "Here but still hmm" << "\n";
+
+		//Rcout << "hm" << hm << "\n";
+		//Rcout << "Sigma" << sigma << "\n";
+
+		//Rcout << "hm" << hm.size() << "\n";
+		//Rcout << "Sigma" << sigma.size() << "\n";
+
+		sigma.elem(find(hm > 1.e-10)) = xBm.elem(find(hm>1.e-10)) / hm.elem(find(hm>1.e-10));
+		sigma.elem(find(hm<= 1.e-10)).fill(1.e30);
+
+		//Rcout << "min idx " << index_min(sigma) << "\n";
+
+		unsigned int q = index_min(sigma); //sort_index(sigma).elem(0);
+	
+		double alpha =  arma::as_scalar(sigma.row(q));
+
+		arma::vec z = xB - alpha * h;
+		arma::uvec Ihm = Ih.row(sortIdx(0));
+		Ih.row(sortIdx(0)) = Ihc.row(q);
+		Ihc(q) = Ihm(0);
+
+		xB = z;
+		xB(q+K) = alpha;
+
+		P(q) = signG + (g == 0);
+		Ih = sort(Ih);
+		arma::uvec IndexIhc = sort_index(Ihc);
+		Ihc = sort(Ihc);
+		P = P(IndexIhc);
+
+		xBm = xB(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K));
+		xBm = xBm(IndexIhc);
+		xB(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K)) = xBm;
+		beta = xB(arma::linspace<arma::uvec>(0,K-1,K));
+	}
+
+	//Rcout << "IX " << IX <<  "\n";
+	//Rcout << "X " << X(kk, IX) <<  "\n";
+	//Rcout << "Pred " << arma::as_scalar(X(kk, IX) * beta) <<  "\n";
+
+	double Ypred = arma::as_scalar(newX(kk, IX) * beta);
+	double rny = arma::as_scalar(newX(kk, Iy)) - Ypred;
+
+	//Rcout << "rny" << rny << "\n";
+
+	int szR = Rny.size();
+	Rny.resize(szR+1);
+	Rny(szR) = rny;
+
+
+	int szP = P.size();
+	P.resize(szP+1);
+	int szB = xB.size();
+	xB.resize(szB+1);
+
+
+	if(rny < 0){
+		P(szP) = -1;
+		xB(szB) = -rny;
+	} else {
+		P(szP) = 1;
+		xB(szB) = rny;
+	}
+
+	int sz = Ihc.size();
+	Ihc.resize(sz+1);
+	Ihc(sz) = n;
+
+
+	//Rcout << "Ihc" << Ihc << "\n";
+
+	if(X.n_rows == n_in_bin){
+		i += 1;
+		arma::uvec stay = arma::linspace<arma::uvec>(1,P.size()-1,P.size()-1);
+
+		P = P(stay);
+		Ihc = Ihc(stay);
+		//arma::vec Pnew = P(stay);
+		//arma::uvec Ihcnew = Ihc(stay);
+
+		X = X.rows(sort(join_cols(Ihc,Ih)));
+
+
+		xBm = xB.elem(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K));
+		xBm = xBm(stay);
+		xB = xB.elem(arma::linspace<arma::uvec>(0,xB.size()-2,xB.size()-1));
+		xB.elem(arma::linspace<arma::uvec>(K,xB.size()-1,xB.size()-K)) = xBm;
+		
+		//Ihc = Ihc(stay);
+
+		Ihc = Ihc - 1;
+		Ih = Ih - 1;
+
+
+	} else{
+		n += 1;
+	}
+
+
+
+	return List::create(Named("Ih") = Ih,
+		Named("Ihc") = Ihc,
+		Named("xB") = xB,
+		Named("X") = X,
+		Named("R") = Rny,
+		Named("P") = P,
+		Named("n") = n,
+		Named("i") = i,
+		Named("Ypred") = Ypred);
+
+}
\ No newline at end of file