### Functions

library(SemiPar)
library(lokern)
##############################################################
myplot.spm <- function (spm.obj, x.out=NULL, ...) 
{
    object <- spm.obj
	if (is.null(x.out)) x.out <- c(object$info$pen$x)
    plot.params <- plotControl(...)
    drv <- plot.params$drv
    se <- plot.params$se
    shade <- plot.params$shade
    default.ylim <- FALSE
    if (is.null(plot.params$ylim)) 
        default.ylim <- TRUE
    pen.present <- !is.null(object$info$pen)
    random.present <- !is.null(object$info$random)
    const.only <- FALSE    
     block.inds <- object$aux$block.inds

        coefs <- c(object$fit$coef$fixed, object$fit$coef$random)
    if (se == TRUE) 
        cov.mat <- object$aux$cov.mat
    if (random.present) {
        random.inds <- block.inds[[length(block.inds)]]
        coefs <- coefs[-random.inds]
        if (se == TRUE) 
            cov.mat <- cov.mat[-random.inds, -random.inds]
    }
    basis <- object$info$pen$basis
    if (drv == 0) 
        ave.vals <- 1
    if (drv > 0) 
        ave.vals <- 0

        num.pen <- ncol(as.matrix(object$info$pen$x))
        for (ipen in 1:num.pen) {
            deg.val <- object$info$pen$degree[ipen]
            if (basis == "trunc.poly") 
                ncol.X.val <- deg.val
            if (basis == "tps") 
                ncol.X.val <- (deg.val - 1)/2
            ave.val <- mean(as.matrix(object$info$pen$x)[, ipen])
            ave.vals <- c(ave.vals, ave.val)
            if (ncol.X.val > 1) 
                for (ipow in 2:ncol.X.val) ave.vals <- c(ave.vals, 
                  ave.val^ipow)
        }

     if (pen.present) {
        num.pen <- ncol(as.matrix(object$info$pen$x))
        for (ipen in 1:num.pen) {
            deg.val <- object$info$pen$degree[ipen]
            knots <- object$info$pen$knots[[ipen]]
            ave.val <- mean(as.matrix(object$info$pen$x)[, ipen])
            if (basis == "trunc.poly") {
                ave.val.knots <- ave.val - knots
                ave.val.knots <- (ave.val.knots * (ave.val.knots > 
                  0))^deg.val
            }
            if (basis == "tps") {
                ave.val.knots <- abs(ave.val - knots)^deg.val
                sqrt.Omega <- object$info$trans.mat[[ipen]]
                ave.val.knots <- t(solve(sqrt.Omega, ave.val.knots))
            }
            ave.vals <- c(ave.vals, ave.val.knots)
        }
    }

        curve.type <- NULL
        if (pen.present) 
            curve.type <- c(curve.type, rep("pen", ncol(as.matrix(object$info$pen$x))))
        num.curves <- length(curve.type)

    plot.params <- set.plot.dflts(object, plot.params, num.curves)
    if ((!is.null(plot.params$xlim)) & (!is.list(plot.params$xlim))) {
        if (length(plot.params$xlim) != 2) 
            stop("illegal xlim")
        plot.params$xlim <- list(lower = rep(plot.params$xlim[1], 
            num.curves), upper = rep(plot.params$xlim[2], num.curves))
    }

    if ((!is.null(plot.params$ylim)) & (!default.ylim) & (!is.list(plot.params$ylim))) {
        if (length(plot.params$ylim) != 2) 
            stop("illegal ylim")
        plot.params$ylim <- list(lower = rep(plot.params$ylim[1], 
            num.curves), upper = rep(plot.params$ylim[2], num.curves))
    }

    if (pen.present) {
        x.vals <- NULL
        x.vals <- cbind(x.vals, object$info$pen$x)
        fc.stt.pos <- 2
        pen.num <- 1
        for (j in 1:num.curves) {
#            grid.size <- plot.params$grid.size[j]
            grid.size <- length(x.out)
            if (drv == 0) 
                C.grid <- matrix(rep(ave.vals, grid.size), grid.size, 
                  length(ave.vals), byrow = TRUE)
            if (drv > 0) 
                C.grid <- matrix(0, grid.size, length(ave.vals))
#            x.grid <- seq(plot.params$xlim$lower[j], plot.params$xlim$upper[j], 
#                length = plot.params$grid.size[j])
		x.grid <- x.out
            if (curve.type[j] == "pen") 
                deg.val <- object$info$pen$degree[pen.num]
            X.g.inst <- NULL

            if (curve.type[j] == "pen") {
                if (basis == "trunc.poly") 
                  ncol.X.val <- deg.val
                if (basis == "tps") 
                  ncol.X.val <- (deg.val - 1)/2
                if (drv == 0) {
                  for (pow in 1:ncol.X.val) {
                    new.col.data <- x.vals[, j]^pow
                    new.col <- x.grid^pow
                    X.g.inst <- cbind(X.g.inst, new.col)
                  }
                }
                if (drv > 0) {
                  for (pow in 1:ncol.X.val) {
                    new.col.data <- x.vals[, j]^pow
                    pow.drv <- pow - drv
                    if (pow.drv >= 0) 
                      new.col <- prod(pow:(pow.drv + 1)) * (x.grid^pow.drv)
                    else new.col <- rep(0, length(x.grid))
                    X.g.inst <- cbind(X.g.inst, new.col)
                  }
                }
            }
            C.g.inst <- X.g.inst

            if (curve.type[j] == "pen") {
                if (basis == "trunc.poly") 
                  ncol.X.val <- deg.val
                if (basis == "tps") 
                  ncol.X.val <- (deg.val - 1)/2
                inst.col.inds <- fc.stt.pos:(fc.stt.pos + ncol.X.val - 
                  1)
                fc.stt.pos <- fc.stt.pos + ncol.X.val
                knots <- object$info$pen$knots[[pen.num]]
                num.knots <- length(knots)
                Z.g.inst <- outer(x.grid, knots, "-")
                if (basis == "trunc.poly") {
                  if (drv == 0) 
                    Z.g.inst <- (Z.g.inst * (Z.g.inst > 0))^deg.val
                  else {
                    if (deg.val >= drv) {
                      mfac <- prod((deg.val - drv + 1):deg.val)
                      Z.g.inst <- mfac * (Z.g.inst * (Z.g.inst > 
                        0))^(deg.val - drv)
                    }
                    else Z.g.inst <- matrix(0, length(x.grid), 
                      length(knots))
                  }
                }
                if (basis == "tps") {
                  if (drv == 0) {
                    Z.g.inst <- abs(Z.g.inst)^deg.val
                    sqrt.Omega <- object$info$trans.mat[[pen.num]]
                    Z.g.inst <- t(solve(sqrt.Omega, t(Z.g.inst)))
                  }
                  else {
                    if (deg.val >= drv) {
                      mfac <- prod((deg.val - drv + 1):deg.val)
                      Z.g.inst <- mfac * (Z.g.inst^(deg.val - 
                        drv - 1) * abs(Z.g.inst))
                      sqrt.Omega <- object$info$trans.mat[[pen.num]]
                      Z.g.inst <- t(solve(sqrt.Omega, t(Z.g.inst)))
                    }
                    else Z.g.inst <- matrix(0, length(x.grid), 
                      length(knots))
                  }
                }
                C.g.inst <- cbind(C.g.inst, Z.g.inst)
                inst.col.inds <- c(inst.col.inds, block.inds[[j + 
                  1]][-(1:ncol.X.val)])
                pen.num <- pen.num + 1
            }
            C.grid[, inst.col.inds] <- C.g.inst
            y.grid <- as.vector(C.grid %*% coefs)
            if (default.ylim) {
                plot.params$ylim$lower[j] <- min(y.grid)
                plot.params$ylim$upper[j] <- max(y.grid)
            }
            if (se == TRUE) {
                se.grid <- sqrt(diag(C.grid %*% cov.mat %*% t(C.grid)))
                lower <- y.grid - 2 * se.grid
                upper <- y.grid + 2 * se.grid
                if (default.ylim) {
                  plot.params$ylim$lower[j] <- min(lower)
                  plot.params$ylim$upper[j] <- max(upper)
                }
                ylim.val <- range(c(lower, upper))
            }
            if (plot.params$plot.it[j] == TRUE) {
                plot(x.grid, y.grid, type = "n", bty = plot.params$bty[j], 
                  main = plot.params$main[j], xlab = plot.params$xlab[j], 
                  ylab = plot.params$ylab[j], xlim = c(plot.params$xlim$lower[j], 
                    plot.params$xlim$upper[j]), ylim = c(plot.params$ylim$lower[j], 
                    plot.params$ylim$upper[j]))
                if (se == TRUE) {
                  if (shade == FALSE) {
                    lines(x.grid, lower, lty = plot.params$se.lty[j], 
                      lwd = plot.params$se.lwd[j], col = plot.params$se.col[j])
                    lines(x.grid, upper, lty = plot.params$se.lty[j], 
                      lwd = plot.params$se.lwd[j], col = plot.params$se.col[j])
                  }
                  if (shade == TRUE) {
                    mat <- cbind(x.grid, lower, upper)
                    mat <- mat[sort.list(mat[, 1]), ]
                    x.grid <- mat[, 1]
                    lower <- mat[, 2]
                    upper <- mat[, 3]
                    polygon(c(x.grid, rev(x.grid)), c(lower, 
                      rev(upper)), col = plot.params$shade.col[j], 
                      border = FALSE)
                  }
                }
                if ((drv >= 1) & (plot.params$zero.line == TRUE)) 
                  abline(h = 0, err = -1)
                lines(x.grid, y.grid, lty = plot.params$lty[j], 
                  lwd = plot.params$lwd[j], col = plot.params$col[j])
                if (plot.params$jitter.rug == TRUE) 
                  rug.x <- jitter(x.vals[, j])
                if (plot.params$jitter.rug == FALSE) 
                  rug.x <- x.vals[, j]
                rug(rug.x, quiet = 1, col = plot.params$rug.col[j])
            }
        }
    }

    if (se) return(list(x=x.grid,y=y.grid,se=se.grid))
    if(!se) return(list(x=x.grid,y=y.grid))

    invisible()
}
###########################################################

myspm <- function (form, random = NULL, group = NULL, family = "gaussian",
    spar.method = "REML", omit.missing = NULL) 
{
    require("nlme")
    spm.info <- spmFormRead(form, omit.missing)
    random.info <- NULL
    if (!is.null(random)) {
        random.info <- random.read(random, group)
        spm.info <- c(spm.info, list(random = random.info))
        group <- as.numeric(factor(group))
    }
    if (!is.null(unlist(spm.info$pen$spar))) {
        if (any(unlist(spm.info$pen$spar) == 0)) 
            stop("zero smoothing parameters not supported in current version.")
    }
    design.info <- spmDesign(spm.info)
    X <- design.info$X
    Z <- design.info$Z
    Z.spline <- design.info$Z.spline
    y <- spm.info$y
    trans.mat <- design.info$trans.mat
    spm.info <- c(spm.info, list(trans.mat = trans.mat))
    block.inds <- design.info$block.inds
    re.block.inds <- design.info$re.block.inds
    col.ones <- rep(1, nrow(X))
    auto.spar.select <- FALSE
    if (!is.null(spm.info$pen)) {
        auto.spar <- 0
        for (j in 1:length(spm.info$pen$name)) auto.spar <- (auto.spar + 
            (is.null(spm.info$pen$spar[[j]]) & (spm.info$pen$adf[[j]] == 
                "miss")))
        if (auto.spar > 0) 
            auto.spar.select <- TRUE
    }

    if (auto.spar.select == FALSE) {
        if (!is.null(spm.info$lin)) 
            num.lin <- ncol(as.matrix(spm.info$lin$x))
        if (is.null(spm.info$lin)) 
            num.lin <- 0
        compon.num <- 1
        if (!is.null(spm.info$pen)) {
            basis.type <- spm.info$pen$basis
            for (j in 1:ncol(as.matrix(spm.info$pen$x))) {
                if (!is.null(spm.info$pen$adf[[j]]) & (spm.info$pen$adf[[j]] != 
                  "miss")) {
                  deg.val <- spm.info$pen$degree[j]
                  stt.ind <- block.inds[[compon.num + num.lin + 
                    1]][1]
                  ncol.Xj <- length(block.inds[[compon.num + 
                    num.lin + 1]])
                  ncol.Xj <- ncol.Xj - length(re.block.inds[[compon.num]])
                  Xj <- X[, stt.ind:(stt.ind + ncol.Xj - 1)]
                  Xj <- cbind(col.ones, Xj)
                  Zj <- Z[, re.block.inds[[compon.num]]]
                  adf.val <- spm.info$pen$adf[[j]]
                    spar.val <- df.to.spar(adf.val + 1, Xj, Zj)
                  if (basis.type == "trunc.poly")
                    spm.info$pen$spar[[j]] <- exp(log(spar.val)/(2 * 
                      deg.val))
                  else spm.info$pen$spar[[j]] <- exp(log(spar.val)/deg.val)
                  compon.num <- compon.num + 1
                }
            }
        }

        diag.G <- NULL
        if (!is.null(spm.info$pen)) 
            for (j in 1:ncol(as.matrix(spm.info$pen$x))) {
                deg.val <- spm.info$pen$degree[j]
                spar.val <- spm.info$pen$spar[[j]]
                num <- length(spm.info$pen$knots[[j]])
                if (basis.type == "trunc.poly") 
                  diag.G <- c(diag.G, rep(1/(exp((2 * deg.val) * 
                    log(spar.val))), num))
                else diag.G <- c(diag.G, rep(1/(exp((deg.val) * 
                  log(spar.val))), num))
            }
    }

    dummy.group.vec <- col.ones
    fdummy.group.vec <- factor(dummy.group.vec)
    assign("dummy.group.vec.Handan", dummy.group.vec, pos = 1)
    assign("fdummy.group.vec.Handan", fdummy.group.vec, pos = 1)
    assign("group.Handan", group, pos = 1)
    assign("X.Declan", X, pos = 1)
    if (!is.null(Z)) {
        if (is.null(random)) {
            assign("Z.Jaida", Z, pos = 1)
            data.fr <- groupedData(y ~ -1 + X.Declan | dummy.group.vec.Handan, 
                data = data.frame(y, X.Declan, Z.Jaida, dummy.group.vec.Handan))
            Z.block <- list()
            for (i in 1:length(re.block.inds)) Z.block[[i]] <- as.formula(paste("~Z.Jaida[,c(", 
                paste(re.block.inds[[i]], collapse = ","), ")]-1"))
            if (length(re.block.inds) == 1) {
                  lme.fit <- lme(y ~ -1 + X.Declan, random = pdIdent(~-1 +
                    Z.Jaida), data = data.fr, method = spar.method)
            }
            if (length(re.block.inds) > 1) {
                  lme.fit <- lme(y ~ -1 + X.Declan, random = list(dummy.group.vec.Handan = pdBlocked(Z.block,
                    pdClass = rep("pdIdent", length(Z.block)))), 
                    data = data.fr, method = spar.method)
            }
        }
        if (!is.null(random)) {
            assign("Z.Jaida", Z, pos = 1)
            assign("Z.spline.Jaida", Z.spline, pos = 1)
            data.fr <- groupedData(y ~ -1 + X.Declan | dummy.group.vec.Handan, 
                data = data.frame(y, X.Declan, Z.Jaida, dummy.group.vec.Handan))
            Z.block <- list()
            for (i in 1:length(re.block.inds)) Z.block[[i]] <- as.formula(paste("~Z.Jaida[,c(", 
                paste(re.block.inds[[i]], collapse = ","), ")]-1"))
            if (length(re.block.inds) == 1) {
                  stop("implement later; use pigs to test")
            }

            if (length(re.block.inds) > 1) {
                  Z.block <- list(list(fdummy.group.vec.Handan = pdIdent(~Z.spline.Jaida -
                    1)), list(group.Handan = pdIdent(~1)))
                  Z.block <- unlist(Z.block, recursive = FALSE)
                  data.fr <- groupedData(y ~ -1 + X.Declan | 
                    fdummy.group.vec.Handan, data = data.frame(y, 
                    X.Declan, Z.spline.Jaida, group))
                  lme.fit <- lme(y ~ -1 + X.Declan, data = data.fr, 
                    random = Z.block, method = spar.method)
            }
        }
        lme.fit <- c(lme.fit, list(sigma = summary(lme.fit)$sigma))
    }
    if (is.null(Z)) {
        data.fr <- cbind(y, X.Declan, dummy.group.vec.Handan)
        G <- NULL
            lm.fit <- lm(y ~ -1 + X.Declan)
            lme.fit <- list(coef = list(fixed = lm.fit$coef), 
                sigma = summary(lm.fit)$sigma)

    }
    if (!is.null(Z)) {
        lme.fit$coef$random <- unlist(lme.fit$coef$random)
        sig.u.hat <- lme.fit$sigma * exp(unlist(lme.fit$modelStruct))
        diag.sqrt.G <- NULL
        for (ib in 1:length(re.block.inds)) diag.sqrt.G <- c(diag.sqrt.G, 
            rep(sig.u.hat[ib], length(re.block.inds[[ib]])))
        G <- diag(diag.sqrt.G^2)
    }
    resid.var <- lme.fit$sigma^2
    if (auto.spar.select == FALSE) {
            G <- resid.var * diag(diag.G)
            qr.out <- lmeFitQr(y, X, Z, G, resid.var = resid.var)
            coef.ests <- qr.out$coefficients[1:(ncol(X) + ncol(Z))]
            lme.fit <- list()
            lme.fit$coef$fixed <- coef.ests[1:ncol(X)]
            lme.fit$coef$random <- coef.ests[(1 + ncol(X)):length(coef.ests)]
    }

    if (auto.spar.select == TRUE) {
        if ((!is.null(spm.info$pen)) | (!is.null(spm.info$krige))) {
            sigu2.hat <- rep(0, length(re.block.inds))
            if (!is.null(Z)) 
                for (ib in 1:length(re.block.inds)) sigu2.hat[ib] <- diag(G)[re.block.inds[[ib]][1]]
            if (is.null(spm.info$krige)) {
                basis.type <- spm.info$pen$basis
                for (ip in 1:ncol(as.matrix(spm.info$pen$x))) {
                  deg.val <- spm.info$pen$degree[ip]
                  if (basis.type == "trunc.poly") 
                    spm.info$pen$spar[[ip]] <- exp(log(resid.var/sigu2.hat[ip])/(2 * 
                      deg.val))
                  else spm.info$pen$spar[[ip]] <- exp(log(resid.var/sigu2.hat[ip])/deg.val)
                }
            }

            if ((!is.null(spm.info$pen)) & (!is.null(spm.info$krige))) {
                basis.type <- spm.info$pen$basis
                for (ip in 1:ncol(as.matrix(spm.info$pen$x))) {
                  deg.val <- spm.info$pen$degree[ip]
                  var.rats <- (resid.var/sigu2.hat[ip])
                  if (basis.type == "trunc.poly") 
                    spm.info$pen$spar[[ip]] <- var.rats^(1/(2 * 
                      deg.val))
                  else spm.info$pen$spar[[ip]] <- var.rats^(1/(deg.val))
                }
                num.pen <- ncol(as.matrix(spm.info$pen$x))
                var.rat <- (resid.var/sigu2.hat[num.pen + 1])
                deg.val <- spm.info$krige$degree
                spm.info$krige$spar <- exp(log(var.rat)/deg.val)
            }
        }
    }

        aux.info <- lmeAux(X, Z, G, resid.var, block.inds)

    if (!is.null(Z)) {
        coef.ests <- c(lme.fit$coef$fixed, lme.fit$coef$random)
        C.mat <- cbind(X, Z)
    }
    if (is.null(Z)) {
        coef.ests <- lme.fit$coef$fixed
        C.mat <- X
    }
    mins <- NULL
    maxs <- NULL
    for (j in 1:length(block.inds)) {
        fitted.j <- as.matrix(C.mat[, block.inds[[j]]]) %*% as.matrix(coef.ests[block.inds[[j]]])
        mins[j] <- min(fitted.j)
        maxs[j] <- max(fitted.j)
    }

    aux.info <- c(aux.info, list(mins = mins, maxs = maxs))
        fitted <- as.vector(C.mat %*% coef.ests)
        resids <- y - fitted
        lme.fit$fitted <- fitted
        lme.fit$residuals <- resids

    names(lme.fit$coef$fixed) <- NULL
    names(lme.fit$coef$random) <- NULL
    spm.fit.obj <- list(fit = lme.fit, info = spm.info, aux = aux.info)
    class(spm.fit.obj) <- "spm"
    return(spm.fit.obj)
}
###########################################################

features.mat <- function(x, ymat, smoother=c("glkerns", "spm", "smooth.spline"), se=FALSE, npts=max(100, 2*length(x)), 
	plot.it=FALSE, c.out=3, fits.return=FALSE, decim.out=2, ...) 
# Feature extraction
# Author:  Ravi Varadhan, June 23, 2009
#
# Function arguments:
# x = vector of independent variable, e.g. time
# ymat = matrix of time-series to be smoothed
# npts = number of points to use in computing features
# smoother = technique to use for automatic smoothing; only two options are allowed
# plot.it = a logical variable indicating whether the smoothed results should be plotted
# se = a logical variable indicating whether standard error of fitted values should be calculated
# c.out = a constant denoting number of standard deviations away from smooth fit for determining whether a point is an #outlier
#
# Extracted features:
# 	frms = Root-mean-squared function value over the range of data
# 	fmean = Mean function value over the range of data
# 	fsd = square-root of the average variance of the function around its mean
# 	noise = standard deviation of residuals after smoothing
# 	snr = signal-to-noise ratio = fsd/noise
#     fwiggle = root-mean-squared second-derivative of function; a measure of "wiggliness" in the function
# 	deriv.range = minimum and maximum of first derivative of smoothed function scaled by fsd
# 	critical.pts = X-locations where the first derivative is zero
# 	curvature = second derivative of smoothed function at critical points ( > 0 for minimum, and < 0 for maximum)
#       scaled by fsd
# 	outliers = points that are potential outliers
{
trapezoid <- function(x,y) sum(diff(x)*(y[-1] + y[-length(y)]))/2

smoother <- match.arg(smoother, c("glkerns", "spm", "smooth.spline"))
ipkgs <- installed.packages()

if (smoother == "spm") {
	if ("SemiPar" %in% ipkgs[,1]) require(SemiPar, quietly=TRUE) 
	else stop("Install package `SemiPar'", call.=FALSE)
	deriv1 <- function(z) myplot.spm(fit, drv=1, plot.it=FALSE, x.out=z, se=se)$y
	deriv2 <- function(z) myplot.spm(fit, drv=2, plot.it=FALSE, x.out=z, se=se)$y
} 
if (smoother == "glkerns") {
	if ("lokern" %in% ipkgs[,1]) require(lokern, quietly=TRUE) 
	else stop("Install package `lokern'", call.=FALSE)
	deriv1 <- function(z) glkerns(x, y, deriv=1, x.out=z, hetero=TRUE)$est
	deriv2 <- function(z) glkerns(x, y, deriv=2, x.out=z, hetero=TRUE)$est
}

if (smoother == "smooth.spline") {
	deriv1 <- function(z) predict(fit, deriv=1, x=z)$y
	deriv2 <- function(z) predict(fit, deriv=2, x=z)$y
}

x.out <- seq(min(x), max(x), length=npts)
x <<- x

if (is.null(dim(ymat))) ymat <- matrix(ymat, nrow=1, ncol=length(ymat))

nr <- nrow(ymat)
nc <- ncol(ymat)
fmat <- matrix(NA, nr, 10)
cpmat <- matrix(NA, nr, max(ceiling(nc/10), 50))
cvmat <- matrix(NA, nr, max(ceiling(nc/10), 50))
olmat <- matrix(NA, nr, max(ceiling(nc/10), 50))

	if (plot.it & nr > 1) par(ask=TRUE)
	if (fits.return) fits <- fits1 <- fits2 <- matrix(NA, nr, nc)

for (k in 1:nr) {
y <- as.numeric(ymat[k, ])

if (smoother == "spm") {
y <<- y

fit <- try(spm(y ~ f(x, basis="trunc.poly", degree=3), omit.missing=TRUE), silent=TRUE)
if (class(fit) == "try-error") fit <- try(spm(y ~ f(x, basis="trunc.poly", degree=2), omit.missing=TRUE), silent=TRUE)
if (class(fit) == "try-error") stop("Fitting error in spm")

fit0 <- myplot.spm(fit, drv=0, plot.it=FALSE, x.out=x.out, se=se)$y
fit1 <- myplot.spm(fit, drv=1, plot.it=FALSE, x.out=x.out, se=se)$y
fit2 <- myplot.spm(fit, drv=2, plot.it=FALSE, x.out=x.out, se=se)$y
resid <- y - fit$fit$fitted
resid.scaled <- abs(scale(resid))
if (fits.return) {
	fits[k, ] <- fit$fit$fitted
	fits1[k, ] <- myplot.spm(fit, drv=1, plot.it=FALSE, x.out=x, se=se)$y
	fits2[k, ] <- myplot.spm(fit, drv=2, plot.it=FALSE, x.out=x, se=se)$y
	}
}

if (smoother == "glkerns") {
fit <- glkerns(x, y, deriv=0, x.out=x, hetero=FALSE)
fit0 <- glkerns(x, y, deriv=0, x.out=x.out, hetero=TRUE)$est
fit1 <- glkerns(x, y, deriv=1, x.out=x.out, hetero=TRUE)$est
fit2 <- glkerns(x, y, deriv=2, x.out=x.out, hetero=TRUE)$est
resid <- y - fit$est
resid.scaled <- abs(scale(resid))
if (fits.return) {
	fits[k, ] <- fit$est
	fits1[k, ] <- glkerns(x, y, deriv=1, x.out=x, hetero=TRUE)$est
	fits2[k, ] <- glkerns(x, y, deriv=2, x.out=x, hetero=TRUE)$est
	}
}

if (smoother == "smooth.spline") {
fit <- smooth.spline(x, y, cv=FALSE)
fit0 <- predict(fit, deriv=0, x=x.out)$y
fit1 <- predict(fit, deriv=1, x=x.out)$y
fit2 <- predict(fit, deriv=2, x=x.out)$y
resid <- y - predict(fit)$y
resid.scaled <- abs(scale(resid))
if (fits.return) {
	fits[k, ] <- predict(fit)$y
	fits1[k, ] <- predict(fit, deriv=1)$y
	fits2[k, ] <- predict(fit, deriv=2)$y
	}
}

fmean <- trapezoid(x.out, fit0) / diff(range(x.out))
fstar <- sqrt(trapezoid(x.out, fit0^2) / diff(range(x.out)))
fsd <- sqrt(fstar^2 - fmean^2)
noise <- attr(resid.scaled, "scaled:scale")
snr <- fsd/noise
fwiggle <- sqrt(trapezoid(x.out, fit2^2) / diff(range(x.out)))

d0 <- (fit1[-1] * fit1[-npts]) < 0
ncrt <- sum(d0)
if (ncrt == 0) crtpts <- curv <- NULL  else {
	crtpts <- curv <- rep(NA, ncrt)
	ind0 <- (1:npts)[d0]
	for (i in 1:ncrt) {
	temp <- try(uniroot(interval=c(x.out[ind0[i]], x.out[1 + ind0[i]]), f=deriv1), silent=TRUE)
	if (class(temp) != "try-error") {
		crtpts[i] <- temp$root
		curv[i] <- deriv2(temp$root)
		}
	}
    }

rm(fit)
outl <- resid.scaled > c.out
if (sum(outl) == 0) outliers <- NULL else outliers <- x[outl] 

if (plot.it) {
	plot(x, y, type="p", ...)
	lines(x.out, fit0, col=2)
	}


if (!is.null(crtpts) ) {
	cpmat[k, 1:length(crtpts)] <- crtpts
	cvmat[k, 1:length(crtpts)] <- curv
	}

if (!is.null(outliers) ) olmat[k, 1:length(outliers)] <- outliers

fmat[k, ] <- as.numeric(c(fmean, as.numeric(range(fit0)), fsd, noise, snr, as.numeric(range(fit1)), fwiggle, length(crtpts)))
}

fmat <- as.data.frame(fmat)
names(fmat) <- c("fmean", "fmin", "fmax", "fsd", "noise", "snr", "d1min", "d1max", "fwiggle", "ncpts")

ncpmat <- max(apply(cpmat, 1, function(x) sum(!is.na(x))))
ncvmat <- max(apply(cvmat, 1, function(x) sum(!is.na(x))))
nolmat <- max(apply(olmat, 1, function(x) sum(!is.na(x))))
if (ncpmat==0) cpmat <- NULL else cpmat <- round(cpmat[, 1:ncpmat, drop=FALSE], decim.out)
if (ncvmat==0) cvmat <- NULL else cvmat <- round(cvmat[, 1:ncvmat, drop=FALSE], decim.out)
if (nolmat==0) olmat <- NULL else olmat <- olmat[, 1:nolmat, drop=FALSE]

ret.obj <- list(fmat=round(fmat, decim.out), cptmat=cpmat, curvmat=cvmat, olmat=olmat)
if (fits.return) attr(ret.obj, "fits") <- list(fits0=fits, fits1=fits1, fits2=fits2)

return(ret.obj)
}
##################################################