"rqss.fit" <-
function (x, y, tau = 0.5, method = "sfn", rhs = NULL, control, ...)
{
    if(is.null(rhs)) rhs <- (1 - tau) * t(x) %*% rep(1,length(y))
    else tau <- 0.5
    fit <- switch(method,
	sfn = rq.fit.sfn(x, y, tau = tau,  rhs = rhs, control = control, ...),
        sfnc = rq.fit.sfnc(x, y, tau = tau,  rhs = rhs, control = control, ...), {
            what <- paste("rq.fit.", method, sep = "")
            if (exists(what, mode = "function"))
                (get(what, mode = "function"))(x, y, ...)
            else stop(paste("unimplemented method:", method))
        })
    fit$contrasts <- attr(x, "contrasts")
    fit$resid <- c(y - x %*% fit$coef)

    fit
}
"untangle.specials" <-
function (tt, special, order = 1)
{
    spc <- attr(tt, "specials")[[special]]
    if (length(spc) == 0)
        return(list(vars = character(0), terms = numeric(0)))
    facs <- attr(tt, "factor")
    fname <- dimnames(facs)
    ff <- apply(facs[spc, , drop = FALSE], 2, sum)
    list(vars = (fname[[1]])[spc], terms = seq(ff)[ff & match(attr(tt,
        "order"), order, nomatch = 0)])
}

"qss" <-
function (x, constraint = "N", lambda = 1, ndum = 0, dummies = NULL, 
    w = rep(1, length(x))) 
{
    if (is.matrix(x)) {
        if (ncol(x) == 2) 
            qss <- qss2(x, constraint = constraint, dummies = dummies, 
                lambda = lambda, ndum = ndum, w = w)
        else if (ncol(x) == 1) 
            x <- as.vector(x)
        else stop("qss objects must have dimension 1 or 2")
    }
    if (is.vector(x)) 
        qss <- qss1(x, constraint = constraint, lambda = lambda, 
            dummies = dummies, ndum = ndum, w = w)
    qss
}

"qss2" <-
function(x, y, constraint = "N", lambda = 1, ndum= 0, dummies = NULL, w=rep(1,length(x))){
#
# Sparse Additive Quantile Smoothing Spline Models - Bivariate (Triogram) Module
#
# This function returns a structure intended to make model.matrix for a bivariate
# nonparametric component of a model formula specified by a call to rqss().  A sparse form
# of the Frisch Newton algorithm is eventually called to compute the estimator.
# An optional convexity/concavity constraint can be specified.  If
# the formula consists of a single qss component then the estimator solves the
# following variational problem:
#
#       min sum rho_tau (y_i - g(x_i)) + lambda V(grad(g))
#
# where V(f) denotes the total variation of the function f.  The solution is a piecewise
# linear function on the Delaunay triangulation formed by the observed (x_i,y_i) points.
# Additive models can consist
# of several components of this form plus partial linear and univariate qss components.
# To resolve the identifiability problem we delete the first column of the qss design
# components.  On return F contains the fidelity portion of the design, A the penalty
# contribution of the design, R the constraint portion, and r the rhs of the constraints.
#
# Constraints are specified by the constraint argument:
#
#	N	none
#	V	convex
#	C	concave
#
# Author:  Roger Koenker   April 2, 2003
#
# For a prototype see triogram in ~roger/projects/tv/cobar/.RData on ysidro.
#
# Warning:   Under development...todo:
#
#       o  weights
#       o  dummy x's
#       o  tau's
#       o  lambda's
#       o  ...
#
   stopifnot(requireNamespace("tripack"))
#
    y <- x[,2]
    x <- x[,1]
    n <- length(x)
    if (n != length(y))
        stop("xy lengths do not match")
    f <- triogram.fidelity(x, y, ndum = ndum, dummies = dummies)
    F <- f$F
    A <- triogram.penalty(f$x, f$y)
    switch(constraint, V = {
        R <- A
        r <- rep(0, nrow(R))
    }, C = {
        R <- -A
        r <- rep(0, nrow(R))
    }, N = {
        R = NULL
        r = NULL
    })
    list(x = list(x = f$x, y = f$y), F = F[, -1], dummies = f$dummies,
        lambda = lambda, A = A[, -1], R = R[, -1], r = r)
}


"qss1" <-
function (x, constraint = "N", lambda = 1, dummies = dummies,
        ndum = 0, w = rep(1, length(x))){
#
# Sparse Additive Quantile Smoothing Spline Models - Univariate Module
#
# This function returns a structure intended to make model.matrix for a univariate
# nonparametric component of a model formula specified by a call to rq().  A sparse form
# of the Frisch Newton algorithm is eventually called to compute the estimator.
# Optional monotonicity and/or convexity/concavity constraints can be specified.  If
# the formula consists of a single qss component then the estimator solves the
# following variational problem:
#
#       min sum rho_tau (y_i - g(x_i)) + lambda V(g')
#
# where V(f) denotes the total variation of the function f.  The solution is a piecewise
# linear function with "knots" at the observed x_i points.  Additive models can consist
# of several components of this form plus partial linear and triogram components.
# To resolve the identifiability problem we delete the first column of the qss design
# components.  On return F contains the fidelity portion of the design, A the penalty
# contribution of the design, R the constraint portion, and r the rhs of the constraints.
#
# Constraints are specified by the constraint argument:
#
#	N	none
#	I	monotone increasing
#	D	monotone decreasing
#	V	convex
#	C	concave
#	CI	concave and monotone increasing
#	...	etc
#
# Author:  Roger Koenker   February 27, 2003
#
# Warning:   Under development...todo:
#
#       o  weights
#       o  dummy x's
#       o  tau's
#       o  lambda's
#       o  ...
#
#
    xun <- unique(x[order(x)])
    h <- diff(xun)
    nh <- length(h)
    nx <- length(x)
    p <- nh + 1
    B <- new("matrix.csr", ra = c(rbind(-1/h, 1/h)), ja = as.integer(c(rbind(1:nh,
        2:(nh + 1)))), ia = as.integer(2 * (1:(nh + 1)) - 1),
        dimension = as.integer(c(nh, nh + 1)))
    makeD <- function(p) {
        new("matrix.csr", ra = c(rbind(rep(-1, (p - 1)), rep(1,
            (p - 1)))), ja = as.integer(c(rbind(1:(p - 1), 2:p))),
            ia = as.integer(2 * (1:p) - 1), dimension = as.integer(c(p -
                1, p)))
    }
    D <- makeD(nh)
    A <- D %*% B
    if (length(xun) == length(x)){
        F <- new("matrix.csr", ra = rep(1, nx), ja = as.integer(rank(x)),
            ia = 1:(nx + 1), dimension = as.integer(c(nx, nx)))
	}
    else {
        F <- new("matrix.csr", ra = rep(1, nx), ja = as.integer(factor(x)),
            ia = 1:(nx + 1), dimension = as.integer(c(nx, length(xun))))
	}

   switch(constraint,
        V = {   R <- A;
                r <- rep(0,nrow(R))
                },
        C = {   R <- -A;
                r <- rep(0,nrow(R))
                },
        I = {   R <- makeD(p)
                r <- rep(0,p-1)
                },
        D = {   R <- -makeD(p)
                r <- rep(0,p-1)
                },
        VI = {  R <- makeD(p)
		R <- rbind(R,A)
                r <- rep(0,nrow(R))
		},
        VD = {  R <- -makeD(p)
		R <- rbind(R,A)
                r <- rep(0,nrow(R))
		},
        CI = {  R <- makeD(p)
		R <- rbind(R,-A)
                r <- rep(0,nrow(R))
		},
        CD = {  R <- -makeD(p)
		R <- rbind(R,-A)
                r <- rep(0,nrow(R))
		},
	N = { R=NULL; r=NULL}
	)
   list(x = list(x=xun), F=F[,-1], lambda = lambda, A=A[,-1], R=R[,-1], r=r)
}
"plot.qss1" <-
function(x, rug = TRUE, jit = TRUE, add = FALSE, ...)
{
if(!add) plot(x[,1],x[,2],type = "n", ...)
lines(x[,1],x[,2], ...)
if(rug) {
    if(jit)
       rug(jitter(x[,1]))
    else rug(x[,1])
    }
}
"plot.qss2" <-
function (x, render = "contour", ncol = 100, zcol = NULL, ...)
{
    stopifnot(requireNamespace("tripack"))
    y <- x[, 2]
    z <- x[, 3]
    x <- x[, 1]
    tri <- tripack::tri.mesh(x, y)
    if (render == "rgl") {
        if(!requireNamespace("rgl",quietly=TRUE))
		stop("The package rgl is required")
        collut <- terrain.colors(ncol)
        if (!length(zcol))
            zcol <- z
        if (max(z) > max(zcol) || min(z) < min(zcol))
            warning("fitted z values out of range of zcol vector")
        zlim <- range(zcol)
        colz <- ncol * (z - zlim[1])/(zlim[2] - zlim[1]) + 1
        colz <- collut[colz]
        s <- c(t(tripack::triangles(tri)[, 1:3]))
        rgl::rgl.triangles(x[s], y[s], z[s], col = colz[s])
    }
    else {
        stopifnot(requireNamespace("akima"))
        if(render == "contour"){
                plot(x, y, type = "n", ...)
                contour(akima::interp(x, y, z), add = TRUE, frame.plot = TRUE, ...)
                tripack::convex.hull(tri, plot.it = TRUE, add = TRUE)
                }
        else if(render == "persp")
                persp(akima::interp(x, y, z, ), theta = -40, phi = 20, xlab = "x",
                        ylab = "y", zlab = "z", ...)
        else stop(paste("Unable to render: ",render))
    }
}
plot.rqss <-
function (x, rug = TRUE, jit = TRUE, bands = NULL, coverage = 0.95, add = FALSE, 
	  shade = TRUE, select = NULL, pages = 0, titles = NULL, bcol = NULL, ...) 
{
    SetLayout <- function(m, p) {
            # Shamelessly cribbed from mgcv
            # m is the number of plots
            # p is the number of pages
        if (p > m) 
            p <- m
        if (p < 0) 
            p <- 0
        if (p != 0) {
            q <- m%/%p
            if ((m%%p) != 0) {
                q <- q + 1
                while (q * (p - 1) >= m) p <- p - 1
            }
            c <- trunc(sqrt(q))
            if (c < 1) 
                c <- 1
            r <- q%/%c
            if (r < 1) 
                r <- 1
            while (r * c < q) r <- r + 1
            while (r * c - q > c && r > 1) r <- r - 1
            while (r * c - q > r && c > 1) c <- c - 1
            oldpar <- par(mfrow = c(r, c))
        }
        else oldpar <- par()
        return(oldpar)
    }
    m <- length(x$qss)
    if (m == 0) 
        stop("No qss object to plot")
    if(length(select)) {
        if(all(select %in% 1:m))
            oldpar <- SetLayout(length(select), pages)
        else
            stop(paste("select must be in 1:",m,sep="")) 
        }
    else
        oldpar <- SetLayout(m, pages)
    if ((pages == 0 && prod(par("mfrow")) < m && dev.interactive()) || 
        pages > 1 && dev.interactive()) 
        ask <- TRUE
    else ask <- FALSE
    if (ask) {
        oask <- devAskNewPage(TRUE)
        on.exit(devAskNewPage(oask))
    }
    qssnames <- names(x$qss)
    if(length(titles)){
        if(length(titles) != length(qssnames))
            stop("Length of titles doesn't match length of qssnames")
	}
    else{
        titles <- paste("Effect of ", qssnames)
	}
    if (length(bands)) {
        rdf <- x$n - x$edf
        if (any(unlist(lapply(x$qss, function(x) ncol(x$xyz) == 3)))) 
            warning("Can't plot confidence bands in 3D (yet)")
        band <- as.list(rep(NA, m))
        V <- summary(x, cov = TRUE, ...)$Vqss
        "summary.qss1" <- function(object, V, ngrid = 400, ...) {
            x <- object$xyz[, 1]
            eps <- 0.01
            newd <- data.frame(x = seq(min(x) + eps, max(x) - 
                eps, length = ngrid))
            G <- predict(object, newd, ...)
            ones <- as.matrix.csr(matrix(1, nrow(G$D), 1))
            D <- cbind(ones, G$D)
            S <- as.matrix(D %*% V %*% t(D))
            se <- sqrt(diag(S))
            cv <- qt(1 - (1 - coverage)/2, rdf)
            if (bands %in% c("uniform","both")) {
                E <- eigen(as.matrix(V))
                B <- E$vectors %*% diag(sqrt(pmax(0,E$values))) %*% t(E$vectors)
                D <- as.matrix(D)
                BX1 <- B %*% t(D[-1, ])
                BX1 <- BX1/sqrt(apply(BX1^2, 2, sum))
                BX0 <- B %*% t(D[-nrow(D), ])
                BX0 <- BX0/sqrt(apply(BX0^2, 2, sum))
                kappa <- sum(sqrt(apply((BX1 - BX0)^2, 2, sum)))
                cvu <- critval(kappa, alpha = 1 - coverage, rdf = rdf)
                if(bands == "both") cv <- c(cvu,cv)
                else cv <- cvu
                }
            list(pred = data.frame(x = G$x, y = G$y, se = se), 
                cv = cv)
        }
    }
    if(!length(select)) 
        select <- 1:m
    for (i in select) {
        qss <- x$qss[[i]]$xyz
        if (ncol(qss) == 3) {
            qss[, 3] <- x$coef[1] + qss[, 3]
            plot.qss2(qss, ...)
        }
        else if (ncol(qss) == 2) {
            if (length(bands)) {
                if (is.na(x$coef["(Intercept)"])) 
                  stop("rqss confidence bands require an intercept parameter")
                B <- summary(x$qss[[i]], V[[i]], ...)
                cv <- B$cv
                B <- B$pred
                B$y <- B$y + x$coef["(Intercept)"]
                if(!length(bcol)) bcol <- c("grey85","grey65")
		for(k in 1:length(cv)){
                 if (add || k > 1) 
                  if(shade){
                     polygon(c(B$x,rev(B$x)),
                        c(B$y - cv[k] * B$se,rev(B$y + cv[k] * B$se)),
			col = bcol[k], border = FALSE)
			}
                  else
                    matlines(B$x, cbind(B$y, B$y + cv[k] * cbind(-B$se, 
                        B$se)), lty = c(1, 2, 2), col = c("black", "blue", "blue"))
                else {
                  matplot(B$x, B$y + cv[k] * cbind(-B$se, B$se), 
                    xlab = paste(qssnames[i]), ylab = "Effect", 
                    type = "n", ...)
                  if(shade){
                     polygon(c(B$x,rev(B$x)),
                        c(B$y - cv[k] * B$se,rev(B$y + cv[k] * B$se)),
			col = bcol[k], border = FALSE)
			}
                  else{
                     lines(B$x, B$y + cv * B$se, lty = 2, ...)
                     lines(B$x, B$y - cv * B$se, lty = 2, ...)
		     }
	           }
                  lines(B$x, B$y, ...)
                }
                band[[1]] <- list(x = B$x, blo = B$y - B$se %o% cv, 
                  bhi = B$y + B$se %o% cv)
                if (rug) {
                  if (jit) 
                    rug(jitter(qss[, 1]))
                  else rug(qss[, 1])
                }
            }
            else {
                qss[, 2] <- x$coef[1] + qss[, 2]
                plot.qss1(qss, xlab = paste(qssnames[i]), ylab = "Effect", 
                  rug = rug, jit = jit, add = add, ...)
            }
            title(titles[i])
        }
        else stop("invalid fitted qss object")
    }
    if (pages > 0) 
        par(oldpar)
    if (length(bands)) 
        class(band) <- "rqssband"
    else
        band <- NULL
    invisible(band)
}

"triogram.fidelity" <- function (x, y, ndum=0, dummies = NULL)
{
#Make fidelity block of the triogram design in sparse matrix.csr form
#The rather esoteric match call identifies and handles duplicated xy points
n <- length(x)
A <- as.data.frame(cbind(x,y))
dupA <- duplicated(A)
if(any(dupA)){
  x <- x[!dupA]
  y <- y[!dupA]
  J <- match(do.call("paste",c(A,"\r")),do.call("paste",c(A[!dupA,],"\r")))
  z <- new("matrix.csr",ra=rep(1,n), ja=J, ia=1:(n+1),dimension=as.integer(c(n,max(J))))
  }
else{
  z <- as(n,"matrix.diag.csr")
  z <- as(z,"matrix.csr")
 }
#Augment with dummy vertices, if any...
if(length(dummies)){
	if (is.list(dummies)){
        	if (all(!is.na(match(c("x", "y"), names(dummies))))){
        		ndum <- length(dummies$x)
			if(length(dummies$y) == ndum){
        			x <- c(x,dummies$x)
        			y <- c(y,dummies$y)
        			zdum <- as.matrix.csr(0,n,ndum)
        			z <- cbind(z,zdum)
				}
        		else stop("dummies x and y components differ in length")
			}
        	else stop("dummies list lacking x and y elements")
		}
    	else stop("dummies argument invalid (not a list) in triogram.fidelity")
	}
else if(ndum > 0){
        u <- runif(ndum); v <- runif(ndum)
        xd <- min(x) + u * (max(x)-min(x))
        yd <- min(y) + v * (max(y)-min(y))
        T <- tripack::tri.mesh(x,y)
        s <- tripack::in.convex.hull(T,xd,yd)
        x <- c(x,xd[s])
        y <- c(y,yd[s])
        ndum <- sum(s)
        zdum <- as.matrix.csr(0,n,ndum)
        z <- cbind(z,zdum)
	dummies <- list(x = xd[s],y = yd[s])
        }
list(x=x,y=y,F=z, dummies = dummies)
}
"triogram.penalty" <- function (x, y, eps = .Machine$double.eps)
{
    n <- length(x)
    tri <- tripack::tri.mesh(x, y)
    bnd <- tripack::on.convex.hull(tri,x,y)
    q <- length(tri$tlist)
    m <- 13 * n
    z <- .Fortran("penalty", as.integer(n), as.integer(m), as.integer(q),
        as.double(x), as.double(y), as.integer(bnd),as.integer(tri$tlist),
        as.integer(tri$tlptr), as.integer(tri$tlend), rax = double(m),
	jax = integer(m), ned = integer(1), as.double(eps), ierr = integer(1),
	PACKAGE = "quantreg")[c("rax", "jax", "iax", "ned", "ierr")]
    if (z$ierr == 1)
        stop("collinearity in ggap")
    nnz <- 4 * z$ned
    ra <- z$rax[1:nnz]
    ja <- z$jax[1:nnz]
    ia <- as.integer(1 + 4 * (0:z$ned))
    dim <- as.integer(c(z$ned, n))
    new("matrix.csr",ra=ra,ja=ja,ia=ia,dimension=dim)
}

predict.rqss <-
function (object, newdata, interval = "none",  level = 0.95, ...) 
{
    ff <- object$fake.formula
    Terms <- delete.response(terms(object$formula, "qss"))
    Names <- all.vars(parse(text = ff))
    if (any(!(Names %in% names(newdata)))) 
        stop("newdata doesn't include some model variables")
    ff <- reformulate(ff)
    nd <- eval(model.frame(ff, data = newdata), parent.frame())
    qssterms <- attr(Terms, "specials")$qss
    if (length(qssterms)) {
        tmp <- untangle.specials(Terms, "qss")
        dropv <- tmp$terms
        m <- length(dropv)
        if (length(dropv)) 
            PLTerms <- Terms[-dropv]
        attr(PLTerms, "specials") <- tmp$vars
    }
    else {
        PLTerms <- Terms
        m <- 0
    }
    if(requireNamespace("MatrixModels") && requireNamespace("Matrix"))
        X <- as(MatrixModels::model.Matrix(PLTerms, data = nd, sparse = TRUE),"matrix.csr")
    else
        X <- model.matrix(PLTerms, data = nd)
    p <- ncol(X)
    y <- X %*% object$coef[1:p]
    X <- as.matrix.csr(X)
    if (m > 0) {
        for (i in 1:m) {
            qss <- object$qss[[i]]
            names <- all.vars(Terms[dropv[i]])
            names <- names[names %in% Names]
            dimnames(qss$xyz)[[2]] <- c(names, "zfit")
            newd <- nd[names]
            if (ncol(qss$xyz) == 3) {
                g <- predict.qss2(qss$xyz, newdata = newd, ...)
                y <- y + g$z
                if(interval == "confidence")
                    X <- cbind(X,g$D)
            }
            else if (ncol(qss$xyz) == 2) {
                g <- predict(qss, newdata = newd, ...)
                y <- y + g$y
                if(interval == "confidence")
                    X <- cbind(X,g$D)
            }
            else stop("invalid fitted qss object")
        }
    }
    if(interval == "confidence"){
        v <- sqrt(diag(X %*% summary(object, cov = TRUE)$V %*% t(X)))
        calpha <- qnorm(1 - (1-level)/2)
        y <- cbind(y,y - v*calpha,y + v*calpha)
        dimnames(y)[[2]] <- c("yhat","ylower","yupper")
	}
    y
}

"predict.qss1" <-
function (object, newdata, ...)
{
    x <- object$xyz[, 1]
    y <- object$xyz[, 2] 
    if(ncol(newdata)==1)
        newdata <- newdata[,1]
    else
        stop("newdata should have only one column for predict.qss1")
    if (any(diff(x) < 0))
        stop("x coordinates in qss1 object not monotone")
    if (max(newdata) > max(x) || min(newdata) < min(x))
        stop("no extrapolation allowed in predict.qss")
    bin <- cut(newdata, unique(x), label = FALSE, include.lowest = TRUE)
    p <- length(x)
    m <- length(newdata)
    V <- cbind(bin, bin + 1)
    B <- cbind(x[bin + 1] - newdata, newdata - x[bin])/(x[bin +
        1] - x[bin])
    ra <- c(t(B))
    ja <- as.integer(c(t(V)))
    ia <- as.integer(c(2 * (1:m) - 1, 2 * m + 1))
    dim <- c(m, p)
    D <- new("matrix.csr", ra = ra, ja = ja, ia = ia, dimension = dim)
    list(x = newdata, y = D %*% y, D = D[, -1])
}

predict.qss2 <- function (object, newdata, ...) 
{
    x <- object[, 1]
    y <- object[, 2]
    z <- object[, 3]
    tri.area <- function(v) {
        0.5 * ((v[2, 1] - v[1, 1]) * (v[3, 2] - v[1, 2]) - (v[3, 
            1] - v[1, 1]) * (v[2, 2] - v[1, 2]))
    }
    barycentric <- function(v) {
        b <- rep(0, 3)
        Area <- tri.area(v[1:3, ])
        b[1] <- tri.area(v[c(4, 2, 3), ])/Area
        b[2] <- tri.area(v[c(1, 4, 3), ])/Area
        b[3] <- tri.area(v[c(1, 2, 4), ])/Area
        if (any(b < 0 || b > 1)) 
            stop("barycentric snafu")
        b
    }
    if (is.list(newdata)) {
	fnames <- (dimnames(object)[[2]])[1:2]
        if (all(!is.na(match(fnames, names(newdata))))) {
            newx <- newdata[[fnames[1]]]
            newy <- newdata[[fnames[2]]]
	   }
        else (stop("qss object and newdata frame names conflict"))
        }
    else if (is.matrix(newdata)) 
        if (ncol(newdata) == 2) {
            newx <- newdata[, 1]
            newy <- newdata[, 2]
        }
        else (stop("newdata matrix must have 2 columns"))
    tri <- tripack::tri.mesh(x, y)
    if (!all(tripack::in.convex.hull(tri, newx, newy))) 
        stop("some newdata points outside convex hull")
    p <- length(x)
    m <- length(newx)
    V <- matrix(0, m, 3)
    B <- matrix(0, m, 3)
    for (i in 1:m) {
        V[i, ] <- unlist(tripack::tri.find(tri, newx[i], newy[i]))
        v <- rbind(cbind(x[V[i, ]], y[V[i, ]]), c(newx[i], newy[i]))
        B[i, ] <- barycentric(v)
    }
    ra <- c(t(B))
    ja <- as.integer(c(t(V)))
    ia <- as.integer(3 * (0:m) + 1)
    D <- new("matrix.csr", ra = ra, ja = ja, ia = ia, dimension = c(m, 
        p))
    list(x = newx, y = newy, z = c(D %*% z), D = D[,-1])
}

fitted.rqss  <-  function(object, ...) (object$X %*% object$coef)[1:object$n]
resid.rqss <- function(object, ...) object$resid[1:object$n]

"rqss" <- function (formula, tau = 0.5, data = parent.frame(), weights, 
    na.action, method = "sfn", lambda = NULL, contrasts = NULL, 
    ztol = 1e-05,  control = sfn.control(), ...) 
{
    call <- match.call()
    m <- match.call(expand.dots = FALSE)
    temp <- c("", "formula", "data", "weights", "na.action")
    m <- m[match(temp, names(m), nomatch = 0)]
    m[[1]] <- as.name("model.frame")
    special <- "qss"
    Terms <- if (missing(data)) 
        terms(formula, special)
    else terms(formula, special, data = data)
    qssterms <- attr(Terms, "specials")$qss
    dropx <- NULL
    if(length(tau) > 1){
        tau <- tau[1]
        warning("multiple taus not supported, using first element")
        }
    if (length(qssterms)) {
        tmpc <- untangle.specials(Terms, "qss")
        ord <- attr(Terms, "order")[tmpc$terms]
        if (any(ord > 1)) 
            stop("qss can not be used in an interaction")
        dropx <- tmpc$terms
        if (length(dropx)) 
            Terms <- Terms[-dropx]
        attr(Terms, "specials") <- tmpc$vars
        fnames <- function(x){
           fy <- all.names(x[[2]])
           if(fy[1] == "cbind") fy <- fy[-1]
           fy
           }
        fqssnames <- unlist(lapply(parse(text = tmpc$vars), fnames))
        qssnames <- unlist(lapply(parse(text = tmpc$vars), function(x) deparse(x[[2]])))
    }
    m$formula <- Terms
    ff <- delete.response(terms(formula(Terms)))
    if(exists("fqssnames")){
       ffqss <- paste(fqssnames, collapse = "+")
       ff <- paste(deparse(formula(ff)), "+", ffqss)
    }
    m <- eval(m, parent.frame())
    weights <- model.extract(m, weights)
    process <- (tau < 0 || tau > 1)
    Y <- model.extract(m, "response")
    if(requireNamespace("MatrixModels") && requireNamespace("Matrix")){
         X <- MatrixModels::model.Matrix(Terms, m, contrasts, sparse = TRUE)
         vnames <- dimnames(X)[[2]]
         X <- as(X ,"matrix.csr")
         }
    else{
         X <- model.matrix(Terms, m, contrasts)
         vnames <- dimnames(X)[[2]]
         }
    p <- ncol(X)
    pf <- environment(formula)
    nrL <- 0
    if (method == "lasso") {
        if (!length(lambda)) 
            stop("No lambda specified for lasso constraint")
        if (length(lambda) == 1) 
            lambda <- c(0, rep(lambda, p - 1))
        if (length(lambda) != p) 
            stop("lambda must be either of length p, or length 1")
        if (any(lambda < 0)) 
            stop("negative lambdas disallowed")
        L <- diag(lambda, nrow = length(lambda))
        L <- L[which(lambda != 0), , drop = FALSE]
        L <- as.matrix.csr(L)
        nrL <- nrow(L)
        ncL <- ncol(L)
    }
    if (length(qssterms) > 0) {
        F <- as.matrix.csr(X)
        qss <- lapply(tmpc$vars, function(u) eval(parse(text = u), 
            data, enclos = pf))
        mqss <- length(qss)
        ncA <- rep(0, mqss + 1)
        nrA <- rep(0, mqss + 1)
        nrR <- rep(0, mqss + 1)
        for (i in 1:mqss) {
            F <- cbind(F, qss[[i]]$F)
            ncA[i + 1] <- ncol(qss[[i]]$A)
            nrA[i + 1] <- nrow(qss[[i]]$A)
            nrR[i + 1] <- ifelse(is.null(nrow(qss[[i]]$R)), 0, 
                nrow(qss[[i]]$R))
            vnames <- c(vnames, paste(qssnames[i], 1:ncA[i + 
                1], sep = ""))
        }
        A <- as.matrix.csr(0, sum(nrA), sum(ncA))
        if (sum(nrR) > 0) {
            R <- as.matrix.csr(0, sum(nrR), sum(ncA))
            nrR <- cumsum(nrR)
        }
        ncA <- cumsum(ncA)
        nrA <- cumsum(nrA)
        lambdas <- rep(0, mqss)
        for (i in 1:mqss) {
            lambdas[i] <- qss[[i]]$lambda
            Arows <- (1 + nrA[i]):nrA[i + 1]
            Acols <- (1 + ncA[i]):ncA[i + 1]
            A[Arows, Acols] <- qss[[i]]$lambda * qss[[i]]$A
            if (nrR[i] < nrR[i + 1]) 
                R[(1 + nrR[i]):nrR[i + 1], (1 + ncA[i]):ncA[i + 
                  1]] <- qss[[i]]$R
        }
        A <- cbind(as.matrix.csr(0, nrA[mqss + 1], p), A)
        if (nrR[mqss + 1] > 0) {
            R <- cbind(as.matrix.csr(0, nrR[mqss + 1], p), R)
            r <- rep(0, nrR[mqss + 1])
        }
        else {
            R <- NULL
            r <- NULL
        }
        if (method == "lasso") {
            A <- rbind(cbind(L, as.matrix.csr(0, nrL, ncol(F) - ncL)), A)
        }
        X <- rbind(F, A)
        Y <- c(Y, rep(0, nrow(A)))
        rhs <- t(rbind((1 - tau) * F, 0.5 * A)) %*% rep(1, nrow(X))
        XpX <- t(X) %*% X
        nnzdmax <- XpX@ia[length(XpX@ia)] - 1
	if(is.null(control[["nsubmax"]]))
           control[["nsubmax"]] <- max(nnzdmax, floor(1000 + exp(-1.6) * nnzdmax^1.2))
        if(is.null(control[["nnzlmax"]]))
           control[["nnzlmax"]] <- floor(2e+05 - 2.8 * nnzdmax + 7e-04 * nnzdmax^2)
        if(is.null(control[["tmpmax"]]))
           control[["tmpmax"]] <- floor(1e+05 + exp(-12.1) * nnzdmax^2.35)
        fit <- if (length(r) > 0) 
            rqss.fit(X, Y, tau = tau, rhs = rhs, method = "sfnc", 
                R = R, r = r, control = control, ...)
        else rqss.fit(X, Y, tau = tau, rhs = rhs, method = "sfn", control = control, ...)
        for (i in 1:mqss) {
            ML <- p + 1 + ncA[i]
            MU <- p + ncA[i + 1]
            qss[[i]] <- list(xyz = cbind(qss[[i]]$x$x, qss[[i]]$x$y, 
                c(0, fit$coef[ML:MU])), dummies = qss[[i]]$dummies)
            if (ncol(qss[[i]]$xyz) == 2) 
                class(qss[[i]]) <- "qss1"
            else class(qss[[i]]) <- "qss2"
        }
        names(qss) <- qssnames
        fit$qss <- qss
    }
    else {
        X <- as.matrix.csr(X)
        nrA <- 0
        if (method == "lasso") { #Pure lasso case
            rhs <- t(rbind((1 - tau) * X, 0.5 * L)) %*% rep(1, nrow(X) + nrow(L))
            X <- rbind(X, L)
            Y <- c(Y, rep(0, nrL))
            #nrA <- c(nrA, nrL) # Why was this here?
        }
        else
            rhs <- NULL
        if (length(weights)) {
            if (any(weights < 0)) 
                stop("negative weights not allowed")
            X <- X * weights
            Y <- Y * weights
        }
        fit <- rqss.fit(X, Y, tau = tau,  rhs = rhs, control = control, ...)
        fit$nrA <- nrA
    }
    names(fit$coef) <- vnames
    n <- length(fit$resid) - nrL - nrA[length(nrA)]
    uhat <- fit$resid[1:n]
    Rho <- function(u, tau) sum(u * (tau - (u < 0)))
    fit$fidelity <- Rho(uhat, tau)
    fit$edf <- sum(abs(uhat) < ztol)
    fit$X <- X
    fit$n <- n
    fit$nrL <- nrL
    fit$terms <- Terms
    fit$fake.formula <- ff
    fit$formula <- formula
    fit$method <- method
    fit$call <- call
    fit$tau <- tau
    if (length(qssterms)) {
        fit$lambdas <- lambdas
        fit$qssnames <- qssnames
        fit$nrA <- nrA
        fit$ncA <- cumsum(c(p, diff(ncA)))
    }
    else fit$ncA <- p
    attr(fit, "na.message") <- attr(m, "na.message")
    class(fit) <- "rqss"
    fit
}


"summary.rqss" <- function(object, cov = FALSE, ztol = 1e-5, ...){
    resid <- object$resid
    coef <- object$coef
    lambdas <- object$lambdas
    formula <- object$formula
    fidelity <- object$fidelity
    edf <- object$edf
    nrA <- object$nrA
    ncA <- object$ncA
    nrL <- object$nrL
    tau <- object$tau
    X <- object$X
    p <- ncol(object$X)
    m <- length(ncA)
    n <- length(resid) - nrL - nrA[m]
    uhat <- resid[1:n]
    h <- bandwidth.rq(tau, n, hs = TRUE)
    if(tau + h > 1)
       stop("tau + h > 1:  error in summary.rq")
    if(tau - h < 0)
       stop("tau - h < 0:  error in summary.rq")
    h <- (qnorm(tau + h) - qnorm(tau - h)) * max(ztol,min(sqrt(var(uhat)),
        (quantile(uhat, 0.75) - quantile(uhat, 0.25))/1.34))
    f <- c(dnorm(uhat/h)/h,rep(1, length(resid) - n))
    D <- t(X) %*% (f * X)
    D <- chol(.5 * (D + t(D)), ...)
    D <- backsolve(D,diag(p))
    D0 <- tau * (1 - tau) * t(X) %*% X
    V <- D %*% D0 %*% D
    scale <- mean(f)
    serr <- sqrt(diag(V))
    ptab <- array(coef[1:ncA[1]], c(ncA[1],4))
    pnames <- names(coef[1:ncA[1]])
    dimnames(ptab) <- list(pnames, c("Estimate", "Std. Error", "t value", "Pr(>|t|)"))
    ptab[, 2] <- serr[1:ncA[1]]
    ptab[, 3] <- ptab[, 1]/ptab[, 2]
    ptab[, 4] <- 2 * (1 - pt(abs(ptab[, 3]), n - edf))
    if(cov) {
        Vcov <- V[1:ncA[1],1:ncA[1]]
        Vqss <- as.list(1:(m-1))
        }
    if(m > 1) {
        penalty <- rep(NA, m - 1)
        qssedfs <- rep(NA, m - 1)
        ntab <- matrix(0, m - 1, 5)
        vnames <- sub(".$","",names(coef[ncA[-length(ncA)]+1]))
        dimnames(ntab) <- list(vnames, c("EDF", "Lambda", "Penalty", "F value", "Pr(>F)"))
        for(i in 2:m) {
           ntab[i - 1, 3] <- sum(abs(resid[n + nrL + ((nrA[i - 1] + 1):nrA[i])]))
           ntab[i - 1, 1] <- sum(abs(resid[n + nrL + ((nrA[i - 1] + 1):nrA[i])]) > ztol)
           v <- V[c(1, (ncA[i-1] + 1):ncA[i]), c(1, (ncA[i-1] + 1):ncA[i])]
           v <- .5 * (v + t(v))
           if(cov) Vqss[[i-1]] <- v
           b <- coef[(ncA[i-1] + 1):ncA[i]]
           ntab[i-1, 4]  <- t(b) %*% solve(v[-1,-1],b)/ntab[i-1,1]
           ntab[i-1, 5]  <- 1 - pf(ntab[i-1,4],ntab[i-1,1], n - edf)
           }
        ntab[,2] <- lambdas
        ntab[,3] <- ntab[,3]/lambdas
        if(cov)
           z <- list(coef = ptab, qsstab = ntab, fidelity = fidelity, tau = tau,
		formula = formula, edf = edf, n = n, Vcov = Vcov, Vqss = Vqss, V = V)
        else
           z <- list(coef = ptab, qsstab = ntab, fidelity = fidelity, tau = tau,
		formula = formula, edf = edf, n = n)
        }
    else
        if(cov)
           z <- list(coef = ptab, fidelity = fidelity, formula = formula, tau = tau,
		edf = edf, n = n, Vcov = Vcov, V = V)
        else
           z <- list(coef = ptab, fidelity = fidelity, formula = formula, tau = tau,
		edf = edf, n = n)

    class(z) <- "summary.rqss"
    return(z)
}
"plot.summary.rqss" <- function(x, ...)
    warning("No plot method for summary.rqss objects:   plot the rqss object instead")

"print.rqss" <- function(x, ...) {
    cat("Formula:\n")
    print(x$formula)
    sx <- summary(x)
    cat("Quantile fidelity at tau = ", x$tau, "is", sx$fidelity, "\n")
    cat("Estimated Model Dimension is", sx$edf, "\n")
}
 
"logLik.rqss" <- function(object, ...){
	n <- object$n
	tau <- object$tau
	val <- n * (log(tau * (1-tau)) - 1 - log(object$fidelity/n))
	attr(val,"n") <- n
	attr(val,"df") <- object$edf
	class(val) <- "logLik"
	val
	}
"AIC.rqss" <- function(object, ... , k = 2){
	v <- logLik(object)
	if(k < 0) 
		k <- log(attr(v,"n"))
	val <- AIC(logLik(object), k = k)
	attr(val,"edf") <- attr(v,"df")
	val
	}
"dither"  <- function(x, type = "symmetric", value = NULL) {
	if(length(x) == 0)
		return(x)
	if(!is.numeric(x))
		stop("'x' must be numeric")
	if(!length(value))
		value <- min(diff(sort(unique(x))))
	if(type == "symmetric")
		v <- x + runif(length(x), -value/2, value/2)
	else if(type == "right")
		v <- x + runif(length(x), 0, value)
	else
		stop("invalid type")
	v
}
	
print.summary.rqss <-
function (x, digits = max(3, getOption("digits") - 3), signif.stars = getOption("show.signif.stars"), 
    ...) 
{
    cat("Formula:\n")
    print(x$formula)
    if (length(x$coef) > 0) {
        cat("\nParametric coefficients:\n")
        printCoefmat(x$coef, digits = digits, signif.stars = signif.stars, 
            na.print = "NA", ...)
    }
    cat("\n")
    if (length(x$qsstab) > 0) {
        cat("Approximate significance of qss terms:\n")
        printCoefmat(x$qsstab, digits = digits, signif.stars = signif.stars, 
            has.Pvalue = TRUE, na.print = "NA", cs.ind = 1, ...)
    }
    cat("\n")
    if (length(x$fidelity) > 0) 
        cat("  Quantile Fidelity at tau = ", x$tau, "  is  ", formatC(x$fidelity, 
		digits = 6, width = 11), "\n", sep = "")
    cat("  Effective Degrees of Freedom = ", formatC(x$edf, digits = 5, width = 8, 
        flag = "-"), "  Sample Size = ", x$n, "\n", sep = "")
    invisible(x)
}
critval <- function(kappa, alpha = 0.05, rdf = 0){
# Hotelling tube critical value for uniform confidence bands
# This should agree to about 6 digits with the following call to locfit:
# crit(const = c(kappa,1), cov = 1-alpha,rdf = rdf)$crit.val
   tube <- function(x,alpha,rdf){
      if(rdf <= 0)
         kappa * exp(-x^2/2)/pi + 2 * (1 - pnorm(x)) -  alpha
      else
         kappa * (1+x^2/rdf)^(-rdf/2)/pi + 2 * (1 - pt(x,rdf)) -  alpha
      }
   uniroot(tube,c(1,5),alpha = alpha, rdf = rdf)$root
   }