.packageName <- "cobs"
####-*- mode: R; kept-old-versions: 12;  kept-new-versions: 20; -*-

n1000cut <- function(n) floor(ifelse(n > 1000, 671+log(n)^3, n))
## original had  "670 +" but that was *not* monotone

if(!exists("is.R", mode="function"))
    is.R <- function()
    exists("version") && !is.null(vl <- version$language) && vl == "R"

## S+ does not allow "cut(*, labels = FALSE)" -- use cut00() for compatibility:
if(is.R()) {
    .First.lib <- function(lib, pkg) library.dynam("cobs",pkg,lib)
    cut00 <- function(x, breaks)
        cut.default(x, breaks, labels = FALSE, include.lowest = TRUE)
} else { ## S-plus  (tested only with S+ 6.0):
    cut00 <- function(x, breaks)
        as.integer(cut.default(x, breaks, include.lowest = TRUE))
}

cobs <- function(x, y, constraint = c("none", "increase", "decrease",
                       "convex", "concave", "periodic"),
                 knots, nknots, method = "quantile",
                 degree = 2, tau = 0.5, lambda = 0, ic = "aic",
                 n.sub = n1000cut(n),
                 knots.add = FALSE, pointwise = NULL,
                 print.warn = TRUE, print.mesg = TRUE, trace = print.mesg,
                 coef = rep(0,nvar), w = rep(1,n),
                 maxiter = 20*n, lstart = 7872, toler.kn = 1e-6,
                 eps = .Machine$double.eps, factor = 1)
{
    ##=########################################################################
    ##
    ## S interface for He and Ng (1997), ``COBS-Qualitatively Constrained
    ##   Smoothing via Linear Programming''.
    ##
    ##=########################################################################

    ## lstart had default = log(big)^2 , where big = .Machine$single.xmax
    ##        this is 7871.74216945922  for IEEE
    ## preamble
    ##
    cl <- match.call()
    if(!is.R() && !is.loaded(symbol.For("drqssbc")))# keep former S+ setup
        dyn.load("cobs_l.o")
    constraint <- match.arg(constraint)
    na.idx <- is.na(x) | is.na(y)
    x <- x[!na.idx]
    y <- y[!na.idx]
    minx <- min(x)
    maxx <- max(x)
    n <- nrq <- length(x)
    tmin <- 1/n
    ox <- order(x)
    xo <- x[ox]
    nj0 <- 1
    lam <- 1
    Tlambda <- lambda
    if(length(unique(y)) == 2 && print.warn)
        ## warn(7)
        cat("\n It looks like you are fitting a binary choice model.",
            "We recommend pre-smoothing your data using smooth.spline(),",
            "loess() or ksmooth() before calling COBS\n", sep="\n ")

    if(is.null(pointwise)) {
        equal <- greater <- smaller <- gradient <- NULL
        n.equal <- n.greater <- n.smaller <- n.gradient <- 0
    }
    else { ## `pointwise'
        if(!is.matrix(pointwise)|| dim(pointwise)[2] != 3)
            stop("  Argument `pointwise' has to be a three-column matrix.")
        kind <- pointwise[,1] # .Alias
        equal   <- pointwise[kind ==  0, , drop = FALSE]
        greater <- pointwise[kind ==  1, , drop = FALSE]
        smaller <- pointwise[kind == -1, , drop = FALSE]
        gradient<- pointwise[kind ==  2, , drop = FALSE]
        n.equal   <- nrow(equal) # maybe 0 in R and Sv4 (i.e. S+ (>=5)
        n.greater <- nrow(greater)
        n.smaller <- nrow(smaller)
        n.gradient<- nrow(gradient)
    }
    ##
    ## generate default knots sequence
    ##
    Tnknots <- if(lambda == 0) 6 else 20
    if(missing(knots)) {
        mk.flag <- TRUE
        if(missing(nknots))
            nknots <- Tnknots
        if(method == "quantile") {
            lux <- length(ux <- unique(xo))
            if(lux <= nknots) {
                nknots <- lux
                knots <- ux
            } else { # nknots < lux: take ``rounded'' quantiles
                knots  <- ux[seq(1, lux, len = nknots)]
            }
        }
        else ## "equidistant" :
            knots <- seq(xo[1],xo[n], len = nknots)
    }
    else {
        knots <- sort(knots)
        mk.flag <- missing(nknots) || nknots != length(knots)
        names(knots) <- NULL
        nknots <- length(knots)
        if(knots[1] > minx || knots[nknots] < maxx)
            stop("  The range of knots should cover the range of x.")
    }
    if(nknots < 2) stop("  A minimum of two knots is needed.")
    if(nknots == 2) {
        if(lambda == 0 && mk.flag)
            stop("  Can't perform automatic knot selection with nknots == 2.")
        else if(degree == 1)
            stop("  You need at least 3 knots when lambda!=0 and degree==1")
    }
    if(nknots != Tnknots && print.warn)
        cat("\n You are using",nknots,
            "knots instead of the default number of",Tnknots,"knots.\n")
    knots[1] <- knots[1]-toler.kn
    knots[nknots] <- knots[nknots]+toler.kn
    ## make sure that there is at least one observation between any pair of
    ## adjacent knots
    if(length(unique(cut00(x, knots))) != nknots - 1)
        stop(" There is at least one pair of adjacent knots that contains no observation.")
    kmax <- nknots
    ##
    ## set up proper dimension for the pseudo design matrix
    ##
    dim.o <- getdim(degree,nknots,constraint)
    neqc <- dim.o$n.eqc + n.equal + n.gradient
    nl1 <- 0
    ks <- dim.o$ks
    n.iqc <- dim.o$n.iqc
    niqc  <- n.iqc + n.greater + n.smaller
    if(lambda == 0) { ## quantile B-splines without penalty
        pw <- 0
        nvar <- dim.o$nvar
        ##
        ## compute B-spline coefficients for quantile B-spline with stepwise
        ## knots selection, quantile B-spline with fixed knots
        ##
        rr <- qbsks(x,y,w,pw, knots,nknots, degree,Tlambda,constraint, n.sub,
                    equal,smaller,greater,gradient, coef,maxiter,
                    trace, n.equal,n.smaller,n.greater,n.gradient,
                    nrq,nl1, neqc,nj0, tau,lam,tmin,kmax,lstart,
                    ks,mk.flag, knots.add, ic,
                    print.mesg = print.mesg, factor=factor,
                    tol.kn = toler.kn, eps = eps, print.warn = print.warn)
        knots <- rr$knots
        nknots <- rr$nknots
    }
    else { ## lambda !=0 : quantile smoothing B-Splines with penalty
        if(degree == 1) {
            nl1 <-  nknots - 2
            nvar <- dim.o$nvar
            pw <- rep(1, nl1)
        }
        else {
            nl1 <- 1
            nvar <- dim.o$nvar + 1      #one more parameter for sigma
            niqc <- niqc + 2 * (nknots - 1)
            pw <- rep(1, nknots-1)
        }
        if(lambda < 0) {                # lambda is chosen by sic
            lam <- -1
            Tlambda <- 1
            nj0 <- 50*n ## = nsol <<-- determines size of sol[,] !
        }
        ##
        ## compute B-spline coefficients for quantile smoothing B-spline
        ##
        if(lambda < 0 && print.mesg)
            cat("\n Searching for optimal lambda. This may take a while.\n",
                "  While you are waiting, here is something you can consider\n",
                "  to speed up the process:\n",
                "      (a) Use a smaller number of knots;\n",
                "      (b) Increase `factor' to 2 or 3; \n",
                "      (c) Set lambda==0 to exclude the penalty term.\n")# 3

        rr <- drqssbc(x,y,w,pw, knots, degree,Tlambda,constraint, n.sub,
                      equal,smaller,greater,gradient, coef,maxiter,
                      trace, n.equal,n.smaller,n.greater,n.gradient,
                      nrq,nl1, neqc,niqc,nvar,nj0, tau,lam,tmin,kmax,lstart,
                      factor, eps, print.warn)
    }
    if(rr$ifl != 1) { # had problem
        if(rr$ifl < 1)
            warning(" ifl = ",rr$ifl," < 1 -- should NOT happen !!")
        else
            switch(rr$ifl,

        1, # 1
        stop("At least one of the additional pointwise constraints is not feasible.\n Please check your `pointwise' argument and rerun cobs."), # 2
        ## warn(6,maxiter)
        cat("\n WARNING! The algorithm has not converged after",maxiter,
            "iterations.\n",
            "Increase the `maxiter' counter and restart cobs with both\n",
            "`coef' and `knots' set to the values at the last iteration.\n"), # 3
        stop("  The problem is ill-conditioned."), # ifl = 4
        stop(" ifl = 5 -- shold not have happened"), # 5
        warning(" ifl = 6 --- nj0 = nsol = 50*n  is too small!"), # 6
        stop(" ifl = 7 : lambda < 0  *and*  tau outside [0,1]")  # 7
                   )
    }
    nvar <- rr$nvar
    Tcoef <- rr$coef[1:nvar]
    ##
    ## compute the residual
    ##
    y.hat <- .splValue(degree, knots, Tcoef, xo)
    y.hat <- y.hat[order(ox)]# original (unsorted) ordering

    r <- list(call = cl,
              tau = tau, degree = degree, constraint = constraint,
              pointwise = pointwise,
              coef = Tcoef, knots = knots, ifl = rr$ifl, icyc = rr$icyc,
              k = min(rr$k, nknots-2+ks), k0 = rr$k,
              x.ps = rr$pseudo.x,
              resid = y - y.hat, fitted = y.hat,
              SSy = sum((y - mean(y))^2),
              lambda = rr$lambda,
              pp.lambda = if(lambda < 0) rr$pp.lambda,
              sic       = if(lambda < 0) log(rr$sic))
    class(r) <- "cobs"
    r
}## cobs()

print.cobs <- function(x, digits = getOption("digits"), ...) {
    if(!is.numeric(lam <- x$lambda))
        stop("`x' is not a valid \"cobs\" object")
    cat("COBS ", if(lam == 0) "regression" else "smoothing",
        " spline (degree = ", x$degree, ") from call:\n  ", deparse(x$call),
        "\n{tau=",format(x$tau,digits),"}-quantile",
        ";  dimensionality of fit: ",x$k," (",x$k0,")\n", sep="")
    nkn <- length(x$knots)
    cat("knots[1 .. ", nkn,"]: ", sep = "")
    chk <- format(x$knots[if(nkn <= 5) 1:nkn else c(1:4, nkn)], digits=digits)
    if(nkn > 5) chk[4] <- "... "
    cat(chk, sep = ", "); cat("\n")
    if(lam != 0) {
        cat("lambda =", format(lam, digits = digits))
        if((nlam <- length(x$pp.lambda))) {
            cat(", selected via SIC, out of", nlam, "ones.")
        }
        cat("\n")
    }
    invisible(x)
} # print

summary.cobs <- function(object, digits = getOption("digits"), ...) {
    if(!is.numeric(lam <- object$lambda))
        stop("`object' is not a valid \"cobs\" object")
    print(object, digits = digits, ...)# includes knots
    if(!is.null(pw <- object$pointwise)) {
        cat("with",nrow(pw),"pointwise constraints\n")
    }
    cat("coef  :\n"); print(object$coef, digits = digits, ...)
    tau <- object$tau
    if(abs(tau - 0.50) < 1e-6)
        cat("R^2 = ", round(100 * (1 - sum(object$resid^2) / object$SSy), 2),
            "% ;  ", sep="")
    k <- sum((r <- resid(object)) <= 0)
    n <- length(r)
    cat("empirical tau (over all): ",k,"/",n,"  = ", format(k/n,digits),
        " (target tau : ",tau,")\n", sep="")

    ## add more -- maybe finally *return* an object and define
    ## print.summary.cobs <- function(x, ...)

} # summary

residuals.cobs <- function (object, ...) object$resid
fitted.cobs <- function (object, ...) object$fitted

predict.cobs <-
    function(object, z, minz = knots[1], maxz = knots[nknots], nz = 100,
             interval = c("none", "confidence", "simultaneous", "both"),
             level = 0.95, ...)
{
    if(is.null(knots <- object$knots) ||
       is.null(coef  <- object$coef)  ||
       is.null(degree<- object$degree)  ||
       is.null(tau   <- object$tau)) stop("not a valid `cobs' object")

    interval <- match.arg(interval)

    big        <- if(is.R()) 3.4028234663852886e+38 else .Machine$single.xmax
    ##IN
    single.eps <- if(is.R())1.1920928955078125e-07 else .Machine$single.eps

    nknots <- length(knots)
    ord <- as.integer(degree + 1)
    nvar   <- length(coef)

    ##DBG cat("pr..cobs(): (ord, nknots, nvar) =",ord, nknots, nvar, "\n")

    ##
    ## compute fitted value at z
    ##
    ## MM: why should z be *inside* (even strictly) the knots interval? _FIXME_
    if(missing(z)) {
        if(minz >= maxz) stop("minz >= maxz")
        ##NOT YET (for "R CMD check" compatibility):
        ## zo <- seq(minz, maxz, len = nz)
        zo <- seq(max(minz,knots[1]     + single.eps),
                  min(maxz,knots[nknots]- single.eps), len = nz)
    }
    else {
        zo <- sort(z)
        ##IN zo <- zo[zo > knots[1] & zo < knots[nknots]]
        nz <- length(zo)
    }

    fit <- .splValue(degree, knots, coef, zo)

    if(interval != "none") {
        ##
        ## compute confidence bands
        ## both (pointwise and simultaneous : as cheap as only one !
        ##
        z3 <- .splBasis(ord = ord, knots, ncoef = nknots + degree - 1, xo = zo)
        idx <- cbind(rep(1:nz, rep(ord, nz)),
                     c(outer(1:ord, z3$offsets, "+")))
        X <- matrix(0, nz, nvar)
        X[idx] <- z3$design
        if(any(ibig <- abs(X) > big)) { ## MM: no sense here!
            X[ibig] <- X[ibig] / big^0.5
            warning("re-scaling ", sum(ibig), "values of spline basis `X'")
        }

        Tqr <- qr(crossprod(object$x.ps))
        if(Tqr$rank != dim(object$x.ps)[2])
            stop("The pseudo design matrix is singular; this can most likely be solved by using a smaller lambda")# when obtaining qsbs.out
        ## Improved way of computing
        ## xQx = diag(X %*%  solve(Tqr)  %*% t(X)) :
        tX <- t(X)
        xQx <- colSums(qr.coef(Tqr, tX) * tX)

        res <- object$resid
        n <- length(res)

        s <- shat(res, tau, 1 - level, hs = TRUE)
        cn <- sqrt(xQx * tau * (1 - tau))
        sde <- cn * s
        an <- sqrt(qchisq(level, object$k0))   * sde
        bn <- qt((1 + level)/2, n - object$k0) * sde

        cbind(z = zo, fit = fit,
              cb.lo = fit - an, cb.up = fit + an,
              ci.lo = fit - bn, ci.up = fit + bn)
    }# interval
    else
        cbind(z = zo, fit = fit)
} # predict



getdim <- function(degree, nknots,
                   constraint = c("none", "increase", "decrease",
                   "convex", "concave", "periodic"))
{
    ##=########################################################################
    ##
    ## Compute the appropriate dimension for the pseudo design matrix
    ##
    ##=########################################################################
    if(degree == 1)	ks <- 2
    else if(degree == 2)ks <- 3
    else stop("degree has to be either 1 or 2")
    constraint <- match.arg(constraint)
    nvar <- nknots - 2 + ks
    n.eqc <- # the number of EQuality Constraints
        if(constraint == "periodic") 2 else 0
    n.iqc <- # the number of IneQuality Constraints
    if(constraint == "increase" || constraint == "decrease")
        nknots - as.integer(degree == 1)
    else if(constraint == "concave" || constraint == "convex")
        nknots - 1 - as.integer(degree == 1)
    else if(constraint == "periodic" || constraint == "none")
         0
    list(n.iqc = n.iqc, n.eqc = n.eqc, ks = ks, nvar = nvar)
} ## getdim()

### These are (only) used for confidence intervals :

shat <- function(residual, tau, alpha, hs)
{
    ##=########################################################################
    ##
    ## sparsity estimate from empirical quantile function using residuals
    ##
    ##=########################################################################
    residual <- sort(residual)
    n <- length(residual)
    residual <- c(residual[1], residual, residual[n])
    grid <- c(0, seq(0.5/n, 1 - 0.5/n, 1/n), 1)
    hn <- dn(tau, n, hs = hs, alpha)
    ## for small n,  tau +/- hn might be outside [0,1]
    bound <- pmax(0, pmin(1, c(tau - hn, tau + hn)))
    idx <- cut00(bound, grid)
    lambda <- bound * n - (idx - 1) + 0.5
    return(diff(lambda * residual[idx + 1] +
                (1 - lambda) * residual[idx])/(2 * hn))
}

dn <- function(p, n, hs = FALSE, alpha)
{
    ##=########################################################################
    ##
    ## compute window width for sparsity estimator
    ## at quantile p, level = 1-alpha,  n observations
    ## according to
    ##   Hall and Sheather (1988),  <<-  hs=TRUE,   or
    ##   Bofinger (1975),           <<-  hs=FALSE
    ##=########################################################################
    x0 <- qnorm(p)
    f0 <- dnorm(x0)
    if(as.logical(hs))
        n^(-1/3) * qnorm(1 - alpha/2)^(2/3) *
            ((1.5 * f0^2)/(2 * x0^2 + 1))^(1/3)
    else n^-0.2 * ((4.5 * f0^4)/(2 * x0^2 + 1)^2)^ 0.2
}
cobsOld <- function(x, y, constraint = c("none", "increase", "decrease",
                       "convex", "concave", "periodic"),
                 z, minz = knots[1], maxz = knots[nknots], nz = 100,
                 knots, nknots, method = "quantile",
                 degree = 2, tau = 0.5, lambda = 0, ic = "aic",
                 knots.add = FALSE, alpha = 0.1, pointwise,
                 print.warn = TRUE, print.mesg = TRUE, trace = print.mesg,
                 coef = rep(0,nvar), w = rep(1,n),
                 maxiter = 20*n, lstart = log(big)^2, toler = 1e-6,
                 factor = 1)
{
    ##=########################################################################
    ##
    ## S interface for He and Ng (1997), ``COBS-Qualitatively Constrained
    ##   Smoothing via Linear Programming''.
    ##
    ##=########################################################################

mesg <- function(number,...) {
    ##=########################################################################
    ##
    ## S function to print intermediate messages
    ##
    ##=########################################################################
    switch(number,
           cat("\n"),
           cat("\n Now we are fitting cobs() ...\n"), # 2
           cat("\n Searching for optimal lambda. This may take a while.\n",
               "  While you are waiting, here is something you can consider\n",
               "  to speed up the process:\n",
               "      (a) Use a smaller number of knots;\n",
               "      (b) Increase `factor' to 2 or 3; \n",
               "      (c) Set lambda==0 to exclude the penalty term.\n") # 3
           )
}
    ## preamble
    ##
    if(!is.R() && !is.loaded(symbol.For("drqssbc")))# keep former S+ setup
        dyn.load("cobs_l.o")
    constraint <- match.arg(constraint)
    na.idx <- is.na(x) | is.na(y)
    x <- x[!na.idx]
    y <- y[!na.idx]
    big        <- if(is.R()) 3.4028234663852886e+38 else .Machine$single.xmax
    single.eps <- if(is.R()) 1.1920928955078125e-07 else .Machine$single.eps
    ## FIXME: `single.eps' is only used for cutting off log(<very small>)
    ## double.eps <- .Machine$double.eps
    minx <- min(x)
    maxx <- max(x)
    n <- nrq <- length(x)
    tmin <- 1/n
    ox <- order(x)
    xo <- x[ox]
    yo <- y[ox]
    nj0 <- 1
    lam <- 1
    Tlambda <- lambda
    if(length(unique(y)) == 2 && print.warn)
        ## warn(7)
        cat("\n It looks like you are fitting a binary choice model.",
            "We recommend pre-smoothing your data using smooth.spline(),",
            "loess() or ksmooth() before calling COBS\n", sep="\n ")

    if(missing(pointwise)) {
        equal <- greater <- smaller <- gradient <- NULL
        n.equal <- n.greater <- n.smaller <- n.gradient <- 0
    }
    else { ## `pointwise'
        if(!is.matrix(pointwise)|| dim(pointwise)[2] != 3)
            stop("  Argument `pointwise' has to be a three-column matrix.")
        kind <- pointwise[,1] # .Alias
        equal   <- pointwise[kind ==  0, , drop = FALSE]
        greater <- pointwise[kind ==  1, , drop = FALSE]
        smaller <- pointwise[kind == -1, , drop = FALSE]
        gradient<- pointwise[kind ==  2, , drop = FALSE]
        n.equal   <- nrow(equal) # maybe 0 in R!
        n.greater <- nrow(greater)
        n.smaller <- nrow(smaller)
        n.gradient<- nrow(gradient)
    }
    ##
    ## generate default knots sequence
    ##
    Tnknots <- if(lambda == 0) 6 else 20
    if(missing(knots)) {
        mk.flag <- TRUE
        if(missing(nknots))
            nknots <- Tnknots
        if(method == "quantile") {
            lux <- length(ux <- unique(xo))
            if(lux < nknots) {
                nknots <- lux
            }
            knots <- ux[seq(1, lux, len = nknots)]
            names(knots) <- NULL
        }
        else
            knots <- seq(xo[1],xo[n], len = nknots)
    }
    else {
        knots <- sort(knots)
        mk.flag <- missing(nknots) || nknots != length(knots)
        names(knots) <- NULL
        nknots <- length(knots)
        if(knots[1] > minx || knots[nknots] < maxx)
            stop("  The range of knots should cover the range of x.")
    }
    if(nknots < 2) stop("  A minimum of two knots is needed.")
    if(nknots == 2) {
        if(lambda == 0 && mk.flag)
            stop("  Can't perform automatic knot selection with nknots == 2.")
        else if(degree == 1)
            stop("  You need at least 3 knots when lambda!=0 and degree==1")
    }
    if(nknots != Tnknots && print.warn)
        cat("\n WARNING! It looks like you are using",nknots,
            "knots instead of the\n   default number of",Tnknots,"knots.\n")
    knots[1] <- knots[1]-toler
    knots[nknots] <- knots[nknots]+toler
    ## make sure that there is at least one observation between any pair of
    ## adjacent knots
    if(length(unique(cut(x, knots, labels=FALSE, include.lowest = TRUE))) != nknots - 1)
        stop(" There is at least one pair of adjacent knots that contains no observation.")
    kmax <- nknots
    ##
    ## set up proper dimension for the pseudo design matrix
    ##
    dim.o <- getdim(degree,nknots,constraint)
    neqc <- dim.o$n.eqc + n.equal + n.gradient
    nl1 <- 0
    ks <- dim.o$ks
    n.iqc <- dim.o$n.iqc
    if(lambda == 0) { ## quantile B-splines without penalty
        pw <- 0
        niqc <- n.iqc + n.greater + n.smaller
        nvar <- dim.o$nvar

        ##
        ## compute B-spline coefficients for quantile B-spline with stepwise
        ## knots selection, quantile B-spline with fixed knots
        ##
        qsbs.o <- qbsks(x,y,w,pw, knots,nknots, degree,Tlambda,constraint,
                        n.sub = n1000cut(n),
                        equal,smaller,greater,gradient, coef,maxiter,
                        trace, n.equal,n.smaller,n.greater,n.gradient,
                        nrq,nl1, neqc,nj0, tau,lam,tmin,kmax,lstart,
                        ks,mk.flag, knots.add, ic, print.mesg,## method,
                        factor = factor, print.warn = print.warn)
        knots <- qsbs.o$knots
        nknots <- qsbs.o$nknots
    }
    else { ## lambda !=0 : quantile smoothing B-Splines with penalty
        if(degree == 1) {
            nl1 <-  nknots - 2
            nvar <- dim.o$nvar
            niqc <- n.iqc + n.greater + n.smaller
            pw <- rep(1, nl1)
        }
        else {
            nl1 <- 1
            nvar <- dim.o$nvar + 1      #one more parameter for sigma
            niqc <- 2 * (nknots - 1) + n.iqc + n.greater + n.smaller
            pw <- rep(1, nknots-1)
        }
        if(lambda < 0) {                # lambda is chosen by sic
            lam <- -1
            Tlambda <- 1
            nj0 <- 50*n ## = nsol <<-- determines size of sol[,] !
        }
        ##
        ## compute B-spline coefficients for quantile smoothing B-spline
        ##
        if(lambda < 0 && print.mesg) mesg(3)
        qsbs.o <- drqssbc(x,y,w,pw, knots, degree,Tlambda,constraint,
                          n.sub = n1000cut(n),
                          equal,smaller,greater,gradient, coef,maxiter,
                          trace, n.equal,n.smaller,n.greater,n.gradient,
                          nrq,nl1, neqc,niqc,nvar,nj0, tau,lam,tmin,kmax,lstart,
                          factor, print.warn = print.warn)
    }
    if(qsbs.o$ifl != 1) { # had problem
        if(qsbs.o$ifl < 1)
            warning(" ifl = ",qsbs.o$ifl," < 1 -- should NOT happen !!")
        else
            switch(qsbs.o$ifl,

        1, # 1
        stop("  At least one of the additional pointwise constraints is not feasible. Please check your `pointwise' argument and rerun cobs."), # 2
        ## warn(6,maxiter)
        cat("\n WARNING! The algorithm has not converged after",maxiter,
            "iterations.\n",
            "Increase the `maxiter' counter and restart cobs with both\n",
            "`coef' and `knots' set to the values at the last iteration.\n"), # 3
        stop("  The problem is ill-conditioned."), # ifl = 4
        stop(" ifl = 5 -- shold not have happened"), # 5
        warning(" ifl = 6 --- nj0 = nsol = 50*n  is too small!"), # 6
        stop(" ifl = 7 : lambda < 0  *and*  tau outside [0,1]")  # 7
                   )
    }
    nvar <- qsbs.o$nvar
    ##
    ## compute fitted value at z
    ##
    if(missing(z)) {
        zo <- seq(max(minz,knots[1]     + single.eps),
                  min(maxz,knots[nknots]- single.eps), len = nz)
    }
    else {
        zo <- z[order(z)]
        zo <- zo[zo > knots[1] & zo < knots[nknots]]
        nz <- length(zo)
    }
    oz <- order(zo)

    new.knots <- c(rep(knots[1], ks-1), knots, rep(knots[nknots], ks-1))
    nk <- length(new.knots)
    derivs <- as.integer(0)
    Tcoef <- qsbs.o$coef[1:nvar]
    Tncoef <- length(Tcoef) ## ==!== nvar
    fit <- .C("spline_value",
              as.double(new.knots),
              as.double(Tcoef),
              Tncoef,
              as.integer(ks),
              as.double(zo),
              as.integer(nz),
              derivs,
              y = double(nz), PACKAGE = "cobs") $ y
    ##
    ## compute the residual
    ##
    z2 <- .C("spline_value",
             as.double(new.knots),
             as.double(Tcoef),
             Tncoef,
             as.integer(ks),
             as.double(xo),
             as.integer(n),
             derivs,
             y = double(n), PACKAGE = "cobs")
    resid <- (y[ox] - z2$y)[order(ox)]
    ##
    ## compute the confidence band
    ##
    X <- matrix(0, nz, nvar)
    z3 <- .C("spline_basis",
             as.double(new.knots),
             ncoef = as.integer(nk - ks),
             as.integer(ks),
             as.double(zo),
             derivs = integer(nz),# 0
             as.integer(nz),
             design = array(0, c(ks, nz)),
             offsets = integer(nz), PACKAGE = "cobs")
    idx <- cbind(rep(oz, rep(ks, nz)),
                 c(outer(1:ks, z3$offsets, "+")))
    X[idx] <- z3$design
    if(any(ibig <- abs(X) > big)) { ## make sure that drqssbc won't overflow
        X[ibig] <- X[ibig] / big^0.5
        warning("re-scaling ", sum(ibig), "values of spline basis `X'")
    }
    chisq.alpha <- qchisq(1 - alpha, qsbs.o$k)
    z.alpha <- qt(1 - alpha/2, n - qsbs.o$k)
    ## MM: not clear if this is good: QR( X' X ) ; later inverse
    Tqr <- qr(t(qsbs.o$pseudo.x) %*% qsbs.o$pseudo.x)
    if(Tqr$rank != dim(qsbs.o$pseudo.x)[2])
        stop("The pseudo design matrix is singular; this can most likely be solved by using a smaller lambda when obtaining qsbs.out"
             )
    Q <- solve(Tqr)                     #note: no n here
    xQx <- diag(X %*% Q %*% t(X))

##Dbg cat("after ` Q <- solve(Tqr) ' and `xQx <- diag(X %*% Q %*% t(X))'")
##Dbg browser()

    s <- shat(resid, tau, alpha, hs = TRUE)
    cn <- sqrt(xQx * tau * (1 - tau))
    an <- sqrt(chisq.alpha) * cn * s
    bn <- z.alpha * cn * s

    list(coef = Tcoef, fit = fit, resid = resid,
         z = zo, knots = knots, ifl = qsbs.o$ifl, icyc = qsbs.o$icyc,
         k = min(qsbs.o$k, nknots-2+ks), lambda = qsbs.o$lambda,
         pp.lambda = if(lambda < 0) qsbs.o$pp.lambda,
         sic       = if(lambda < 0) log(qsbs.o$sic),
         cb.lo = fit - an, cb.up = fit + an,
         ci.lo = fit - bn, ci.up = fit + bn)
}## cobsOld()
### used to be part of ./cobs.R

drqssbc <- function(x,y, w = rep(1,n), pw, knots, degree,Tlambda, constraint,
                    n.sub = n1000cut(nrq),
		    equal,smaller, greater,gradient, coef, maxiter = 20*n,
		    trace = 1,
                    n.equal = nrow(equal), n.smaller = nrow(smaller),
                    n.greater = nrow(greater), n.gradient = nrow(gradient),
		    nrq = length(x), nl1, neqc,niqc, nvar,nj0,
		    tau = 0.50, lam, tmin, kmax, lstart, factor = 1,
                    eps = .Machine$double.eps, print.warn = TRUE)
{
    ##=########################################################################
    ##
    ## Estimate the B-spline coefficients for quantile *smoothing* spline, using
    ##		Ng (1996)  `An Algorithm for Quantile Smoothing Splines',
    ##		Computational Statistics & Data Analysis, 22, 99-118.
    ##
    ##=########################################################################
    big	       <- if(is.R()) 3.4028234663852886e+38 else .Machine$single.xmax
    single.eps <- if(is.R()) 1.1920928955078125e-07 else .Machine$single.eps
    toler.kn <- 1e-6
    ## Note   nrq != length(x) , e.g., in case of sub.sampling+ fit full
    if(lam >= 0)
        n <- nrq
    else if((n.old <- nrq) != (n <- as.integer(n.sub))) {
	##
	## sub-sampling for smoothing B-splines parametric programming
        ## select a sub-sample of size n.sub
	##
	sub.idx <- seq(1,n.old, length = n)
	x.old <- x; x <- x[sub.idx]
	y.old <- y; y <- y[sub.idx]
	w.old <- w; w <- w[sub.idx]
    }
    if(degree == 1) {
	X <- l1.design(x,w,constraint,equal,smaller,greater,gradient,knots,
		       pw,n.equal,n.smaller,n.greater,n.gradient,
		       nrq=n, nl1,neqc,niqc,nvar,Tlambda)
	niqc1 <- 0
    }
    else {
	X <- loo.design(x,w,constraint,equal,smaller,greater,gradient,knots,
			pw,n.equal,n.smaller,n.greater,n.gradient,
			nrq=n, nl1,neqc,niqc,nvar,Tlambda)
	niqc1 <- if(Tlambda == 0) 0 else 2*(length(knots) - 1)
    }
    if(any(ibig <- abs(X) > big)) { ## make sure that drqssbc won't overflow
	X[ibig] <- X[ibig] / big^0.5
	warning("re-scaling ", sum(ibig), "values of Lp-design `X'")
    }
    Tnobs <- nrow(X)
    Tequal   <- if(n.equal    > 0)    equal[,3] # else NULL
    Tsmaller <- if(n.smaller  > 0) -smaller[,3]
    Tgreater <- if(n.greater  > 0)  greater[,3]
    Tgradient<- if(n.gradient > 0) gradient[,3]

    Y <- c(y*w, rep(0,nl1), Tequal, Tgradient,
	   rep(0,Tnobs-n-nl1-n.equal-n.gradient-n.smaller-n.greater),
	   Tsmaller,Tgreater)
    ##storage.mode(X) <- "single" # would round to ~ 7 digits in S+, not in R
    d <- matrix(0.0, Tnobs + 5, nvar + 2) # double
    sol <- matrix(0.0, nvar + 6, nj0) # double -- to contain "sol"ution
    z0 <- .Fortran("drqssbc",
		   as.integer(n),
		   as.integer(nl1),
		   as.integer(neqc),
		   as.integer(niqc),
		   as.integer(niqc1),
		   as.integer(nvar),
		   integer(1),		# nact
		   ifl = integer(1),
		   as.integer(maxiter), # mxs
		   as.integer(trace),
		   X = as.double(t(X)), # e
		   as.integer(nvar),	# ner
		   coef = as.double(coef),# x
		   as.double(Y),	# f
		   obj = double(1),	# erql1n
		   resid = double(Tnobs),
		   integer(Tnobs),	# indx
		   double(((3 * nvar + 13) * nvar + 2)/2 + 2 * Tnobs),# w
		   nt = integer(1),	# nt
		   as.integer(nj0),	# nsol
		   sol = sol,
		   as.double(c(tau, lam)),# == tl[1:2] == (t, lam)
		   as.double(toler.kn),
		   as.double(big),
		   as.double(eps),
		   icyc = integer(2),
		   as.double(tmin),
		   k = integer(1),
		   as.integer(kmax),	# k0
		   as.double(lstart),
		   as.double(factor),
                   PACKAGE = "cobs")
    sol <- z0$sol[,1:z0$nt]
    names(z0$icyc) <- c("icyc", "tot.cyc")
    if(lam < 0) {
        ##
        ## search for optimal lambda
        ##
        ifl.idx <-
            if(maxiter > 20*n) sol[3,] != 3 # trap ifl=3
            else rep(TRUE, ncol(sol))
        fidel <- sol[4,ifl.idx]
        eff.k <- sol[6,ifl.idx]
        sic <- fidel/n * n^(eff.k/(2*n))
        mlam.idx <- min((1:ncol(sol[,ifl.idx]))[sic == min(sic)])
        ##
        ## trap infeasible solution when lam<0
        ##
        if(sol[,ifl.idx][3,mlam.idx] == 0)
            return(list(ifl = 2))

        Tcoef <- sol[,ifl.idx][7:nrow(sol),mlam.idx]
        if(print.warn && sol[3,mlam.idx] != 3) {
            if(sol[,ifl.idx][6,mlam.idx] >= length(knots))
                ## warn(2)
                cat("\n WARNING! Since the optimal lambda chosen by SIC corresponds to the",
                    "   roughest possible fit, you might want to consider doing one of",
                    "   the following:",
                    "   (1) plot the components $sic against $pp.lambda of cobs to see",
                    "   if a bigger lambda value at another local minimum of $sic will",
                    "   yield a more reasonable fit;",
                    "   (2) increase the number of knots.\n", sep="\n ")
            else if(abs(sic[mlam.idx]-sic[length(sic)]) < single.eps * sic[length(sic)])
                ## warn(3)
                cat("\n WARNING! Since the optimal lambda chosen by SIC corresponds to the",
                    "   roughest possible fit, you might want to plot the components",
                    "   $sic against $pp.lambda of cobs to see if a bigger lambda value",
                    "   at another local minimum of $sic will yield a more reasonable fit.\n",
                    sep="\n ")

            if(sol[,ifl.idx][2,mlam.idx] == lstart)
                ## warn(4)
                cat("\n WARNING!  Since the optimal lambda chosen by SIC reached the smoothest",
                    "   possible fit allowed by `lstart', you might want to rerun cobs with",
                    "   a larger `lstart' value to see if it makes a difference if you haven't",
                    "   done so.\n", sep="\n ")
        }
        if(n.old != n) { ## was sub sampling; refit full sample for one Tlambda
            Tlambda <- sol[,ifl.idx][2,mlam.idx]
	    rqss <- drqssbc(x.old,y.old,w.old, pw,knots,degree, Tlambda,
			    constraint, n.sub,
			    equal,smaller,greater,gradient,
			    Tcoef,maxiter,trace,
			    n.equal,n.smaller,n.greater,n.gradient,
			    nrq = n.old, nl1,neqc,niqc,nvar, nj0 = 1,
			    tau,lam = 1, tmin,kmax,lstart,
			    factor, eps, print.warn)
	    list(coef = rqss$coef, ifl = sol[,ifl.idx][3,mlam.idx],
		 icyc = rqss$icyc, nvar = nvar, lambda = Tlambda,
		 pp.lambda = sol[,ifl.idx][2,], sic = sic,
		 k = min(rqss$k, length(knots)-2+degree+1), pseudo.x = X)
	}
	else
	    list(coef = Tcoef, ifl = sol[,ifl.idx][3,mlam.idx],
		 icyc = z0$icyc, nvar = nvar,
		 lambda	   = sol[,ifl.idx][2,mlam.idx],
		 pp.lambda = sol[,ifl.idx][2,], sic = sic,
		 k = min(sol[,ifl.idx][6,mlam.idx], length(knots)-2+degree+1),
		 pseudo.x = X)
    }
    else # `lam >= 0'
	list(coef = z0$coef,
	     fidel = sum((tau-(-z0$resid[1:n] < 0))*(-z0$resid[1:n]))*2,
	     k = min(z0$k, length(knots)-2+degree+1), ifl = z0$ifl,
	     icyc = z0$icyc, nvar = nvar, lambda = Tlambda, pseudo.x = X)
}## drqssbc
####  Create B-Spline Design matrices for COBS :
####
####  l1.design()  -- L_1         penalty
####  loo.design() -- L_Infinity  penalty

### used to be part of ./cobs.R

l1.design <-
function(x, w, constraint, equal, smaller, greater, gradient,
         knots, pw, n.equal, n.smaller, n.greater, n.gradient,
         nrq, nl1, neqc, niqc, nvar, lambda)
{
    ##=########################################################################
    ##
    ## Generate the pseudo design matrix for L1 penalty
    ##
    ##=########################################################################

    ## create the pseudo design
    ##
    ##x <- as.matrix(x, ncol = 1)

    ks <- 2 ## <==> degree == 1
    nk <- length(knots) + 2 # *(ks - 1)
    nkm3 <- nk - 3
    nkm4 <- nk - 4
    ncoef <- nk - ks
    ox <- order(x)
    sortx <- x[ox]
    nobs <- nrq + nl1 + neqc + niqc
    X <- matrix(0, nobs, nvar)

    z1 <- .splBasis(ord = ks, knots, ncoef, xo = sortx)
    idx1 <- cbind(rep(ox, rep(ks, nrq)),
                  c(outer(1:ks, z1$offsets, "+")))
    X[idx1] <- t(t(z1$design) * rep(w[ox], ks))
    ##
    ## formulate inequality constraints for the pseudo X
    ##
    if(lambda != 0 || constraint != "none") {
        z2 <- .splBasis(ord = ks, knots, ncoef, xo = knots[1:nkm3],
                        derivs = rep(1, nkm3))
        if(lambda != 0) {
            idx2 <- array(rep(1:nkm4, 2), c(nkm4, 2))
            idx2[, 1] <- idx2[, 1] + nrq
            X[idx2]   <-  - z2$design[1, 1:nkm4] * lambda
            idx2[, 2] <- idx2[, 2] + 2
            X[idx2]   <- z2$design[2, 2:nkm3] * lambda
            idx2[, 2] <- idx2[, 2] - 1
            X[idx2] <- (z2$design[1, 2:nkm3] -
                        z2$design[2, 1:nkm4]) * lambda
            ##
            ## assign different weight to roughness penalty
            ##
            X[nrq + 1:nkm4, ] <- X[nrq + 1:nkm4, ] * pw
        }
        if(constraint == "increase" || constraint == "decrease") {
            niqc1 <- nkm3
            idx3 <- cbind(rep(nrq+nl1+neqc + 1:niqc1, rep(ks, niqc1)),
                          c(outer(1:ks, z2$offsets, "+")))
            X[idx3] <- if(constraint == "increase") z2$design else -z2$design
        }
        else if (constraint == "periodic") {
            ##
            ## this portion corresponds to equality constraints
            ##
            niqc1 <- 0
            neqc3 <- 2
            z1.3 <- .splBasis(ord = ks, knots, ncoef,
                              xo = sortx[c(1,nrq)], derivs = c(1,1))
            idx1.3 <- cbind(rep(1:neqc3, rep(ks, neqc3)),
                            c(outer(1:ks, z1.3$offsets,"+")))
            X.temp <- matrix(0,neqc3,nvar)
            X.temp[idx1.3] <- z1.3$design
            X[nrq+nl1+neqc,  ] <- X.temp[2,] - X.temp[1,]
            X[nrq+nl1+neqc-1,] <- X[ox[nrq],]- X[ox[1],]
        }
        else if(constraint == "convex" || constraint == "concave") {
            niqc1 <- nkm4
            idx3 <- array(rep(1:nkm4, 2), c(nkm4, 2))
            idx3[, 1] <- idx3[, 1] + nrq + nl1 + neqc
            sgn <- if(constraint == "convex") +1 else -1
            X[idx3] <- -sgn* z2$design[1, 1:nkm4]
            idx3[,2] <- idx3[,2] + 2
            X[idx3] <- sgn * z2$design[2, 2:nkm3]
            idx3[,2] <- idx3[,2] - 1
            X[idx3] <- sgn * (z2$design[1, 2:nkm3] - z2$design[2, 1:nkm4])
        }
    }
    niqc2 <- n.smaller
    niqc3 <- n.greater
    if(n.smaller > 0) {
        o.smaller <- order(smaller[,2])
        smaller.o <- smaller[,2][o.smaller]
        z3.1 <- .splBasis(ord = ks, knots, ncoef, xo = smaller.o)
        idx3.1 <- cbind(rep(nrq+nl1+neqc+niqc1 + o.smaller,rep(ks,niqc2)),
                        c(outer(1:ks,z3.1$offsets,"+")))
        X[idx3.1] <- -z3.1$design
    }
    if(n.greater > 0) {
        o.greater <- order(greater[,2])
        greater.o <- greater[,2][o.greater]
        z3.2 <- .splBasis(ord = ks, knots, ncoef, xo = greater.o)
        idx3.2 <- cbind(rep(nrq+nl1+neqc+niqc1+niqc2 + o.greater,rep(ks,niqc3)),
                        c(outer(1:ks,z3.2$offsets,"+")))
        X[idx3.2] <- z3.2$design
    }
    neqc1 <- n.equal
    neqc2 <- n.gradient

    if(n.equal > 0) { ## formulate equality constraints for the pseudo X

        o.equal <- order(equal[,2])
        equal.o <- equal[,2][o.equal]
        z1.1 <- .splBasis(ord = ks, knots, ncoef, xo = equal.o)
        idx1.1 <- cbind(rep(nrq+nl1+ o.equal, rep(ks, neqc1)),
                        c(outer(1:ks, z1.1$offsets,"+")))
        X[idx1.1] <- z1.1$design
    }

    if(n.gradient > 0) { ## gradient constraints for the pseudo X

        o.gradient <- order(gradient[,2])
        gradient.o <- gradient[,2][o.gradient]
        z1.2 <- .splBasis(ord = ks, knots, ncoef, xo = gradient.o,
                          derivs = rep(1, neqc2))
        idx1.2 <- cbind(rep(nrq+nl1+neqc1+ o.gradient, rep(ks, neqc2)),
                        c(outer(1:ks, z1.2$offsets,"+")))
        X[idx1.2] <- z1.2$design
    }
    return(X)
} ## l1.design

loo.design <- function(x, w, constraint, equal, smaller, greater, gradient,
                       knots, pw, n.equal, n.smaller, n.greater, n.gradient,
                       nrq, nl1, neqc, niqc, nvar, lambda)
{
    ##=########################################################################
    ##
    ## Generate the pseudo design matrix for L_oo penalty
    ##
    ##=########################################################################

    ##x <- as.matrix(x, ncol = 1)

    ks <- 3 ## <==> degree == 2
    nk <- length(knots) + 2*(ks - 1)# = length(new.knots)
    ncoef <- nk - ks
    nd <- nk - 5
    ox <- order(x)
    sortx <- x[ox]
    nrql1 <- nrq + nl1
    nrleq <- nrql1 + neqc
    nobs <- nrleq + niqc
    X <- matrix(0, nobs, nvar)
    z1 <- .splBasis(ord = ks, knots, ncoef, xo = sortx)
    idx1 <- cbind(rep(ox, rep(ks, nrq)),
                  c(outer(1:ks, z1$offsets, "+")))
    X[idx1] <- t(t(z1$design) * rep(w[ox], ks))
    ##
    ## formulate the inequality constraints -- s''()
    ##
    if(lambda != 0) {
        niqc1 <- 2*nd
        z2 <- .splBasis(ord = ks, knots, ncoef, xo = knots[1:nd],
                        derivs = rep(2, nd))
        X[(nrq+1):nrql1,] <- cbind(c(rep(0,nvar-1), lambda))
        idx2 <- cbind(rep(nrleq + 1:nd, rep(ks,nd)),
                      c(outer(1:ks,z2$offsets,"+")))
        X[idx2] <- z2$design
        idx2[,1] <- idx2[,1]+nd
        X[idx2] <- -z2$design
        ## assign different weight to roughness penalty
        X[nrleq + 1:niqc1, ] <-
            X[nrleq + 1:niqc1, ] * rep(pw, 2)
        X[nrleq + 1:niqc1, nvar ] <- rep(1,niqc1)
    } else niqc1 <- 0
    niqc2 <- n.smaller
    niqc3 <- n.greater
    niqc4 <- niqc - niqc1 -niqc2 -niqc3
    if(constraint != "none") {
        if(constraint == "convex" || constraint == "concave") {
            if(lambda == 0) # z2 not yet above
                z2 <- .splBasis(ord = ks, knots, ncoef, xo = knots[1:niqc4],
                                derivs = rep(2, niqc4))
            idx3 <- cbind(rep(nrleq+niqc1 + 1:niqc4, rep(ks,niqc4)),
                          c(outer(1:ks,z2$offsets,"+")))
            X[idx3] <- if(constraint == "convex") z2$design else -z2$design
        }
        else if(constraint == "periodic") {
            neqc3 <- 2
            z2 <- .splBasis(ord = ks, knots, ncoef,
                            xo = sortx[c(1,nrq)], derivs = c(1,1))
            idx2 <- cbind(rep(1:neqc3, rep(ks,neqc3)),
                          c(outer(1:ks,z2$offsets,"+")))
            X.temp <- matrix(0,neqc3,nvar)
            X.temp[idx2] <- z2$design
            X[nrleq,  ] <- X.temp[2,]-X.temp[1,]
            X[nrleq-1,] <- X[ox[nrq],]-X[ox[1],]
        }
        else {
            z3 <- .splBasis(ord = ks, knots, ncoef, xo = knots[1:niqc4],
                            derivs = rep(1, niqc4))
            idx3 <- cbind(rep(nrleq+niqc1 + 1:niqc4, rep(ks,niqc4)),
                          c(outer(1:ks,z3$offsets,"+")))
            X[idx3] <- if(constraint == "increase") z3$design else -z3$design
        }
    }
    if(n.smaller > 0) {
        o.smaller <- order(smaller[,2])
        smaller.o <- smaller[,2][o.smaller]
        z3.1 <- .splBasis(ord = ks, knots, ncoef, xo = smaller.o)
        idx3.1 <- cbind(rep(nrleq+niqc1+niqc4 + o.smaller,rep(ks,niqc2)),
                        c(outer(1:ks,z3.1$offsets,"+")))
        X[idx3.1] <- -z3.1$design
    }
    if(n.greater > 0) {
        o.greater <- order(greater[,2])
        greater.o <- greater[,2][o.greater]
        z3.2 <- .splBasis(ord = ks, knots, ncoef, xo = greater.o)
        idx3.2 <- cbind(rep(nrleq+niqc1+niqc4+niqc2 + o.greater, rep(ks,niqc3)),
                        c(outer(1:ks,z3.2$offsets,"+")))
        X[idx3.2] <- z3.2$design
    }
    ##
    ## formulate the equality constraints
    ##
    neqc1 <- n.equal
    neqc2 <- n.gradient
    if(n.equal > 0) {
        derivs <- rep(0, neqc1)
        o.equal <- order(equal[,2])
        equal.o <- equal[,2][o.equal]
        z1.1 <- .splBasis(ord = ks, knots, ncoef, xo = equal.o)
        idx1.1 <- cbind(rep(nrql1 + o.equal, rep(ks, neqc1)),
                        c(outer(1:ks, z1.1$offsets,"+")))
        X[idx1.1] <- z1.1$design
    }
    if(n.gradient > 0) {
        o.gradient <- order(gradient[,2])
        gradient.o <- gradient[,2][o.gradient]
        z1.1 <- .splBasis(ord = ks, knots, ncoef, xo = gradient.o,
                          derivs = rep(1, neqc2))
        idx1.2 <- cbind(rep(nrql1+neqc1 + o.gradient, rep(ks, neqc2)),
                        c(outer(1:ks, z1.2$offsets,"+")))
        X[idx1.2] <- z1.2$design
    }
    return(X)
} ## loo.design
### used to be part of ./cobs.R

qbsks <- function(x,y,w,pw, knots,nknots, degree,Tlambda, constraint,
                  n.sub = n1000cut(n),
                  equal,smaller, greater,gradient, coef,maxiter,
                  trace, n.equal,n.smaller,n.greater,n.gradient,
                  nrq,nl1, neqc, nj0, tau,lam,tmin,kmax,lstart,
                  ks,mk.flag, knots.add, ic, print.mesg,
                  factor, tol.kn = 1e-6, eps = .Machine$double.eps, print.warn)
{
    ##=########################################################################
    ##
    ## Compute B-spline coefficients for quantile B-spline with stepwise knots
    ## selection, quantile B-spline with fixed knots (REGRESSION SPLINE), using
    ##		Ng (1996)  `An Algorithm for Quantile Smoothing Splines',
    ## 		Computational Statistics & Data Analysis, 22, 99-118.
    ##
    ##=########################################################################

    ## single.eps <- if(is.R()) 1.1920928955078125e-07 else .Machine$single.eps
    ## double.eps <- .Machine$double.eps
    smll.log <- 50*floor(log(.Machine$double.xmin)/50) # heavily rounding down
    finite.log <- function(x) {
        r <- log(x)
        if(is.na(r) || r > -Inf) r else smll.log # = -750 for IEEE arithmetic
    }
    n <- n.old <- nrq # "n <-" : for default of n.sub
    n <- nrq <- n.sub <- as.integer(n.sub)
    if(n != n.old) {
        ##
        ## select a sub-sample of size n.sub
        ##
        sub.idx <- seq(1,n.old, length = n)
	x.old <- x; x <- x[sub.idx]
	y.old <- y; y <- y[sub.idx]
	w.old <- w; w <- w[sub.idx]
    }
    xo <- x[order(x)]
    logn <- log(n)
    const <- switch(ic, aic = 2/n, sic = logn/n)
    constraint.old <- constraint
    if(mk.flag) { ##  perform first step knots selection
        if(print.mesg) cat("\n Performing general knot selection ...\n")# 4
        Tic <- Tifl <- double(nknots-1)
        for(i in 1:(nknots-1)) {
            Tknots <- knots[seq(1,nknots, len = i+1)]
            Tnknots <- length(Tknots)
            if(Tnknots == 2 && degree == 1 &&
               (constraint == "convex" || constraint == "concave"))
                ## guard against convex fit when only 2 knots are used
                constraint <- "none"
            dim.o <- getdim(degree, Tnknots, constraint)
            ks <- dim.o$ks
            n.iqc <- dim.o$n.iqc
            Tnvar <- dim.o$nvar
            niqc <- n.iqc + n.greater + n.smaller
            rqss <- drqssbc(x,y,w,pw,Tknots,degree,Tlambda,constraint, n.sub,
                            equal,smaller,greater,gradient,coef,maxiter,
                            trace,n.equal,n.smaller,n.greater,n.gradient,
                            nrq,nl1,neqc,niqc,Tnvar,nj0,
                            tau,lam,tmin,kmax,lstart, factor,eps, print.warn)
            constraint <- constraint.old
            Tic[i] <- finite.log(rqss$fidel) -logn + (i-1+ks)*const
            Tifl[i] <- rqss$ifl
        }
        Tic.min <- min(Tic)
        Tifl.final <- Tifl[Tic == Tic.min]
        nknots.min <- min((1:(nknots-1))[Tic == Tic.min])
        if(nknots.min == 1 && Tifl.final != 1) {
            ##
            ## when the chosen nknots=2, guard against anomaly of ifl=5 when
            ## constraint=='periodic', or ifl=2 when the chosen model is infeasible.
            ##
            Tic.min <- min(Tic[2:length(Tic)])
            Tifl.final <- Tifl[Tic == Tic.min]
            nknots.min <- min((1:(nknots-1))[Tic == Tic.min])
        }
        if(Tifl.final == 2 || Tifl.final == 3 || Tifl.final == 4)
            return(list(ifl = Tifl.final))
        if((nknots.min + 1) == nknots && print.warn)
            ## warn(5,nknots,ic)
            cat("\n WARNING! Since the number of ",nknots," knots selected by ",
                ic," reached the\n",
                "  upper bound during general knot selection, you might want to rerun\n",
                "  cobs with a larger number of knots. \n")

        knots <- knots[seq(1,nknots, len = nknots.min+1)]
        names(knots) <- NULL
        ##
        ## perform knots deletion
        ##
        delete <- TRUE
        if(print.mesg) cat("\n Deleting unnecessary knots ...\n")# 5
        while(delete && nknots.min > 1) {
            Tnknots <- length(knots)
            Tic1 <- rep(0,(Tnknots-2))
            Tnknots.1 <- Tnknots - 1
            if(Tnknots.1 == 2 && degree == 1  &&
               (constraint == "convex" || constraint == "concave"))
                ## guard against convex fit when only 2 knots are used
                constraint <- "none"
            dim.o <- getdim(degree,Tnknots.1,constraint)
            ks <- dim.o$ks
            n.iqc <- dim.o$n.iqc
            Tnvar <- dim.o$nvar
            niqc <- n.iqc + n.greater + n.smaller
            Tcoef <- rep(0,Tnvar)
            for(i in 2:(Tnknots-1)) {
                Tknots <- knots[-i]
                rqss <- drqssbc(x,y,w,pw,Tknots, degree,Tlambda,constraint, n.sub,
                                equal,smaller,greater,gradient,
                                Tcoef, maxiter, trace,
                                n.equal,n.smaller,n.greater,n.gradient,
                                nrq,nl1,neqc,niqc,Tnvar,nj0,
                                tau,lam,tmin,kmax,lstart,factor, eps,print.warn)
                constraint <- constraint.old
                Tic1[i-1] <- finite.log(rqss$fidel)-logn+(Tnknots.1-2+ks)*const
                Tcoef <- rqss$coef
            }
            Tic1.min <- min(Tic1)
            idx.del <- min((2:(Tnknots-1))[Tic1 == Tic1.min])
            if((delete <- Tic1.min <= Tic.min)) {
                Tic.min <- Tic1.min
                if(print.mesg >= 3)
                    cat("\n A knot at ",signif(knots[idx.del]),
                        " is deleted.\n")# 6
                knots <- knots[-idx.del]
                nknots.min <- length(knots)-1
            }
        }
        if(print.mesg >= 2) cat("\n No more knot to be deleted.\n") # 7
        ##
        ## perform knots addition
        ##
        if(knots.add) {
            add <- TRUE
            Tnknots <- length(knots)
            if(print.mesg) cat("\n Searching for missing knots ...\n")# 8
            while(add && Tnknots < nknots) {
                Tic2 <- double(Tnknots-1)
                knots.add <- (knots[1:(Tnknots-1)]+knots[2:Tnknots])/2
                Tnknots.1 <- Tnknots + 1
                dim.o <- getdim(degree,Tnknots.1,constraint)
                ks <- dim.o$ks
                n.iqc <- dim.o$n.iqc
                Tnvar <- dim.o$nvar
                niqc <- n.iqc + n.greater + n.smaller
                Tcoef <- double(Tnvar)
                for(i in 1:(Tnknots-1)) {
                    Tknots <- sort(c(knots,knots.add[i]))
                    if(length(unique(cut00(x, Tknots))) != Tnknots)
                        Tic2[i] <- Tic.min+1
                    else {
                        rqss <-
                            drqssbc(x,y,w,pw,Tknots,degree,Tlambda, constraint, n.sub,
                                    equal,smaller,greater,gradient,
                                    Tcoef,maxiter, trace,
                                    n.equal,n.smaller,n.greater,n.gradient,
                                    nrq,nl1,neqc,niqc,Tnvar,nj0,
                                    tau,lam,tmin,kmax,lstart, factor, eps,print.warn)
                        Tic2[i] <- finite.log(rqss$fidel) -logn +
                            (Tnknots.1-2+ks)*const
                        Tcoef <- rqss$coef
                    }
                }
                Tic2.min <- min(Tic2)
                idx.add <- min((1:(Tnknots-1))[Tic2 == Tic2.min])
                if((add <- Tic2.min <= Tic.min)) {
                    Tic.min <- Tic2.min
                    knots <- sort(c(knots,knots.add[idx.add]))
                    if(print.mesg >= 2)
                        cat("\n A knot at ",signif(knots.add[idx.add]),
                            " is added.\n")# 9
                }
                Tnknots <- length(knots)
            }# end while(add ..)
            if(print.mesg >= 2) cat("\n No more knot to be added.\n")# 10
        } # (knots.add)
        if(print.mesg) cat("\n Computing the final fit ...\n")# 11
    } # end if(mk.flag)

    ##
    ## compute the B-spline coefficients for the full sample
    ##
    nknots <- length(knots)
    rk <- diff(range(knots))
    knots[1] <- knots[1] - tol.kn*rk
    knots[nknots] <- knots[nknots] + tol.kn*rk
    if(n != n.old) {
        x <- x.old
        y <- y.old
        w <- w.old
        nrq <- n.old
    }
    if(nknots == 2 && (constraint == "convex" || constraint == "concave") &&
       degree == 1) # guard against convex fit when only 2 knots are used
        constraint <- "none"
    dim.o <- getdim(degree,nknots,constraint)
    ks <- dim.o$ks
    n.iqc <- dim.o$n.iqc
    Tnvar <- dim.o$nvar
    niqc <- n.iqc + n.greater + n.smaller
    rqss <- drqssbc(x,y,w,pw, knots, degree,Tlambda,constraint, n.sub,
                    equal,smaller,greater,gradient, coef,maxiter,
                    trace, n.equal,n.smaller,n.greater,n.gradient,
                    nrq,nl1, neqc,niqc, Tnvar,nj0,
                    tau,lam,tmin,kmax,lstart,
                    factor, eps,print.warn)
    constraint <- constraint.old
    list(coef = rqss$coef, fidel = rqss$fidel,
         k = nknots-2+ks, ifl = rqss$ifl, icyc = rqss$icyc,
         knots = knots, nknots = nknots, nvar = Tnvar, lambda = Tlambda,
         pseudo.x = rqss$pseudo.x)
}# end qbsks()
#### Abstract out the .C() calls into these auxiliaries;
#### Makes it easier later to replace by calls to  (*, PACKAGE = "splines")

###--> ../tests/spline-ex.R  for investigating these
###    ~~~~~~~~~~~~~~~~~~~~

.splBasis <- function(ord, knots, ncoef, xo, derivs = rep(0, n))
{
    ## Purpose: encapsulate .C("spline_basis", ..)
    ## ----------------------------------------------------------------------
    ## Arguments: from result of B-spline fit, see ?cobs , etc
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 20 Feb 2002, 14:46

    ord <- as.integer(ord)
    new.knots <- c(rep(knots[1], ord-1),
                   knots,
                   rep(knots[length(knots)], ord-1))

    if(ord + length(knots) != ncoef + 2)
        warning(".splBasis(): (ord,length(knots),ncoef)=",
                paste(ord,length(knots),ncoef, sep=", "),
                " -- not ``matching'' ?\n")

    n <- length(xo <- as.double(xo))

    .C("spline_basis",
       as.double(new.knots),
       as.integer(ncoef),
       ord, # "order"
       xo, # "xvals"
       derivs = as.integer(derivs),
       n,
       design = array(0, c(ord, n)),# "basis"
       offsets = integer(n),
       PACKAGE = "cobs")[c("design","offsets")]
}

.splValue <- function(degree, knots, coef, xo)
{
    ## Purpose: encapsulate .C("spline_value", ..)
    ## ----------------------------------------------------------------------
    ## Arguments: from result of B-spline fit, see ?cobs , etc
    ## ----------------------------------------------------------------------
    ## Author: Martin Maechler, Date: 20 Feb 2002, 13:48

    degree <- as.integer(degree)
    ord <- as.integer(degree + 1)
    new.knots <- c(rep(knots[1], degree),
                   knots,
                   rep(knots[length(knots)], degree))
    derivs <- as.integer(0)
    n <- length(xo)
    .C("spline_value",
       as.double(new.knots),
       as.double(coef),
       length(coef),
       ord,
       as.double(xo),
       as.integer(n),
       derivs,
       y = double(n), PACKAGE = "cobs")$y
} ## .splValue
