ksmooth <-
  function(x, y, kernel=c("box", "normal"), bandwidth=0.5, range.x=range(x),
           n.points=max(100, length(x)), x.points)
{
# box is [-0.5, 0.5]. normal is sd = 1.4826/4
  if(missing(y))
    stop("y must be supplied.\nFor density estimation use density()")
  kernel <- match.arg(kernel)
  krn <- switch(kernel, "box" = 1, "normal" = 2)
  if(missing(x.points))
    x.points <- seq(range.x[1], range.x[2], len=n.points)
  else {
    n.points <- length(x.points)
    x.points <- sort(x.points)
  }
  ord <- order(x)
  z <- .C("BDRksmooth",
          as.double(x[ord]),
          as.double(y[ord]),
          as.integer(length(x)),
          xp=as.double(x.points),
          yp=double(n.points),
          as.integer(n.points),
          as.integer(krn),
          as.double(bandwidth)
          )
  list(x=z$xp, y=z$yp)
} 

ppr <- function(x, ...) UseMethod("ppr")

ppr.formula <-
function(formula, data=sys.parent(), weights, subset,
	 na.action, contrasts=NULL, ...)
{
  call <- match.call()
  m <- match.call(expand = F)
  m$contrasts <- m$... <- NULL
  m[[1]] <- as.name("model.frame")
  m <- eval(m, sys.parent())
  Terms <- attr(m, "terms")
  attr(Terms, "intercept") <- 0
  X <- model.matrix(Terms, m, contrasts)
  Y <- model.extract(m, response)
  w <- model.extract(m, weights)
  if(length(w) == 0) w <- rep(1, nrow(X))
  fit <- ppr.default(X, Y, w, ...)
  fit$terms <- Terms
  fit$call <- call
  structure(fit, class=c("ppr.form", "ppr"))
}

ppr.default <- 
function(x, y, weights=rep(1,n), ww=rep(1,q), nterms, max.terms=nterms, 
         optlevel=2, sm.method=c("supsmu", "spline", "gcvspline"),
         bass=0, span=0, df=5, gcvpen=1)
{
  call <- match.call()
  sm.method <- match.arg(sm.method)
  ism <- switch(sm.method, supsmu=0, spline=1, gcvspline=2)
  if(missing(nterms)) stop("nterms is missing with no default")
  mu <- nterms; ml <- max.terms
  x <- as.matrix(x)
  y <- as.matrix(y)
  n <- nrow(x)
  if(nrow(y) != n) stop("mismatched x and y")
  p <- ncol(x)
  q <- ncol(y)
  if(!is.null(dimnames(x))) xnames <- dimnames(x)[[2]]
  else xnames <- paste("X", 1:p, sep="")
  if(!is.null(dimnames(y))) ynames <- dimnames(y)[[2]]
  else ynames <- paste("Y", 1:p, sep="")  
  msmod <- ml*(p+q+2*n)+q+7+ml		# for asr
  nsp <- n*(q+15)+q+3*p
  ndp <- p*(p+1)/2+6*p
  .Fortran("bdrsetppr", 
	   as.double(span), as.double(bass), as.integer(optlevel),
	   as.integer(ism), as.double(df), as.double(gcvpen)
	   )
  Z <- .Fortran("bdrsmart",
		as.integer(ml), as.integer(mu),
		as.integer(p), as.integer(q), as.integer(n),
		as.double(weights),
		as.double(t(x)),
		as.double(t(y)),
		as.double(ww),
		smod=double(msmod), as.integer(msmod),
		double(nsp), as.double(nsp),
		double(ndp), as.double(ndp),
		edf=double(ml)
		)
  smod <- Z$smod
  ys <- smod[q+6]
  tnames <- paste("term", 1:mu)
  alpha <- matrix(smod[q+6 + 1:(p*mu)],p, mu,
		  dimnames=list(xnames, tnames))
  beta <- matrix(smod[q+6+p*ml + 1:(q*mu)], q, mu,
		     dimnames=list(ynames, tnames))
  fitted <- drop(matrix(.Fortran("bdrpred",
		       as.integer(nrow(x)),
		       as.double(x),
		       as.double(smod),
		       y = double(nrow(x)*q),
		       double(2*smod[4])
		       )$y, ncol=q, dimnames=dimnames(y)))
  jt <- q + 7 + ml*(p+q+2*n)
  gof <- smod[jt] * n * ys^2
  gofn <- smod[jt+1:ml] * n * ys^2
# retain only terms for the size of model finally fitted
  jf <- q+6+ml*(p+q)
  smod <- smod[c(1:(q+6+p*mu), q+6+p*ml + 1:(q*mu), 
		 jf + 1:(mu*n), jf+ml*n + 1:(mu*n))]
  smod[1] <- mu
  structure(list(call=call, ml=max.terms, p=p, q=q, 
		 gof=gof, gofn=gofn, 
		 df=df, edf=Z$edf[1:mu],
		 xnames=xnames, ynames=ynames,
		 alpha=drop(alpha), beta=ys*drop(beta),
		 yb=smod[5+1:q], ys=ys,
		 fitted.values=fitted, residuals=drop(y-fitted),
		 smod=smod), 
	    class="ppr")
}

print.ppr <- function(x, ...)
{
  if(!is.null(cl <- x$call)) {
    cat("Call:\n")
    dput(cl)
  }
  mu <- x$call$nterms; ml <- x$ml
  cat("\nGoodness of fit:\n")
  gof <- x$gofn; names(gof) <- paste(1:ml, "terms")
  print(format(gof[mu:ml], ...), quote=F)
  invisible(x)
}

summary.ppr <- function(x, ...)
{
  print.ppr(x, ...)
  mu <- x$call$nterms
  cat("\nProjection direction vectors:\n")
  print(format(x$alpha, ...), quote=F)
  cat("\nCoefficients of ridge terms:\n")
  print(format(x$beta, ...), quote=F)
  if(any(x$edf >0)) {
    cat("\nEquivalent df for ridge terms:\n")
    edf <- x$edf; names(edf) <- paste("term", 1:mu)
    print(round(edf,2), ...)
  }
  invisible(x)
}

plot.ppr <- function(fit, ask, type="o", ...)
{
  ppr.funs <- function(obj)
    {
	# cols for each term
      p <- obj$p; q <- obj$q
      sm <- obj$smod
      n <- sm[4]; mu <- sm[5]; m <- sm[1]
      jf <- q+6+m*(p+q)
      jt <- jf+m*n
      f <- matrix(sm[jf+1:(mu*n)],n, mu)
      t <- matrix(sm[jt+1:(mu*n)],n, mu)
      list(x=t, y=f)
    }
  obj <- ppr.funs(fit)
  if(!missing(ask)) {
    oldpar <- par()
    on.exit(par(oldpar))
    par(ask = ask)
  }
  for(i in 1:fit$call$nterms) {
    ord <- order(obj$x[ ,i])
    plot(obj$x[ord, i], obj$y[ord, i], type = type,
	 xlab = paste("term", i), ylab = "", ...)
  }
}

predict.ppr <- function(obj, newdata)
{
  if(missing(newdata)) return(obj$fitted)
  if(!is.null(obj$terms))
    x <- model.matrix(delete.response(obj$terms), newdata)
  else x <- as.matrix(newdata)
  if(ncol(x) != obj$p) stop("wrong number of columns in x")
  drop(matrix(.Fortran("bdrpred",
		       as.integer(nrow(x)),
		       as.double(x),
		       as.double(obj$smod),
		       y = double(nrow(x)*obj$q),
		       double(2*obj$smod[4])
		       )$y,
	      ncol=obj$q,
	      dimnames=list(dimnames(x)[[1]], obj$ynames)
	      ))
}
smooth.spline <-
  function(x, y, w = rep(1, length(x)), df = 5, spar = 0, cv = FALSE,
           all.knots = FALSE, df.offset = 0, penalty = 1)
{
  sknotl <- function(x)
  {
    a1 <- log(50, 2)
    a2 <- log(100, 2)
    a3 <- log(140, 2)
    a4 <- log(200, 2)
    n <- length(x)
    if(n < 50) nk <- n
    else if (n >= 50  && n < 200) nk <- 2^(a1+(a2-a1)*(n-50)/150) 
    else if (n >= 200 && n < 800) nk <- 2^(a2+(a3-a2)*(n-200)/600)
    else if (n >= 800 && n < 3200)nk <-  2^(a3+(a4-a3)*(n-800)/2400)
    else if (n >= 3200) nk <- 200 + (n-3200)^0.2
    c(rep(x[1], 3), x[seq(1,n,len=trunc(nk))], rep(x[n], 3))
  }
  if(missing(y)) {
    if(is.list(x)) {
      if(any(is.na(match(c("x", "y"), names(x)))))
        stop("cannot find x and y in list")
      y <- x$y
      x <- x$x
    } else if(is.complex(x)) {
      y <- Im(x)
      x <- Re(x)
    } else if(is.matrix(x) && ncol(x) == 2) {
      y <- x[, 2]
      x <- x[, 1]
    } else {
      y <- x
      x <- time(x)
    }
  }
  n <- length(x)
  if(n != length(y)) stop("lengths of x and y must match")
  if(missing(w)) w <- rep(1, n)
  else {
    if(length(x) != length(w)) stop("lengths of x and w must match")
    if(any(w < 0)) stop("all weights should be non-negative")
    w <- (w * sum(w > 0))/sum(w)	#this makes the weights sum to n
  }
  if(missing(spar)) ispar <- 0
  else if(spar < 1.01e-15) ispar <- 0 else ispar <- 1
  if(cv)  icrit <- 2 else  icrit <- 1
  dfinfo <- df.offset
  if(!missing(df)) {
    if(df > 1 & df < n) {
      icrit <- 3
      dfinfo <- df
    } else warning("you must supply 1 < df < n")
  }
  n <- as.integer(n)
  x <- signif(x, 6)
  ux <- unique(sort(x))
  ox <- match(x, ux)
  tmp <- matrix(unlist(tapply(seq(along=y), ox,
                              function(x,y,w) c(mean(y[x]), sum(w[x])),
                              y = y, w = w)), , 2, byrow=T)
  ybar <- tmp[,1]
  wbar <- tmp[,2 ]
  nx <- length(ux)
  xbar <- (ux - ux[1])/(ux[nx] - ux[1])
  if(all.knots) {
    knot <- c(rep(xbar[1], 3), xbar, rep(xbar[nx], 3))
    nk <- nx + 2
  } else {
    knot <- sknotl(xbar)
    nk <- length(knot) - 4
  }
  low.parm <- 0
  high.parm <- 1.5
  fit <- .Fortran("qsbart",
                  as.double(penalty),
                  as.double(dfinfo),
                  x = as.double(xbar),
                  y = as.double(ybar),
                  w = as.double(wbar),
                  as.integer(nx),
                  as.double(knot),
                  as.integer(nk),
                  coef = double(nk),
                  ty = double(nx),
                  lev = double(nx),
                  crit = double(1),
                  as.integer(c(icrit, ispar)),
                  spar = as.double(spar),
                  as.double(c(0, 1.5, 0.001)),
                  as.integer(0),
                  double((17 + nk) * nk),
                  as.integer(4),
                  as.integer(1),
                  ier = as.integer(1))
  if(fit$ier > 0) {
    warning("smoothing parameter value too small or too large")
    fit$ty <- rep(mean(y), nx)
  }
  lev <- fit$lev
  df <- sum(lev)
  if(cv) {
    ww <- wbar
    ww[!(ww > 0)] <- 1
    cv.crit <- weighted.mean(((y - fit$ty[ox])/(1 - (lev[ox] * w)/ww[ox]))^2,
                             w)
  } else cv.crit <- weighted.mean((y - fit$ty[ox])^2, w)/
    (1 - (df.offset + penalty * df)/sum(wbar))^2
  pen.crit <- sum(wbar * (ybar - fit$ty) * ybar)
  fit.object <- list(knot = knot, nk = nk, min = ux[1],
                     range = ux[nx] - ux[1], coef = fit$coef)
  class(fit.object) <- "smooth.spline.fit"
  object <- list(x = ux, y = fit$ty, w = wbar, yin = ybar, 
                 lev = lev, cv.crit = cv.crit, pen.crit = pen.crit, df = df, 
                 spar = fit$spar, fit = fit.object, call = match.call())
  class(object) <- "smooth.spline"
  object
}

predict.smooth.spline <- function(object, x, deriv = 0)
{
  if(missing(x)) return(object[c("x", "y")])
  fit <- object$fit
  if(is.null(fit)) stop("not a valid smooth.spline object")
  else predict(fit, x, deriv)
}

print.smooth.spline <- function(x, ...)
{
  if(!is.null(cl <- x$call)) {
    cat("Call:\n")
    dput(cl)
  }
  if(is.null(cv <- cl$cv)) cv <- F
  cat("\nSmoothing Parameter (Spar):", format(x$spar), "\n")
  cat("Equivalent Degrees of Freedom (Df):", format(x$df), "\n")
  cat("Penalized Criterion:", format(x$pen.crit), "\n")
  crss <- if(cv) "PRESS:" else "GCV:"
  cat(crss, format(x$cv.crit), "\n")
  invisible(x)
}

predict.smooth.spline.fit <- function(object, x, deriv = 0)
{
  if(missing(x))
    x <- seq(from = object$min, to = object$min + object$range, 
             length = length(object$coef) - 4)
  xs <- (x - object$min)/object$range
  extrap.left <- xs < 0
  extrap.right <- xs > 1
  interp <- !(extrap.left | extrap.right)
  y <- xs
  if(any(interp))
     y[interp] <- .Fortran("bvalus",
                           as.integer(sum(interp)),
                           as.double(object$knot),
                           as.double(object$coef),
                           as.integer(object$nk),
                           as.double(xs[interp]),
                           s = double(sum(interp)),
                           as.integer(deriv))$s
  if(any(!interp)) {
    xrange <- c(object$min, object$min + object$range)
    if(deriv == 0) {
      end.object <- Recall(object, xrange)$y
      end.slopes <- Recall(object, xrange, 1)$y * object$range
      if(any(extrap.left))
        y[extrap.left] <- end.object[1] + end.slopes[1] * (xs[extrap.left] - 0)
      if(any(extrap.right))
        y[extrap.right] <- end.object[2] + end.slopes[2] * (xs[extrap.right] - 1)
    } else if(deriv == 1) {
      end.slopes <- Recall(object, xrange, 1)$y * object$range
      y[extrap.left] <- end.slopes[1]
      y[extrap.right] <- end.slopes[2]
    }
    else y[!interp] <- 0
  }
  if(deriv > 0) y <- y/(object$range^deriv)
  list(x = x, y = y)
}

supsmu <-
  function(x, y, wt = rep(1, length(y)), span = "cv", periodic = F, bass = 0)
{
  if(span == "cv") span <- 0
  n <- length(y)
  if(length(x) != n)
    stop("number of observations in x and y must match.")
  if(length(wt) != n)
    stop("number of weights must match number of observations.")
  if(span < 0 || span > 1) stop("span must be between 0 and 1.")
  if(periodic) {
    iper <- 2
    xrange <- range(x)
    if(xrange[1] < 0 || xrange[2] > 1)
      stop("x must be between 0 and 1 for periodic smooth")
  } else iper <- 1
  okay <- is.finite(x + y + wt)
  ord <- order(x[okay], y[okay])
  ord <- cumsum(!okay)[okay][ord] + ord
  xo <- x[ord]
  leno <- length(ord)
  if(diff <- n - leno)
    warning(paste(diff, "observation(s) with NAs, NaNs and/or Infs deleted"))
  .Fortran("bdrsetsmu")
  smo <- .Fortran("bdrsupsmu",
                       as.integer(leno),
                       as.double(xo),
                       as.double(y[ord]),
                       as.double(wt[ord]),
                       as.integer(iper),
                       as.double(span),
                       as.double(bass),
                       smo=double(leno),
                       double(n*7), double(1))$smo
  #eliminate duplicate xsort values and corresponding smoothed values
  dupx <- duplicated(xo)
  list(x = xo[!dupx], y = smo[!dupx])
}
.First.lib <- function(lib, pkg) library.dynam("ppr", pkg, lib)
