crq.fit.pow <- function(x, y, yc, tau=0.5, weights=NULL, start, left=TRUE,maxit = 500){

 x <- as.matrix(x)
 y <- as.vector(y)
 n <- nrow(x)
 p <- ncol(x)
 if(missing(yc)) {
     yc <- rep(0,n) 
     if(!left) warning("default right (!?) censoring at zero")
    } 
 yc <- as.vector(yc)
 if(length(weights)){
     if (any(weights < 0))
         stop("negative weights not allowed")
     contr <- attr(x, "contrasts")
     x <- x * weights
     y <- y * weights
     }
 if(left) {
	x <- -x
	y <- -y
	yc <- -yc
	tau <- 1-tau
	}
# Starting value of beta 
if(missing(start)){ 
   tol <- .Machine$double.eps^(2/3)
   r <- rq.fit.br(x,y,tau)$residuals
   h <- which(abs(r) < tol)[1:p]
   }
else if(start[1] == "global"){
# Generate H matrix 
   H <- combos(n,p)
   m <- ncol(H)

# Starting value of beta 
   U <- solve(x[H[,1],])
   bh <- U %*% y[H[,1]]

f <- .Fortran("brutpow",
        as.integer(n),
        as.integer(p),
        as.integer(m),
        as.integer(H),
        as.double(x),
        as.double(y),
        as.double(yc),
        as.double(bh),
        as.double(tau),
        as.double(U),
        double(p),
        double(p),
        kminz = integer(1),
        nflag = as.integer(0),
        PACKAGE = "quantreg")
if(f$nflag!=0)
        warning(switch(f$nflag,
                "Error in pivot:  hout not in h",
                "Error in pivot:  hin in h",
                "Error in pivot:  hin out of bounds",
                "Error in findk:  k not found"
                ))
   k  <- f$kminz
   bh <- solve(x[H[,k],],y[H[,k]])
   residuals <- as.matrix(y - pmin(yc, x %*% bh))

   return(list(coefficients= bh,residuals = residuals))
  }
else if(is.numeric(start) && length(start) == p){
	if(all(start %in% 1:n)) h <- start
	else  h <- order(abs(y - x %*% start))[1:p]
	}
else
   stop("Invalid starting value")

xhinv <- try(solve(x[h,]))
if(class(xhinv) == "try-error")
	stop("Singular basic solution generated by 'start'")

f <- .Fortran("powell", 
	as.integer(n), 
	as.integer(p), 
	as.integer(2*p), 
	as.double(x),
	as.double(y), 
	as.double(yc), 
	coef = double(p), 
	as.double(tau), 
	as.integer(h), 
	double(n), 
	as.double(xhinv), 
	double(n), 
	double(2*p),
	double(p),
	double(p),
	as.integer(maxit), 
	nflag = as.integer(0), 
	PACKAGE = "quantreg")
if(f$nflag!=0) 
	warning(switch(f$nflag, "Max iterations reached",
		"Solution may be nonunique",
		"Error in pivot:  hout not in h",
		"Error in pivot:  hin in h",
		"Error in pivot:  hin out of bounds"
		))
coef <- f$coef
residuals <- as.matrix(y - pmin(yc, x %*% coef))

return (list(coefficients= coef, residuals = residuals))
}


Curv <- function (y,  yc, ctype = c("left", "right")) 
{
    nn <- length(y)
    if(length(yc) != nn)
        stop("Event times and censoring times of different length")
    ctype <- match.arg(ctype)
    if(ctype == "right" && any(y > yc))
        stop("Event times can not exceed ctimes for right censoring")
    if(ctype == "left" && any(y < yc))
        stop("Event times can not be less than ctimes for left censoring")
    if(!(ctype == "left" || ctype == "right"))
        stop("Invalid ctype for method Powell")
    ss <- cbind(y, yc)
    dimnames(ss) <- list(NULL, c("time", "ctime"))
    attr(ss, "ctype") <- ctype
    class(ss) <- "Surv"
    ss
}

QTECox <- function(x, smooth = TRUE){
# compute quantile treatment effect for a Cox PengHuang fit
g <- survfit(x)
if(smooth)
        g <- supsmu(1-g$surv,g$time)
taus <- (g$x[-1] + g$x[-length(g$x)])/2
Qhat <- (g$y[-1] + g$y[-length(g$y)])/2

dQ <- diff(Qhat)/diff(taus)
taus <- taus[-1]; Qhat <- Qhat[-1]
QTE <- outer(dQ * (1-taus) * log(1-taus) / Qhat, coef(x))
list(QTE = QTE , taus = taus)
}
boot.crq <- function(x, y, c, taus, method, ctype = "right", R=100,  
		mboot,  bmethod = "Bose", ...)
        {
        n <- length(y)
        p <- ncol(x)
        if(missing(mboot)) mboot <- n
        A <- array(0,dim=c(p,length(taus),R))

        for (i in 1:R){
                if(bmethod == "xy-pair"){
                        w <- table(sample(1:n,mboot,replace=TRUE))
                        s <- as.numeric(names(w))
                        w <- as.numeric(w)
                        yb <- y[s]
                        xb <- x[s,]
                        cb <- c[s]
                        }
                else if(bmethod == "Bose"){
                        w <- rexp(n)
                        yb <- y
                        xb <- x
                        cb <- c
                        }
                else
                        stop("invalid bmethod for boot.crq")
		if(method == "Portnoy")
                	a <- crq.fit.por(xb,yb,cb, weights = w, ctype = ctype, ... )
		else if(method == "PengHuang")
                	a <- crq.fit.pen(xb,yb,cb, weights = w, ctype = ctype, ... )
		else
			stop("Invalid method for boot.crq")
                if((i %% floor(R/10)) == 0 )
                        cat(paste("bootstrap roughly ",100*(i/R)," percent complete\n"))
                A[,,i] <- coef(a,taus)
                }
        list(A = A, n = length(y), mboot = mboot)
        }

coef.crq <- function(object, taus = 1:4/5, ...)
        {
        # Extract coefficients from the crq solution array 
        
        if(min(taus) < 0 || max(taus) > 1) stop("taus out of range [0,1]")
	if(length(object$tau) == 1){
        	coef <- object$coefficients
        	return(coef)
        	}
        taus <- sort(taus)
        S <- object$sol
        r <- S[1, ]
        r <- c(r[1],r)
        r <- (r[-1]+r[-length(r)])/2
	r <- c(0,r)
	if(is.unsorted(r)) r <- r[-length(r)] #kludge until why this happens is found
        B <- S[-1,]
        J <- length(r)
	B <- t(cbind(B[,1],B))
        ts <- taus[taus < max(r)]
        bin <- findInterval(ts,r)
        wgt <- (ts - r[bin])/(r[bin + 1] - r[bin])
        binlag <- bin - 1
        binlag[binlag == 0] <- 1
        coef <- t(wgt * B[bin, , drop = FALSE] + (1 - wgt) * B[binlag, , drop = FALSE])
        nna <- length(taus) - length(ts)
        if(nna > 0)
                coef <- cbind(coef, matrix(NA,nrow(coef),nna))
        taulabs <- paste("tau=", format(round(taus, 3)))
        dimnames(coef)[[2]] <- taulabs
        coef[-nrow(coef),]  # Delete Qhat entries
        }

crq <- function (formula, taus, data, subset, weights, na.action, 
	method = c("Powell","Portnoy","PengHuang"), contrasts = NULL, ...)
{
    require(survival)
    call <- match.call()
    mf <- match.call(expand.dots = FALSE)
    m <- match(c("formula", "data", "subset", "weights", "na.action"),
        names(mf), 0)
    mf <- mf[c(1, m)]
    mf$drop.unused.levels <- TRUE
    mf[[1]] <- as.name("model.frame")
    mf <- eval.parent(mf)
    if (method == "model.frame")
        return(mf)
    mt <- attr(mf, "terms")
    X <- model.matrix(mt, mf, contrasts)
    weights <- as.vector(model.weights(mf))
    Y <- model.extract(mf, "response")
    eps <- .Machine$double.eps^(2/3)
    if(!inherits(Y,"Surv"))
        stop("Response must be a survival object")
    method <- match.arg(method)
    if(method == "Powell") {
	ctype <- attr(Y, "ctype")
	if(!ctype %in% c("right","left")) 
             stop("Only right or left censoring Surv objects are allowed")
	left <- (ctype == "left")
        if (any(taus < -eps) || any(taus > 1 + eps))
            stop("invalid taus:  taus should be >= 0 and <= 1")
	y <- Y[,1]
	cen <- Y[,2]
	if (length(taus) > 1) {
          coef <- matrix(0, ncol(X), length(taus))
          fitted <- resid <- matrix(0, nrow(X), length(taus))
          for (i in 1:length(taus)) {
            z <- crq.fit.pow(X, y, cen, tau = taus[i], weights, left = left, ...)
            coef[, i] <- z$coefficients
            resid[, i] <- z$residuals
            fitted[, i] <- y - z$residuals
          }
          taulabs <- paste("tau=", format(round(taus, 3)))
          dimnames(coef) <- list(dimnames(X)[[2]], taulabs)
          dimnames(resid) <- list(dimnames(X)[[1]], taulabs)
          fit <- list(coefficients = coef, residuals = resid, fitted.values = fitted)
	  fit$tau <- taus
	  class(fit) <- "crqs"
	}
	else {
          fit <- crq.fit.pow(X, y, cen, tau = taus, weights, left = left, ...)
	  fit$tau <- taus
	  class(fit) <- "crq"
	}
     }
    else if(method == "Portnoy"){
        ctype <- "right"
        if(attr(Y,"type") != "right"){
            if(attr(Y,"type") == "left") ctype <- "left"
            else stop("Only right censoring Surv objects are allowed for Portnoy method")
            }
        y <- Y[,1]
        cen <-  Y[,2]
        fit <- crq.fit.por(X, y, cen, weights, ctype = ctype,   ...)
	class(fit) <- "crq"
	}
    else if(method == "PengHuang"){
        ctype <- "right"
        if(attr(Y,"type") != "right"){
            if(attr(Y,"type") == "left") ctype <- "left"
            else stop("Only right censoring Surv objects are allowed for Portnoy method")
            }
        y <- Y[,1]
        cen <-  Y[,2]
        fit <- crq.fit.pen(X, y, cen, weights,  ctype = ctype, ...)
	class(fit) <- "crq"
	}
    else
	stop("Method not defined for crq")

fit$terms <- mt
fit$call <- call
fit$formula <-  formula(mt)
fit$method <-  method
fit$ctype <-  ctype
attr(fit, "na.message") <- attr(m, "na.message")
fit
}
predict.crqs <-
function (object, newdata, type = NULL, ...) {
pred <-  predict.rqs(object, newdata, ...) 
	}
predict.crq <-
function (object, newdata,  ...) {
method <- object$method
if(method %in% c("Portnoy","PengHuang")){ #Kludge to make crq sol matrix look like rq sol matrix 
	p2 <- nrow(object$sol)
	object$sol <- object$sol[c(1,p2,p2,2:(p2-1)),]
	}
pred <- switch(method,
	"Powell" = predict.rq(object, newdata,  ...), 
	"Portnoy"  = predict.rq.process(object, newdata,  ...), 
	"PengHuang" = predict.rq.process(object, newdata,  ...) 
	)
}
crq.fit.por <- function(x, y, cen, weights = NULL, grid, ctype = "right")
{
      p <- ncol(x)
      n <- length(y) 
      cen <- 1 - cen #NB: Fortran routine wants censoring indicator flipped (!!!)
      mp <- n + max(1, sum(cen))
      eps <- 1e-04
      kmax <- 10000 # perhaps this should be parameter in future
      if(length(weights)){
		if (any(weights < 0)) 
		        stop("negative weights not allowed")
    		contr <- attr(x, "contrasts")
    		x <- weights * x
    		y <- weights * y
		}
    if(ctype == "left") y <- -y
    if(missing(grid))  grid <-  seq(1/n,1-1/n,by=1/(5 + 3 * n^.4))
    if(is.numeric(grid)){
		ginit <- min(grid)
		dgrid <- diff(grid)
		gstep <- median(dgrid)
		if(any(dgrid <0)) stop("grid is not monotonic")
		if(gstep < eps) stop("grid stepsize too small")
		nsol <- 3*n
		mw = -1
		}
      else if(grid  == "pivot"){
		nsol <- 3*n
		ginit <- 1/(2*n)
		gstep <- 1/(2*n)
		mw <- 20
		}
      else 
		stop("Invalid grid")
	z <- .Fortran("crq",
		as.integer(n),
		as.integer(p),
                as.integer(mp), 
                as.integer(mp+5), 
		as.integer(p+2),
		as.double(x),
		as.double(y),
		as.integer(cen),
		as.double(ginit),
                as.integer(mw),
                as.double(gstep),
		ift = integer(1),
                as.integer(kmax),
		h = integer(p),
                xh = double(p*p),
		wa = double(mp*p),
		wb = double(mp),
		wc = double(mp*(p+2)),
		wd = double(mp),
		we = double(mp),
		wf = double(p),
		iflag = integer(mp),
                as.integer(nsol),
		sol = double(nsol*(p+2)),
                lsol = integer(1),
		icen = integer(n),
		tcen = double(n),
		lcen = integer(1),
		PACKAGE = "quantreg")
	nw <- z$h[1]
	flag <- z$ift
	msg <- switch(flag,
		paste("Error in input dimensions, n,p,mw "),
		paste("Error in input dimensions, n,p,mw "),
		paste("Error in input dimensions, n,p,mw "),
		paste("Less than p=",p,"observations above tau = 0 solution"),
		paste("Possible degeneracy at",nw,"tau values.",
          		"$tau.degen: first mp =", n + 5 + sum(cen)," such tau values"), 
		paste("Number of pivots to be saved in sol > nsol.",
          		"Redefine nsol: use nsol < n to save for tau = i/(nsol-1)"),
		paste("Error with partial return: possible degeneracies",
         		"Max number of rq calls exceeded: dither x or increase mw"),
		paste("Premature stop: defective conditional distribution"),
		paste("Simplex iteration limit exceeded -- consider dithering y"))
	if(flag > 0 && flag != 5 && flag < 10)
		ifelse(flag <= 3,stop(msg),warning(msg))
	J <- z$lsol
	B <- matrix(z$sol, nrow=p+2, ncol=nsol, byrow=FALSE)[,1:J]
	dimnames(B) <- list(c("tau",dimnames(x)[[2]],"Qbar"),NULL)
	ic <- z$icen
	sp <- (1:n)[ic == 1]
	tsp <- z$tcen[sp]
        t1 <- z$wd[1:nw]
        if(ctype == "left") {
            sol[1,] <- 1 - sol[1,]
            sol[-1,] <- - sol[-1,]
            }
	a<-list(sol=B, Isplit=sp, tsp = tsp,status=ic)
	class(a) <- "crq"
	return(a) 
	}
crq.fit.pen <- function(x, y, cen, weights=NULL,grid, ctype = "right" ){
      p <- ncol(x)
      n <- length(y)
      if(missing(grid)) grid <-  seq(1/n,1-1/n,by=min(0.01,1/(2*length(y)^.7))) 
      if(!is.numeric(grid)) stop("Invalid grid")
      if(any(grid < 0) || any(grid > 1)) stop("Invalid grid")
      m <- length(grid)
      xbar <- apply(x,2,mean)
      if(length(weights)){
                if (any(weights < 0))
                        stop("negative weights not allowed")
                contr <- attr(x, "contrasts")
                x <- x * weights
                y <- y * weights
                }
        if(ctype == "left") y <- -y
        s <- rep(0,n)
        u <- rep(1,n)
        d <- rep(1,n)
        r <- rep(1,p)
        B <- matrix(0,p,m)
        cc <- as.logical(cen)
        y1 <- y[cc]
	n1 <- length(y1)
        x1 <- x[cc,]
	z <- .Fortran("crqfnb", as.integer(n), as.integer(p), 
		a1 = as.double(t(as.matrix(x1))), c1 = as.double(-y1), n1=as.integer(n1),
		as.double(x), as.double(y),as.double(cen),B =as.double(B),
		g = as.double(grid),m = as.integer(m), as.double(r), 
		as.double(s), as.double(d), as.double(u),
        	wn = double(n1 * 9), wp = double((p + 3) * p),
        	info = integer(1), PACKAGE = "quantreg")
	J <- z$m - 1
	B <- matrix(-z$B, p, m)
	B <- B[,1:J]
        qhat <- t(xbar) %*% B
        B <- rbind(grid[1:J],B,qhat)
        dimnames(B) <- list(c("tau",dimnames(x)[[2]],"Qhat"),NULL)
        if(ctype == "left") {
            B[1,] <- 1 - B[1,]
            B[-1,] <- - B[-1,]
            }
        B  <- list(sol=B)
	class(B) <- "crq"
	B
        }
plot.summary.crqs <-
function (x, nrow = 3, ncol = 3, CoxPHit = NULL, ...) {
    taus <- function(x) x$tau
    xx <- unlist(lapply(x, taus))
    coef <- lapply(x, coefficients)
    p <- nrow(coef[[1]])
    k <- ncol(coef[[1]])
    if(k != 6) stop("summary.crqs object has wrong column dimension")
    m <- length(xx)
    blab <- dimnames(coef[[1]])[[1]]
    a <- array(unlist(coef), c(p, k, m))
    if(length(CoxPHit))
	CoxQTE <- QTECox(CoxPHit)
    oldpar <- par(no.readonly=TRUE)
    par(mfrow = c(nrow, ncol))
    for (i in 2:p) {
            b  <- a[i, 1, ]
            bl <- a[i, 2, ]
            bu <- a[i, 3, ]
        plot(rep(xx, 2), c(bl, bu), xlab = "", ylab = "", type = "n")
        title(paste(blab[i]), cex = 0.75)
        polygon(c(xx, rev(xx)), c(bl, rev(bu)), col = "LightSkyBlue")
        points(xx, b, cex = 0.5, pch = "o", col = "blue")
        lines(xx, b, col = "blue")
        abline(h = 0)
        if(length(CoxPHit)) {
	    lines(CoxQTE$taus,CoxQTE$QTE[,i-1],col="red")
            }
    }
par(oldpar)
}
summary.crq <-
function (object, taus = 1:4/5, alpha = .05, se = "boot", covariance = TRUE, ...) {
    mt <- terms(object)
    m <- model.frame(object)
    x <- model.matrix(mt, m, contrasts = object$contrasts)
    Y <- model.response(m)
    method <- object$method
    ctype <- object$ctype
    y <- Y[,1]
    cen  <- Y[,2]
    wt <- as.vector(model.weights(object$model))
    if (!is.null(wt)) {
        resid <- resid * wt
        x <- x * wt
        y <- y * wt
	cen <- cen * wt
    }
    if(method == "Powell"){
       coef <- coefficients(object)
       vnames <- dimnames(x)[[2]]
       resid <- object$residuals
       tau <- object$tau
       n <- length(resid)
       p <- length(coef)
       rdf <- n - p
       if (se == "boot") {
           if(attr(Y,"type") == "left")
              s <- cen  <  x %*% coef
           else
              s <- cen  >  x %*% coef
           B <- boot.rq(x[s, ], y[s], tau, ...)
           cov <- cov(B)
           serr <- sqrt(diag(cov))
           }
       else stop("Only boot method is implemented for crq inference")
       coef <- array(coef, c(p, 4))
       dimnames(coef) <- list(vnames, c("Value", "Std. Error", "t value", 
           "Pr(>|t|)"))
       coef[, 2] <- serr
       coef[, 3] <- coef[, 1]/coef[, 2]
       coef[, 4] <- if(rdf > 0) 
           2 * (1 - pt(abs(coef[, 3]),rdf)) 
       else NA
       object <- object[c("call", "terms")]
       if (covariance == TRUE) 
           object$cov <- cov
       object$B <- B
       object$coefficients <- coef
       object$rdf <- rdf
       object$tau <- tau
       class(object) <- "summary.crq"
       return(object)
       }
    else if(method == "Portnoy" || method == "PengHuang") {
       coef <- as.matrix(coef(object,taus))
       coef <- coef[,apply(coef,2,function(x) any(!is.na(x))),drop = FALSE] # Delete NA columns if any
       taus <- taus[1:ncol(coef)]
       B <- boot.crq(x, y, cen, taus, method = method, ctype = ctype, ...)
       nas <- apply(is.na(B$A[1,,]),1,sum)
       sqmn <- sqrt(B$mboot/B$n)
       fact <-   qnorm(1 - alpha/2)/qnorm(.75)
       B <- apply(B$A, 1:2, quantile, probs = 1:3/4, na.rm = TRUE)
       D <- .5 * fact *(B[3,,]-B[1,,]) * sqmn
       L <- coef - D
       U <- coef + D
       S <- (U - L)/(2 * qnorm(1 - alpha/2))
       T <- coef/S
       P <- 2 * (1 - pnorm(abs(T)))
       G <- list()
       cnames <- c("Value","Lower Bd","Upper Bd","Std Error","T Value","Pr(>|t|)")
       for(i in 1:length(taus)){
	   tab <- cbind(coef[,i],L[,i],U[,i],S[,i],T[,i],P[,i])
	   dimnames(tab)[[2]] <- cnames
	   G[[i]] <- list(tau = taus[i], coefficients = tab, NAs = nas[i])
	   }
       class(G) <- "summary.crqs"
       return(G)
       }
    else
       stop("Invalid method for summary.crq")
   }
"summary.crqs" <-
function (object, ...) {
	if(!object$method == "Powell")
		stop("Invalid method")
        taus <- object$tau
        xsum <- as.list(taus)
        for(i in 1:length(taus)){
                xi <- object
                xi$coefficients <- xi$coefficients[,i]
                xi$residuals <- xi$residuals[,i]
                xi$tau <- xi$tau[i]
                class(xi) <- "crq"
                xsum[[i]] <- summary(xi, ...)
                }
        class(xsum) <- "summary.crqs"
        xsum
        }

print.crq <- function(x, ...){
	if(!is.null(cl <- x$call)) {
                cat("Call:\n")
                dput(cl)
		}
        coef <- coef(x)
        cat("\nCoefficients:\n")
	print.default(coef)
        invisible(x)
	}

print.summary.crqs <- function(x, ...)
    lapply(x,print.summary.crq)
print.summary.crq <- function (x, digits = max(5, .Options$digits - 2), ...) {
    coef <- x$coefficients
    tau <- x$tau
    NAs <- x$NAs
    cat("\ntau: ")
    print(format(round(tau, digits = digits)), quote = FALSE, ...)
    if(NAs > 0) {
       cat("         Number of NA Bootstrap Replications:  ")
       print(format(NAs), quote = FALSE, ...)
       }
    cat("\nCoefficients:\n")
    print(format(round(coef, digits = digits)), quote = FALSE, ...)
    invisible(x)
}