.packageName <- "klaR"
EDAM <- function(EV0, nzx = 0, iter.max = 10, random = TRUE, standardize = FALSE, wghts = 0, classes = 0,
    sa = TRUE, temp.in = 0.5, temp.fin = 0.0000001, temp.gamma = 0){
    if (is.data.frame(EV0)) EV0 <- as.matrix(EV0)
    diss <- FALSE
    nrEV0 <- nrow(EV0)
    if (ncol(EV0) == nrEV0 && EV0 == t(EV0)){
        diss <- TRUE
        EV.dist <- EV0
        EV0 <- cbind(1:nrEV0, 1:nrEV0)
        rownames(EV0) <- rownames(EV.dist)
    }
    EV.keep <- EV0    
    if (standardize){
        EV.sd <- sqrt(apply(EV0, 2, var))
        EV.sd[!EV.sd] <- 1
        EV0 <- t(t(EV0) / EV.sd)
    } 
    if(wghts[1]) EV0 <- EV0 * kronecker(t(wghts), rep(1, nrEV0))
    EV0.name <- deparse(substitute(EV0))
    if(!nzx){
        flsqrtc <- floor(sqrt(nrEV0))
        nzx <- max(which(nrEV0 / (1:flsqrtc) == floor(nrEV0 / (1:flsqrtc))))
    }
    if (nzx > nrEV0){
        warning("Given argument nzx (", nzx, ") is bigger than rownumber (", nrEV0, 
        ") of argument ", EV0.name, ". So nzx has been set to ", nrEV0, ".",  call. = FALSE)
        nzx <- nrEV0
    }
    fnzy <- floor(nrEV0/nzx)
    if (fnzy != nrEV0/nzx){
        warning("Rownumber (", nrEV0, ") of argument ", EV0.name, " is not a multiple of given argument nzx (",
            nzx,"). ", "\n", "A ", nzx, "x", floor(nrEV0/nzx),"-Map from the first ", nzx * fnzy, 
            " rows of ", EV0.name, " was constructed instead.", call. = FALSE)
            EV0 <- EV0[1:(nzx * fnzy),]
            nrEV0 <- nrow(EV0)
            if (diss) EV.dist <- EV.dist[1:(nzx*fnzy), 1:(nzx*fnzy)]
    }
    nzy <- nrEV0/nzx
    Z0 <- matrix(1:(nzx*nzy),nzy,nzx,byrow=TRUE)
    Cells0 <- cbind(kronecker(1:nzy, rep(1,nzx)), rep(1:nzx, nzy))
    if (!temp.gamma) temp.gamma <- (temp.fin/temp.in)^(1/(max(nzx,nzy)-3))


    TopoS <- function(EV.dist, Cells.dist){
        dim(EV.dist) <- NULL
        beta.est <- (1/(Cells.dist %*% Cells.dist)) %*% Cells.dist %*% EV.dist
        Cells.dist.est <- (Cells.dist * beta.est) - EV.dist
        return(1 - sqrt((Cells.dist.est %*% Cells.dist.est) / (EV.dist %*% EV.dist)))
    }

    sim.ann <- function(EV.old, EV, EV.dist.old, EV.dist, S.old, S.new, temperature){                           
        change.prob <- exp(-(S.old-S.new)/temperature)
        if (runif(1) > change.prob){
            EV <- EV.old 
            EV.dist <- EV.dist.old}               
        return(list(EV=EV, EV.dist=EV.dist))
    }
    Cells0.dist <- distmirr(dist(Cells0,method = "euclidean"))
    dim(Cells0.dist) <- NULL
    EV <- EV0
    EV.best <- EV
    if(random) EV <- EV0[sample(1:(nzy*nzx)),]
    if(!diss) EV.dist <- distmirr(dist(EV, method = "euclidean"))    
    EV.dist.best <- EV.dist
    S.initial <- TopoS(EV.dist,Cells0.dist)
    S.memo <- rep(0,iter.max*(max(nzy,nzx)-1))
    S.best <- S.initial
    ng <- max(nzy,nzx)-1
    step.glob <- 0
    temp.ng.start <- temp.in
    change.prob <- 1
    while(ng > 1){
        iter <- 0
        if (sa) {
            temperature <- temp.ng.start
            temp.iter.gamma <- temp.gamma^(1/iter.max)
        }
        while (iter < iter.max){
            step.glob <- step.glob+1
            iter <- iter +1
            cat(step.glob, "/",step.glob-iter+iter.max*ng, "  ", S.memo[step.glob-1], "\n")
            if(.Platform$OS.type == "windows") flush.console()
            for (i in 1:nzy){
                for (j in 1:nzx){
                    check.rd <- 0
                    check.lu <- 0
                    check.ru <- 0
                    check.ld <- 0
                    if (i<nzy-1){
                        EV.old <- EV
                        EV.dist.old <- EV.dist
                        check.rd <- check.rd+1 
                        check.ld <- check.ld+1   
                        till <- min(nzy,i+ng)
                        rel.cells <- Z0[,j][(i+1):till]
                        mat.to.order <- EV[rel.cells,]
                        vec.dists <- EV.dist[Z0[i,j],rel.cells]
                        ovd <- order(vec.dists)
                        if (any(ovd!= seq(along = vec.dists))){
                            EV[rel.cells,] <- mat.to.order[ovd,]
                            EV.dist[rel.cells,] <- EV.dist[rel.cells[ovd],]
                            EV.dist[,rel.cells] <- EV.dist[,rel.cells[ovd]]
                            rownames(EV.dist)[rel.cells] <- rownames(EV.dist)[rel.cells[ovd]]
                            S.old <- TopoS(EV.dist.old,Cells0.dist)
                            S.new <- TopoS(EV.dist,Cells0.dist)
                            if (S.old > S.new){    
                                if (!sa){EV <- EV.old
                                    EV.dist <- EV.dist.old
                                }
                                else{
                                    simann <- sim.ann(EV.old, EV, EV.dist.old, EV.dist, S.old, S.new, temperature)
                                    EV <- simann$EV
                                    EV.dist <- simann$EV.dist}
                            }
                            if (S.new > S.best){
                                S.best <- TopoS(EV.dist,Cells0.dist)
                                EV.best <- EV
                                EV.dist.best <- EV.dist
                            }
                        }
                    }
                    if (i>2){
                        EV.old <- EV
                        EV.dist.old <- EV.dist   
                        check.lu <- check.lu+1
                        check.ru <- check.ru+1 
                        from <- max(1,i-ng)
                        rel.cells <- Z0[,j][from:(i-1)]
                        mat.to.order <- EV[rel.cells,]
                        vec.dists <- EV.dist[Z0[i,j],rel.cells]
                        ovd <- rev(order(vec.dists))
                        if (any(ovd!=c(1:length(vec.dists)))){
                            EV[rel.cells,] <- mat.to.order[ovd,]
                            EV.dist[rel.cells,] <- EV.dist[rel.cells[ovd],]
                            EV.dist[,rel.cells] <- EV.dist[,rel.cells[ovd]]
                            rownames(EV.dist)[rel.cells] <- rownames(EV.dist)[rel.cells[ovd]]
                            S.old <- TopoS(EV.dist.old,Cells0.dist)
                            S.new <- TopoS(EV.dist,Cells0.dist)
                            if (S.old > S.new){    
                                if(!sa){EV <- EV.old
                                    EV.dist <- EV.dist.old}
                                else{
                                    simann <- sim.ann(EV.old, EV, EV.dist.old, EV.dist, S.old, S.new, temperature)
                                    EV <- simann$EV
                                    EV.dist <- simann$EV.dist}
                            }
                            if (S.new > S.best){
                                S.best <- TopoS(EV.dist,Cells0.dist)
                                EV.best <- EV
                                EV.dist.best <- EV.dist
                            }
                        }
                    }
                    if (j<nzx-1){
                        EV.old <- EV
                        EV.dist.old <- EV.dist
                        check.rd <- check.rd+1
                        check.ru <- check.ru+1
                        till <- min(nzx,j+ng)
                        rel.cells <- Z0[i,][(j+1):till]
                        mat.to.order <- EV[rel.cells,]
                        vec.dists <- EV.dist[Z0[i,j],rel.cells]
                        ovd <- order(vec.dists)                        
                        if (any(ovd!=c(1:length(vec.dists)))){
                            EV[rel.cells,] <- mat.to.order[ovd,]
                            EV.dist[rel.cells,] <- EV.dist[rel.cells[ovd],]
                            EV.dist[,rel.cells] <- EV.dist[,rel.cells[ovd]]
                            rownames(EV.dist)[rel.cells] <- rownames(EV.dist)[rel.cells[ovd]]
                            S.old <- TopoS(EV.dist.old,Cells0.dist)
                            S.new <- TopoS(EV.dist,Cells0.dist)
                            if (S.old > S.new){    
                                if (!sa){EV <- EV.old
                                    EV.dist <- EV.dist.old
                                }
                                else{
                                    simann <- sim.ann(EV.old, EV, EV.dist.old, EV.dist, S.old, S.new, temperature)
                                    EV <- simann$EV
                                    EV.dist <- simann$EV.dist}
                            }
                            if (S.new > S.best){
                                S.best <- TopoS(EV.dist,Cells0.dist)
                                EV.best <- EV
                                EV.dist.best <- EV.dist
                            }
                        }
                    }
                    if (j>2){
                        EV.old <- EV
                        EV.dist.old <- EV.dist
                        check.lu <- check.lu+1
                        check.ld <- check.ld+1
                        from <- max(1,j-ng)
                        rel.cells <- Z0[i,][from:(j-1)]
                        mat.to.order <- EV[rel.cells,]
                        vec.dists <- EV.dist[Z0[i,j],rel.cells]
                        ovd <- rev(order(vec.dists))
                        if (any(ovd!=c(1:length(vec.dists)))){
                            EV[rel.cells,] <- mat.to.order[ovd,]
                            EV.dist[rel.cells,] <- EV.dist[rel.cells[ovd],]
                            EV.dist[,rel.cells] <- EV.dist[,rel.cells[ovd]]
                            rownames(EV.dist)[rel.cells] <- rownames(EV.dist)[rel.cells[ovd]]
                            S.old <- TopoS(EV.dist.old,Cells0.dist)
                            S.new <- TopoS(EV.dist,Cells0.dist)
                            if (S.old > S.new){    
                                if(!sa){EV <- EV.old
                                    EV.dist <- EV.dist.old}
                                else{
                                    simann <- sim.ann(EV.old, EV, EV.dist.old, EV.dist, S.old, S.new, temperature)
                                    EV <- simann$EV
                                    EV.dist <- simann$EV.dist}
                            }
                            if (S.new > S.best){
                                S.best <- TopoS(EV.dist,Cells0.dist)
                                EV.best <- EV
                                EV.dist.best <- EV.dist
                            }
                        }
                    }
                    if (check.rd==2){
                        EV.old <- EV
                        EV.dist.old <- EV.dist
                        count.to.margin <- min((nzy-i),(nzx-j))
                        till <- min(count.to.margin,ng)
                        rel.cells <- diag(Z0[(i+1):(i+till),][,(j+1):(j+till)])
                        mat.to.order <- EV[rel.cells,]
                        vec.dists <- EV.dist[Z0[i,j],rel.cells]
                        ovd <- order(vec.dists)                        
                        if (any(ovd!=c(1:length(vec.dists)))){
                            EV[rel.cells,] <- mat.to.order[ovd,]
                            EV.dist[rel.cells,] <- EV.dist[rel.cells[ovd],]
                            EV.dist[,rel.cells] <- EV.dist[,rel.cells[ovd]]
                            rownames(EV.dist)[rel.cells] <- rownames(EV.dist)[rel.cells[ovd]]
                            S.old <- TopoS(EV.dist.old,Cells0.dist)
                            S.new <- TopoS(EV.dist,Cells0.dist)
                            if (S.old > S.new){    
                                if(!sa){EV <- EV.old
                                    EV.dist <- EV.dist.old}
                                else{
                                    simann <- sim.ann(EV.old, EV, EV.dist.old, EV.dist, S.old, S.new, temperature)
                                    EV <- simann$EV
                                    EV.dist <- simann$EV.dist}
                            }
                            if (S.new > S.best){
                                S.best <- TopoS(EV.dist,Cells0.dist)
                                EV.best <- EV
                                EV.dist.best <- EV.dist
                            }
                        }
                    }
                    if (check.lu==2){
                        EV.old <- EV
                        EV.dist.old <- EV.dist
                        count.to.margin <- min((i-1),(j-1))
                        till <- min(count.to.margin,ng)
                        rel.cells <- diag(Z0[(i-till):(i-1),][,(j-till):(j-1)])
                        mat.to.order <- EV[diag(Z0[(i-till):(i-1),][,(j-till):(j-1)]),]
                        vec.dists <- EV.dist[Z0[i,j],rel.cells]
                        ovd <- rev(order(vec.dists))
                        if (any(ovd!=c(1:length(vec.dists)))){
                            EV[rel.cells,] <- mat.to.order[ovd,]
                            EV.dist[rel.cells,] <- EV.dist[rel.cells[ovd],]
                            EV.dist[,rel.cells] <- EV.dist[,rel.cells[ovd]]
                            rownames(EV.dist)[rel.cells] <- rownames(EV.dist)[rel.cells[ovd]]
                            S.old <- TopoS(EV.dist.old,Cells0.dist)
                            S.new <- TopoS(EV.dist,Cells0.dist)
                            if (S.old > S.new){    
                                if(!sa){EV <- EV.old
                                    EV.dist <- EV.dist.old}
                                else{
                                    simann <- sim.ann(EV.old, EV, EV.dist.old, EV.dist, S.old, S.new, temperature)
                                    EV <- simann$EV
                                    EV.dist <- simann$EV.dist}
                            }
                            if (S.new > S.best){
                                S.best <- TopoS(EV.dist,Cells0.dist)
                                EV.best <- EV
                                EV.dist.best <- EV.dist
                            }
                        }
                    }
                    if (check.ru==2){
                        EV.old <- EV
                        EV.dist.old <- EV.dist
                        count.to.margin <- min((i-1),(nzx-j))
                        till <- min(count.to.margin,ng)
                        rel.cells <- diag(Z0[(i-till):(i-1),][,(j+till):(j+1)])
                        mat.to.order <- EV[rel.cells,]
                        vec.dists <- EV.dist[Z0[i,j],rel.cells]
                        ovd <- rev(order(vec.dists))
                        if (any(ovd!=c(1:length(vec.dists)))){
                            EV[rel.cells,] <- mat.to.order[ovd,]
                            EV.dist[rel.cells,] <- EV.dist[rel.cells[ovd],]
                            EV.dist[,rel.cells] <- EV.dist[,rel.cells[ovd]]
                            rownames(EV.dist)[rel.cells] <- rownames(EV.dist)[rel.cells[ovd]]
                            S.old <- TopoS(EV.dist.old,Cells0.dist)
                            S.new <- TopoS(EV.dist,Cells0.dist)
                            if (S.old > S.new){    
                                if (!sa){EV <- EV.old
                                    EV.dist <- EV.dist.old}
                                else{
                                    simann <- sim.ann(EV.old, EV, EV.dist.old, EV.dist, S.old, S.new, temperature)
                                    EV <- simann$EV
                                    EV.dist <- simann$EV.dist}
                            }
                            if (S.new > S.best){
                                S.best <- TopoS(EV.dist,Cells0.dist)
                                EV.best <- EV
                                EV.dist.best <- EV.dist
                            }
                        }
                    }
                    if (check.ld==2){
                        EV.old <- EV
                        EV.dist.old <- EV.dist
                        count.to.margin <- min((nzy-i),(j-1))
                        till <- min(count.to.margin,ng)
                        rel.cells <- diag(Z0[(i+1):(i+till),][,(j-1):(j-till)])
                        mat.to.order <- EV[rel.cells,]
                        vec.dists <- EV.dist[Z0[i,j],rel.cells]
                        ovd <- order(vec.dists)                        
                        if (any(ovd!=c(1:length(vec.dists)))){
                            EV[rel.cells,] <- mat.to.order[ovd,]
                            EV.dist[rel.cells,] <- EV.dist[rel.cells[ovd],]
                            EV.dist[,rel.cells] <- EV.dist[,rel.cells[ovd]]
                            rownames(EV.dist)[rel.cells] <- rownames(EV.dist)[rel.cells[ovd]]
                            S.old <- TopoS(EV.dist.old,Cells0.dist)
                            S.new <- TopoS(EV.dist,Cells0.dist)
                            if (S.old > S.new){    
                                if(!sa){EV <- EV.old
                                    EV.dist <- EV.dist.old}
                                else{
                                    simann <- sim.ann(EV.old, EV, EV.dist.old, EV.dist, S.old, S.new, temperature)
                                    EV <- simann$EV
                                    EV.dist <- simann$EV.dist}
                            }
                            if (S.new > S.best){
                                S.best <- TopoS(EV.dist,Cells0.dist)
                                EV.best <- EV
                                EV.dist.best <- EV.dist
                            }
                        }
                    }
                }
            }    
            S.memo[step.glob] <- TopoS(EV.dist,Cells0.dist)
            if (iter > 1 || step.glob > 1){
                for (pruf in (1:(step.glob-1))){
                    if (S.memo[pruf]==TopoS(EV.dist,Cells0.dist)){
                        iter <- iter.max
                    }
                }
            }
            if (sa) temperature <- temperature*temp.iter.gamma
        }
        ng <-ng-1
        if (sa) temp.ng.start <- temp.ng.start*temp.gamma
    }
    EV <- EV.best
    EV.dist <- EV.dist.best
    
    Z.old.terms <- Z0

    for (k in 1:(nzx*nzy)){
        mat.test <- matrix(rep(EV0[k,], nzx*nzy), nzx*nzy, ncol(EV0), byrow = TRUE)
        index <- !diag((EV - mat.test) %*% t(EV - mat.test))
        Z.old.terms[Cells0[index, 1], Cells0[index, 2]] <- k
    }
    row.names(EV) <- rownames(EV0)[as.vector(t(Z.old.terms))]
    S.series <- c(S.initial, S.memo[S.memo])    
    preimages <- EV.keep[as.vector(t(Z.old.terms)),]
    if (diss) {
        preimages <- EV.dist
    }
    if (classes[1] == 0) classes <- rep(1, nzx*nzy)
    cl.ord <- classes[as.vector(t(Z.old.terms))]
    results <- list(preimages = preimages, Z = Z0, Z.old.terms = Z.old.terms, 
        cl.ord = cl.ord, S = TopoS(EV.dist, Cells0.dist))
    class(results) <- "EDAM"
    return(results)
}
NaiveBayes <- function (x, ...) 
    UseMethod("NaiveBayes")

NaiveBayes.formula <- function (formula, data, ..., subset, na.action = na.pass) 
{
    call <- match.call()
    Yname <- as.character(formula[[2]])
    if (is.data.frame(data)) {
        m <- match.call(expand = FALSE)
        m$... <- NULL
        m$na.action <- na.action
        m[[1]] <- as.name("model.frame")
        m <- eval(m, parent.frame())
        Terms <- attr(m, "terms")
        if (any(attr(Terms, "order") > 1)) 
            stop("NaiveBayes cannot handle interaction terms")
        Y <- model.extract(m, "response")
        X <- m[ , -attr(Terms, "response")]
        return(NaiveBayes(X, Y, ...))
    }
    else stop("NaiveBayes formula interface handles data frames only")
}


NaiveBayes.default <- function (x, grouping, prior = NULL, usekernel = FALSE, ...) 
{
    x <- data.frame(x)
    if (is.null(prior)) apriori <- table(grouping) / length(grouping)
    else apriori <- as.table(prior / sum(prior))
    call <- match.call()
    #Yname <- deparse(substitute(grouping))
    Yname<- "grouping"
    est <- function(var) 
      if(is.numeric(var)) {
        if (usekernel)
            lapply(split(var, grouping), FUN = function(xx) density(xx, ...))
        else
            cbind(tapply(var, grouping, mean), tapply(var, grouping, sd))
      }
    else {
        tab <- table(grouping, var)
        tab / rowSums(tab)
    }
    tables <- lapply(x, est)
    #for (i in 1:length(tables)) 
    #    {names(dimnames(tables[[i]])) <- c(Yname, colnames(x)[i])}
    names(dimnames(apriori)) <- Yname
    structure(list(apriori = apriori, tables = tables, levels = levels(grouping), 
        call = call, x = x, usekernel = usekernel, varnames = colnames(x)), 
        class = "NaiveBayes")
}
TopoS <- function(EV.dist, Cells.dist){
    dim(Cells.dist) <- NULL
    dim(EV.dist) <- NULL
    beta.est <- (1/(Cells.dist %*% Cells.dist)) %*% Cells.dist %*% EV.dist
    Cells.dist.est <- (Cells.dist * beta.est) - EV.dist
    return(1 - sqrt((Cells.dist.est %*% Cells.dist.est) / (EV.dist %*% EV.dist)))
}

distmirr <- function(dis){
    n <- attr(dis, "Size")
    mirr <- matrix(0, n, n)
    mirr[lower.tri(mirr)] <- dis
    mirr + t(mirr)
}
b.scal <- function(member, grouping, dis = FALSE, eps = 0.0001)
{
betaregion <- function(betaobj, lev, dis = FALSE, eps = 0.0001)
    {
    if(is.null(betaobj)) 
        return(list(CR=0, N=NULL, pm=NULL, NM=NULL, S=0))
    member <- data.matrix(betaobj[,1:length(lev)])
    colnames(member) <- lev
    kdach <- lev[max.col(member)]
    CR <- mean(kdach == betaobj$grouping)
    k <- which(kdach[1] == lev)
    memvec <- member[ , k]
    NT <- length(memvec)
    pm <- mean(memvec)
    S <- if(length(memvec) == 1) 0 else var(memvec)
    memberneu <- member
    if(S < eps)
        {
           NM <- min((pm * (1-pm) / S) - 1, 1e4)
           N <- floor(min(NT, NM))
        }
    else 
        {
        NM <- (pm * (1-pm) / S) - 1
        NU <- floor(min(NT, NM))
        NO <- floor(max(NT, NM))
        if (!dis || CR == 1) N <- NU 
        else # choose optimal N
            {
            N <- NU:NO
            DC <- numeric(length(N))
            for (i in seq(along = N))
                {
                memvecneu <- qbeta(pbeta(memvec, pm * NM, (1-pm) * NM), CR * N[i] , (1-CR) * N[i]) 
                memberneu[,-k] <- ifelse(rep(memvec, length(lev) - 1) > (1 - 1e-6), 
                                        (1 - memvecneu) / (length(lev) - 1), 
                                        member[ , -k] * (1 - memvecneu) / (1 - memvec))
                memberneu[,k] <- memvecneu
                ergi <- ucpm(memberneu, betaobj$grouping)
                DC[i] <- ergi$CR * NT + ergi$AC ### optimize CR+AC!!!!
                }
            N <- N[which.max(DC)]
            }
        }
    return(list(CR=CR, N=N, pm=pm, NM=NM, S=S))
    }  
lev <- levels(grouping)
membercheck(member)
memberdata <- data.frame(member, grouping = grouping)
memberlist <- vector(length(lev), mode = "list")
names(memberlist) <- as.character(1:length(lev))
mcol <- max.col(member)
if(!all(lev %in% mcol))
    stop("At least one observation is required for each predicted class.")
temp <- split(memberdata, mcol)
for(i in names(temp))
    memberlist[i] <- temp[i]
res <- list()
res$model <- lapply(memberlist, betaregion, lev=lev, dis=dis, eps=eps)
res$eps <- eps
res$member <- betascale(res, member)
return(res)
}



betascale <- function(betaobj, member)
{
if (missing(member)) member <- betaobj$member
else 
    {
    membercheck(member)
    eps <- betaobj$eps
    betaobj <- betaobj$model
    lev <- names(betaobj)
    member <- as.matrix(member)
    memberalt <- member
    for (i in seq(along = lev))
        {
        ti <- which(max.col(memberalt) == lev[i])
        memvec <- member[ti, i]
        if(!length(memvec))
            next
        bets <- betaobj[[i]]
        if (bets$S < eps || bets$CR == 1) {
            memvecneu <- memvec - bets$pm + bets$CR
            memvecneu[memvecneu > 1] <- 1
        }
        else{
            memvecneu <- with(bets, qbeta(pbeta(memvec, pm * NM, (1 - pm) * NM), 
                            CR * N, (1 - CR) * N))
        }
        member[ti,-i] <- matrix(ifelse(rep(memvec, length(lev)-1) > (1-1e-6), 
                                (1-memvecneu)/(length(lev)-1), 
                                member[ti,-i]*(1-memvecneu)/(1-memvec)),
                         ncol = length(lev)-1)
        member[ti,i] <- memvecneu
        }
    }
return(member)
}
calc.trans <- function(x)
{
# calculates the transition matrix given the time series of states
# ----------------------------------------------------------------
# x: vector of states
    x <- factor(x)
    if(length(x) > 1){
        tbl <- table(x[1:(length(x)-1)], x[2:length(x)])
        return(tbl / rowSums(tbl))
    }
    else return(table(x[1], x[1]) / sum(table(x[1])))
}
centerlines<-function(n)
{
li<-NULL
for (i in 2:n)
    for (j in 1:(i-1))
        {
        dummy<-numeric(n)
        dummy[c(i,j)]<-0.5
        li<-rbind(li,rep(1/n,n),dummy,rep(NA,n))
        }
rownames(li)<-NULL
return(li)
}
dkernel<-function(x, kernel=density(x), interpolate=FALSE)
{
foo<-function(x,kernel,n)
{
if (x <= kernel$x[1]) return(kernel$y[1])
else if (x >= kernel$x[n]) return(kernel$y[n])
else 
    {
        pos<-which.min((kernel$x-x)^2)
        if (kernel$x[pos]>x) return(mean(c(kernel$y[pos],kernel$y[(pos-1)])))
        else return(mean(c(kernel$y[pos],kernel$y[(pos+1)])))
    }
}
n<-length(kernel$x)
if (interpolate) y<-sapply(x,foo,kernel=kernel,n=n)
else y<-sapply(x,FUN=function(y){kernel$y[(which.min((kernel$x-y)^2))]})
return(y)
}
e.scal<-function(x, k = 1, tc = NULL)
{
    scal <- function(k, x)
    {
        xx <- exp(k * x)
        return(xx / rowSums(xx))
    }
    scal2 <- function(k, x, tc)
    {
        xx <- scal(k, x)
        ec <- factor(max.col(xx), levels = seq(along = colnames(xx)), 
            labels = colnames(xx))
        return(mean(ec != tc))
    }
    if(!is.null(tc)) 
        k <- optimize(scal2, c(0, 1000), x = x, tc = tc)$minimum
    x <- scal(k, x)
    return(list(sv = x, k = k))
}
"errormatrix" <-
function(true, predicted, relative=FALSE)
# if `relative=TRUE', rows are nomalized by their sum,
# the last row is is divided by the total number of misclassifications,
# and the lower right cell is the total misclassification rate.
{
  stopifnot(length(true)==length(predicted))
  tnames <- 
    if(is.factor(true)) levels(true) 
    else unique(true)
  pnames <- 
    if(is.factor(predicted)) levels(predicted) 
    else unique(predicted)
  allnames <- sort(union(tnames, pnames))
  n <- length(allnames)
  true <- factor(true, levels = allnames)
  predicted <- factor(predicted, levels = allnames)
  tab <- table(true, predicted)
  mt <- tab * (matrix(1, ncol = n, nrow = n) - diag( , n, n))
  rowsum <- rowSums(mt)
  colsum <- colSums(mt)
  result <- rbind(cbind(tab, rowsum), c(colsum, sum(colsum)))
  dimnames(result) <- list("true" = c(allnames, "-SUM-"), 
    "predicted" = c(allnames, "-SUM-"))
  if(relative){
    total <- sum(result[1:n, 1:n])
    # normalize last row:
    n1 <- n + 1
    result[n1, 1:n] <- 
        if(result[n1, n1] != 0) result[n1, 1:n] / result[n1, n1] 
        else 0
    # normalize remaining matrix:
    rownorm <- function(Row,Length)
    { return( if(any(Row[1:Length]>0)) Row/sum(Row[1:Length])
              else rep(0,Length+1) ) 
    }
    result[1:n,] <- t(apply(result[1:n,], 1, rownorm, Length=n))
    # normalize lower right cell:
    result[n1, n1] <- result[n1, n1] / total
  }
  return(result)
}
hmm.sop <- function(sv, trans.matrix, prob.matrix)
{
# sv: Start value
# trans.matrix: probability matrix of transitions
# prob.matrix: p(c|x)
  rekurs <- prob.matrix
  rekurs[1,] <- prob.matrix[1, ] * trans.matrix[sv, ]
  for(z in 2:nrow(prob.matrix))
    for(j in 1:ncol(prob.matrix))
        rekurs[z,j] <- 
          prob.matrix[z, j] * trans.matrix[, j] %*% rekurs[(z-1), ]
  rekurs <- rekurs / rowSums(rekurs)
  return(rekurs)
}
membercheck <- function(member){
    member <- as.matrix(member)
    if(any(member) < 0 || any(member) > 1) 
        stop("Membership values (posterior probabilities) must be in [0,1].")
    if(!identical(all.equal(rowSums(member), rep(1, nrow(member))), TRUE))
        stop("Membership values (posterior probabilities) must sum up to 1.")
}
partimat<-function (x, ...) 
    UseMethod("partimat")


partimat.default <- function(x, grouping, method = "lda", prec = 100, 
    nplots.vert, nplots.hor, main = "Partition Plot", name, mar,
    plot.matrix = FALSE, ...){
  
    nvar <- ncol(x)
    if(nvar < 2) stop("at least 2 variables required")
    if(nlevels(grouping) < 2) stop("at least two classes required")
    nobs <- nrow(x)
    if(missing(name)) name <- colnames(x)
    # plot in scatterplot matrix
    if(plot.matrix){
        plot.new()
        if(missing(mar)) mar <- rep(0, 4)
        opar <- par(mfrow = c(nvar, nvar), mar = mar, oma = rep(3, 4), xpd = NA)
        on.exit(par(opar))
        for (i in 2:nvar)
            for (j in 1:(i-1))
            {
                par(mfg = c(i, j))
                drawparti(grouping, x[,j], x[,i], method = method, 
                    prec = prec, legend.err = plot.matrix, xaxt="n", 
                    yaxt="n", xlab="", ylab="", ...)
                if(j == 1) axis(2)
                if(i == nvar) axis(1)
                par(mfg = c(j, i))
                drawparti(grouping, x[,i], x[,j], method = method, 
                    prec = prec, legend.err = plot.matrix, xaxt="n", 
                    yaxt="n", xlab="", ylab="", ...)
                if(j == 1) axis(3) 
                if(i == nvar) axis(4)
            }
        for (i in 1:nvar)
            {
                par(mfg = c(i, i))
                plot(x[,i], x[,i], type = "n", xaxt="n", yaxt="n", xlab="", ylab="")
                if(i == 1){
                    axis(2); axis(3)
                }
                else if(i == nvar){
                    axis(1); axis(4)
                }
                mxi <- mean(range(x[,i]))
                text(mxi, mxi, name[i], ...)
            }
    } # Plotting fuzzy in rows and columns to save space:
    else{
        ncomb <- round(0.5 * nvar * (nvar-1))
        if (missing(nplots.hor) && missing(nplots.vert)){
            nplots.hor<-ceiling(sqrt(ncomb))
            nplots.vert<-floor(sqrt(ncomb))
        }
        else if (missing(nplots.hor)) nplots.hor<-ceiling(ncomb/nplots.vert)
        else if (missing(nplots.vert)) nplots.vert<-ceiling(ncomb/nplots.hor)
        vars <- matrix(ncol=ncomb,nrow=2*nobs)
        varname <- matrix(ncol=ncomb,nrow=2)
        k <- 1
        for (i in 2:nvar)
            for (j in 1:(i-1))
            {
                vars[,k] <- c(x[,i], x[,j])
                varname[,k] <- c(name[i], name[j])
                k <- k + 1
            }

        if(missing(mar)) mar <- c(5.1, 4.1, 2.1, 1.1)
        opar <- par(mfrow = c(nplots.vert, nplots.hor), mar = mar, 
            oma = c(0, 0, !is.null(main), 0))
        on.exit(par(opar))
        
        sapply(1:ncomb, function(k) 
            drawparti(grouping = grouping, x = vars[(1:nobs), k], 
                y = vars[(nobs+1):(2*nobs), k], method = method, 
                xlab = varname[1,k], ylab = varname[2,k], prec = prec, 
                legend.err = plot.matrix, ...)
        )
        par(mfrow=c(1,1))
        title(main = main, outer = TRUE)
    }
}

partimat.formula<-function (formula, data = NULL, ..., subset, na.action = na.fail) 
{
    m <- match.call(expand.dots = FALSE)
    if (is.matrix(eval.parent(m$data))) 
        m$data <- as.data.frame(data)
    m$... <- NULL
    m[[1]] <- as.name("model.frame")
    m <- eval.parent(m)
    Terms <- attr(m, "terms")
    grouping <- model.response(m)
    x <- model.matrix(Terms, m)
    xvars <- as.character(attr(Terms, "variables"))[-1]
    if ((yvar <- attr(Terms, "response")) > 0) 
        xvars <- xvars[-yvar]
    xlev <- if (length(xvars) > 0) {
        xlev <- lapply(m[xvars], levels)
        xlev[!sapply(xlev, is.null)]
    }
    xint <- match("(Intercept)", colnames(x), nomatch = 0)
    if (xint > 0) 
        x <- x[, -xint, drop = FALSE]
    res <- partimat.default(x, grouping, ...)
    res$terms <- Terms
    cl <- match.call()
    cl[[1]] <- as.name("partimat")
    res$call <- cl
    res$contrasts <- attr(x, "contrasts")
    res$xlevels <- xlev
    attr(res, "na.message") <- attr(m, "na.message")
    if (!is.null(attr(m, "na.action"))) 
        res$na.action <- attr(m, "na.action")
    res
    invisible()
}

partimat.matrix<-function (x, grouping, ..., subset, na.action = na.fail) 
{
    if (!missing(subset)) {
        x <- x[subset, , drop = FALSE]
        grouping <- grouping[subset]
    }
    if (!missing(na.action)) {
        dfr <- na.action(structure(list(g = grouping, x = x), 
            class = "data.frame"))
        grouping <- dfr$g
        x <- dfr$x
    }
    res <- partimat.default(x, grouping, ...)
    cl <- match.call()
    cl[[1]] <- as.name("partimat")
    res$call <- cl
    res
    invisible()
}

partimat.data.frame<-function (x, ...) 
{
   res <- partimat.matrix(structure(data.matrix(x), class = "matrix"), 
        ...)
    cl <- match.call()
    cl[[1]] <- as.name("partimat")
    res$call <- cl
    res
    invisible()    
}

drawparti <- function(grouping, x, y, method = "lda", prec = 100, 
    xlab=NULL, ylab=NULL, col.correct = "black", col.wrong = "red", 
    col.mean = "black", col.contour = "darkgrey", gs = as.character(grouping), 
    pch.mean = 19, cex.mean = 1.3, print.err = 0.7, legend.err = FALSE,
    legend.bg = "white", imageplot = TRUE, image.colors = cm.colors(nc), ...){
    #grouping: class vec.
    #x: first data vec.
    #y: second data vec.
    #prec: nr. of hor/vert splits.
    
    z <- switch(method,
        lda = lda(grouping ~ x + y,...),
        qda = qda(grouping ~ x + y,...),
        svmlight = svmlight(grouping ~ x + y,...),
        rda = rda(grouping~ x + y, 
            data = cbind.data.frame("grouping" = grouping, "x" = x, "y" = y), ...),
        sknn = sknn(grouping ~ x + y,...),
        rpart = rpart(grouping~ x + y,...),
        disco = disco(grouping ~ x + y,
            data = cbind.data.frame("grouping" = grouping, "x" = x, "y" = y), ...),
        naiveBayes = naiveBayes(grouping~ x + y, 
            data = cbind.data.frame("grouping" = grouping, "x" = x, "y" = y), ...),
        error("method not yet supported"))

    # Build a grid on the 2 coordinates
    xg <- seq(min(x), max(x), length = prec)
    yg <- seq(min(y), max(y), length = prec)
    grd <- expand.grid(x = xg, y = yg)
    # Calcultate posterior Probabilities on grid points
    temp <- switch(method,
        lda = predict(z, grd,...)$post,
        qda = predict(z, grd,...)$post,
        svmlight = e.scal(predict(z, grd,...)$post)$sv,
        rda = predict(z, grd, posterior=TRUE, aslist=TRUE)$post,
        rpart = predict(z, grd, ...),
        sknn = predict(z, grd, ...)$post,
        disco = predict(z, grd, ...)$post,
        naiveBayes = predict(z, grd , type="raw", ...),
        error("method not yet supported"))
    khead<- switch(method,
        lda = predict(z, data.frame(cbind(x,y)),...)$class,
        qda = predict(z, data.frame(cbind(x,y)),...)$class,
        svmlight = predict(z, data.frame(cbind(x,y)),...)$class,
        rda = predict(z, data.frame(cbind(x,y)), posterior=TRUE, aslist=TRUE)$class,
        rpart = predict(z, data.frame(cbind(x,y)), type="class", ...),
        sknn = predict(z, data.frame(cbind(x,y)),...)$class,
        disco = predict(z, data.frame(cbind(x,y)),...)$class,
        naiveBayes = predict(z, data.frame(cbind(x,y)), ...),
        error("method not yet supported"))

    colorw <- grouping != khead
    err <- round(mean(colorw), 3)
    color <- ifelse(colorw, col.wrong, col.correct)

    nc <- ncol(temp)
    if(imageplot){
        image(xg, yg, matrix(apply(temp, 1, which.max), ncol = prec), 
            main = NULL, col = image.colors, breaks = (0:nc) + .5, 
            xlab = xlab, ylab = ylab, ...)
        points(x, y, pch = gs, col = color, ...)
        box()
    }
    else 
        plot(x, y, pch = gs, col = color, main = NULL, xlab = xlab, ylab = ylab, ...)
    if((method=="lda") || (method=="qda")) 
        points(z$means, pch = pch.mean, cex = cex.mean, col = col.mean)
    
    # For each class calculate the difference between prob. and max(prob) for other class,
    # so, the obs is assigned to class iff diff>0
    if(!imageplot)
        for(i in 1:ncol(temp)){
            dummy <- temp[,i] - apply(temp[ , -i, drop = FALSE], 1, max)            
        # Draw contour line at hight=0, i.e. class border
            contour(xg, yg, matrix(dummy, ncol = prec), levels = 0, 
                add = TRUE, drawlabels = FALSE, col = col.contour)
        }

    if(print.err){
        if(legend.err)
            legend(par("usr")[1], par("usr")[4], 
                legend = paste("Error:", err), bg = legend.bg, cex = print.err)
        else
            mtext(paste("app. error rate:", err), 3, cex = print.err)
    }

}
plot.NaiveBayes<-function(x, vars, n=1000, legendplot=TRUE,
                            lty=1:length(x$apriori), col=rainbow(length(x$apriori)),
                            ylab="Density",main="Naive Bayes Plot",...)
{
if (missing(vars)) vars<-names(x$tables)
vars<-vars[is.element(vars,names(x$tables))]
if (length(vars>0))
for (j in 1:length(vars))
{
    dummy<-(x$tables[names(x$tables)==as.name(vars[j])])
    if(class(dummy[[1]])=="matrix")
    {
        dummy<-data.frame(dummy)
        plotvector<-seq(min(x$x[,vars[j]]),max(x$x[,vars[j]]),len=n)
        pv<-matrix(0,nrow=nrow(dummy),ncol=n)
        for (i in 1:nrow(dummy))
        pv[i,]<-dnorm(plotvector,mean=dummy[i,1],sd=dummy[i,2])*x$apriori[i]
    plot(plotvector,pv[1,],type="l",lty=lty[1],ylim=c(0,max(pv)),xlab=vars[j],ylab=ylab,col=col[1], main=main,...)
    for(i in 2:nrow(dummy))
        lines(plotvector,pv[i,],lty=lty[i], col=col[i],...)
    if (legendplot) legend(min(plotvector),max(pv),legend=rownames(dummy),lty=lty,col=col)
    }
    if(class(dummy[[1]])=="table")
    {
    mosaicplot(dummy[[1]], main=main, ...)
    }
    if(class(dummy[[1]])=="list")
    {        
        plotvector<-seq(min(x$x[,vars[j]]),max(x$x[,vars[j]]),len=n)
        pv<-matrix(0,nrow=length(dummy[[1]]),ncol=n)
        for (i in 1:length(dummy[[1]]))
            pv[i,]<-dkernel(plotvector,kernel=dummy[[1]][[i]])*x$apriori[i]
        plot(plotvector,pv[1,],type="l",lty=lty[1],ylim=c(0,max(pv)),xlab=vars[j],ylab=ylab,col=col[1], main=main,...)
        for(i in 2:nrow(pv))
            lines(plotvector,pv[i,],lty=lty[i], col=col[i],...)
    if (legendplot) legend(min(plotvector),max(pv),legend=names(dummy[[1]]),lty=lty,col=col)
    }
    }
}
plot.EDAM <- function(...) shardsplot(...)
predict.NaiveBayes<-function (object, newdata, 
    threshold = 0.001, ...) 
{
    if (missing(newdata)) newdata<-object$x
    if (!any(is.null(colnames(newdata)),is.null(object$varnames))) { # (both colnames & varnames are given) 
        if(all(is.element(object$varnames,colnames(newdata)))){        # (varnames is a subset of colnames)   
            newdata <- data.frame(newdata[,object$varnames])
        }
    }
    
    nattribs <- ncol(newdata)
    isnumeric <- sapply(newdata, is.numeric)
    L <- sapply(1:nrow(newdata), function(i) {
        ndata <- as.numeric(newdata[i, ])
        L <- object$apriori * apply(sapply(1:nattribs, function(v) {
            nd <- ndata[v]
            if (is.na(nd)) 
                rep(1, length(object$apriori))
            else {
                prob <- if (isnumeric[v]) {
                  msd <- object$tables[[v]]
                  if (object$usekernel) sapply(msd,FUN=function(y){dkernel(x=nd,kernel=y,...)})
                  else dnorm(nd, msd[, 1], msd[, 2])
                }
                else object$tables[[v]][, nd]
                prob[prob == 0] <- threshold
                prob
            }
        }), 1, prod)
        L/sum(L)
    })
    posterior<-t(L)
    classdach<-factor(object$levels[apply(L, 2, which.max)], levels = object$levels)
    colnames(posterior)<-object$levels
    rownames(posterior)<-names(classdach)<-rownames(newdata)
    return(list(class=classdach,posterior=posterior))
}
quadtrafo<-function(e, f=NULL, g=NULL, h=NULL)
{
require(scatterplot3d)
if (is.matrix(e)) dat <- e
else    if (is.null(f)) dat <- t(e)
        else dat <- cbind(e,f,g,h)
result <- dat %*%   
    matrix(c(1, 0.5, 0.5, 0, 
            0, sqrt(3)/2, sqrt(3)/6, 0, 
            0, 0, sqrt(6)/3, 0), nrow = 4)
            
colnames(result)<-c("x","y","z")
return(result)
}


quadlines <- function(e, f=NULL, g=NULL, h=NULL, sp=s3d, ...)
{
  result <- quadtrafo(e,f,g,h)
  sp$points3d(result[,1], result[,2], result[,3], type="l", ...)
  invisible(result)
}


quadpoints <- function(e, f=NULL, g=NULL, h=NULL, sp=s3d, ...)
{
  result <- quadtrafo(e,f,g,h)
  sp$points3d(result[,1], result[,2], result[,3], type="p", ...)
  invisible(result)
}




quadplot <- function(e=NULL, f=NULL, g=NULL, h=NULL, angle=75, scale.y=0.6, 
                    label=1:4, labelcol=rainbow(4), labelpch=19, 
                    labelcex=1.5, main="", s3d.control = list(), 
                    simplex.control = list(), legend.control = list(),
                    ...)
{
corners <- quadtrafo(diag(4))
s3d <- do.call("scatterplot3d", 
    c(list(0.5, 0.2886751, 0.2041241, type="n", 
        xlim=range(corners[,1]), ylim=range(corners[,2]), zlim=range(corners[,3]),
        axis = FALSE, grid=FALSE, angle = angle, scale.y = scale.y, main=main),
      s3d.control)
)

# We need to get the correct "mar" values set in the first scatterplot3d() call:
par(mar = get("mar", pos=environment(s3d$xyz.convert)))
# We don't want to waste too much space:
usr <- as.vector(sapply(s3d$xyz.convert(corners), range))
par(usr = usr, xpd = NA)

# outer simplex:
do.call("quadlines", 
    c(list(e = diag(4)[c(1:4,1,3,2,4),], sp = s3d), 
      simplex.control)
)
do.call("quadpoints", 
    c(list(e = diag(4), sp = s3d, pch = labelpch, col = labelcol, cex = labelcex)))
# legend:
do.call("legend", 
    c(list(usr[1], usr[4], legend = label, col = labelcol, pch = labelpch, cex = labelcex),
      legend.control)
)
# points:
if (!is.null(e)) 
    quadpoints(e, f, g, h, sp = s3d, ...)
return(s3d)
}
rda <- function(x, ...) 
{
  UseMethod("rda")
}



rda.formula <- function(formula, data = NULL, ...) 
{
    m <- match.call(expand.dots = FALSE)
    if (is.matrix(eval.parent(m$data))) 
        m$data <- as.data.frame(data)
    m$... <- NULL
    m[[1]] <- as.name("model.frame")
    m <- eval.parent(m)
    Terms <- attr(m, "terms")
    grouping <- model.response(m)
    x <- model.matrix(Terms, m)
    xvars <- as.character(attr(Terms, "variables"))[-1]
    if ((yvar <- attr(Terms, "response")) > 0) 
        xvars <- xvars[-yvar]
    xlev <- if (length(xvars) > 0) {
        xlev <- lapply(m[xvars], levels)
        xlev[!sapply(xlev, is.null)]
    }
    xint <- match("(Intercept)", colnames(x), nomatch = 0)
    if (xint > 0) 
        x <- x[, -xint, drop = FALSE]
    res <- rda.default(x, grouping, ...)
    res$terms <- Terms
    cl <- match.call()
    cl[[1]] <- as.name("rda")
    res$call <- cl
    res$contrasts <- attr(x, "contrasts")
    res$xlevels <- xlev
    attr(res, "na.message") <- attr(m, "na.message")
    if (!is.null(attr(m, "na.action"))) 
        res$na.action <- attr(m, "na.action")
    res
}




rda.default <- function(x, grouping=NULL, prior=NULL, 
                        gamma=NA, lambda=NA, regularization=c("gamma"=gamma, "lambda"=lambda),
                        crossval=TRUE, fold=10, train.fraction=0.5, 
                        estimate.error=TRUE, output=FALSE,
                        startsimplex=NULL, max.iter=100, trafo=TRUE,
                        simAnn=FALSE, schedule=2, T.start=0.1, halflife=50, zero.temp=0.01, alpha=2, K=100, ...)
{
  classify <- function(dataset, mu, sigma, pooled, gamma, lambda, g, p, n, prior)
  # mu            : matrix with group means as columns 
  # sigma         : (p x p x g)-Array of group covariances 
  # pooled        : pooled covariance matrix 
  # gamma, lambda : regularization parameters 
  # g, p, n       : numbers of classes, variables & observations
  {
    # compute likelihood of each datum for each group: 
    likelihood <- matrix(NA, nrow=n, ncol=g)
    dataset <- matrix(dataset,ncol=p)
    i <- 1
    singu <- FALSE
    while ((i <= g) & (!singu)){ # compute likelihoods group-wise
      reg.cov <- (1-lambda)*sigma[,,i] + lambda*pooled                   # shift towards pooled cov. 
      reg.cov <- if (is.matrix(reg.cov)) # cov is a matrix
                   (1-gamma)*reg.cov + gamma*mean(diag(reg.cov))*diag(p) # shift towards identity 
                 else # one variable  =>  cov is 1x1-matrix:
                   (1-gamma)*reg.cov + gamma*reg.cov                     # shift towards identity 
      # try to invert covariance matrix: 
      inv.trial <- try(solve(reg.cov), silent=TRUE)
      singu <- ((! is.matrix(inv.trial)) || (!is.finite(ldc<-log(1/det(inv.trial)))))
      if (!singu) likelihood[,i] <- prior[i]*.dmvnorm(dataset, mu=mu[,i], inv.sigma=inv.trial, logDetCov=ldc)
      i <- i+1
    }
    if (!singu) result <- apply(likelihood,1,function(x){order(x)[g]}) # return classifications
    else result <- rep(0,nrow(dataset)) # return definitely false classifications 
    return(result)
  }

  goalfunc <- function(paramvec=c(0.5,0.5)) 
  # depends on a (2-dim.) parameter vector, so it can be handled by minimization function. 
  # first element: gamma,  second element: lambda. 
  # returns mean misclassification rate for the bootstrap samples. 
  {
    error.rates <- numeric(fold)
    for (i in 1:fold) {
      siggi <- array(covariances[,,,i],c(p,p,g))
      mumu <- array(means[,,i], c(p,g))
      prediction <- classify(data[-train[,i],], 
                             mu=mumu, sigma=siggi, pooled=covpooled[,,i],
                             gamma=paramvec[1], lambda=paramvec[2], g=g, p=p, n=n-sum(train[,i]!=0), prior=prior)
      errors <- prediction != grouping[-train[,i]]
      group.rates <- tabulate(grouping[-train[,i]][errors],g) / test.freq[,i] #error rate by group
      group.rates[test.freq[,i]==0] <- (g-1)/g # conservative estimate for "empty" groups 
      error.rates[i] <- t(prior)%*%group.rates
    }
    return(mean(error.rates))
  }
  
  nelder.mead <- function(func, startsimplex, mini=NULL, maxi=NULL, link=function(x)(x),
                          a=1, b=0.5, g=3, r=0.5, e=.Machine$double.eps^0.5, 
                          max.iter=1e6, best.possible=-Inf, output=FALSE, 
                          simAnn=TRUE, schedule=1, T.start=NULL, halflife=200, zero.temp=0.01, alpha=2, K=500, degen=1)
  # minimizes the function `func' with several real parameters. 
  # Parameters: 
  #  func          :  the function to be minimized; must require a vector of (at least two) real-value-parameters 
  #  startsimplex  :  either a starting simplex or just the number of dimensions (for the real parameters) 
  #  mini, maxi    :  vectors of lower & upper bounds for the parameters 
  #  a             :  a > 0      `reflection coefficient' 
  #  b             :  0 < b < 1  `contraction coefficient' 
  #  g             :  g > 1      `expansion coefficient' 
  #  r             :  0 < r < 1  `reduction coefficient' 
  #  e             :  e > 0      `epsilon' (for stop criterion) 
  #  max.iter      :  (<= Inf) maximum number of iterations 
  #  best.possible :  best possible value; algorithm stops if reached. 
  #  output        :  flag for text output during calculation 
  # Parameters for Simulated Annealing: 
  #  simAnn        :  indicates whether Simulated Annealing should be used 
  #  T.start       :  starting temperature for Simulated Annealing 
  #  schedule      :  annealing schedule 1 or 2 (exponential or polynomial) 
  #  halflife      :  number of iterations until temperature is reduced to a half    (schedule I) 
  #  zero.temp     :  temperature at which it is set to zero                         (schedule I) 
  #  alpha         :  power of temperature reduction (linear, quadratic, cubic,...)  (schedule II) 
  #  K             :  number of iterations until temperature = 0                     (schedule II) 
  #  degen         :  number to substract or add to mini, maxi in degnerated simplices 
  {
    rmunif<-function(n,min=0,max=1)
    {
        dim<-length(min)
        erg<-matrix(0,nrow=n,ncol=dim)
        for (i in 1:dim) erg[,i]<-runif(n,min[i],max[i])
        if (n==1) erg<-as.vector(erg)
        return(erg)
    }
    restrict <- function(x) # forces x to be within its bounds (given by `mini' & `maxi') 
    {                       # (necessary for reflexion & expansion) 
      new.x <- x
      new.x[x<mini] <- mini[x<mini]
      new.x[x>maxi] <- maxi[x>maxi]
      return(new.x)
    }
    max.sample <- function(f) # if there is more than one worst point, one of these is sampled. 
    {
      i <- which(f==max(f))
      if (length(i)>1) i <- sample(i,1)
      return(i)
    }
    min.sample <- function(f)
    {
      i <- which(f==min(f))
      if (length(i)>1) i <- sample(i,1)
      return(i)
    }
    rlog <- function(n=1) # returns a "logarithmically distributed" random number 
    {                     # (actually equals an Exp(1)-distributed RV)            
      return(-log(runif(n)))
    }

    if (is.vector(startsimplex) && length(startsimplex)==1 && startsimplex>0 && startsimplex==trunc(startsimplex))
      startsimplex <- rbind(rep(-sqrt(((sqrt(1.5)-sqrt(0.5))^2)/2),startsimplex),diag(startsimplex))   
    p <- ncol(startsimplex)                  #  p = number of parameters 
    if (is.null(mini)) mini <- rep(-Inf,p)
    if (is.null(maxi)) maxi <- rep(Inf,p)
   
    minir<-mini
    maxir<-maxi
    for (i in 1:length(mini))
        {
        if ((mini[i]==-Inf) && (maxi[i]==Inf)) {minir[i]<- -degen; maxir[i]<- degen}
        else if ((mini[i]==-Inf) && (maxi[i]!=Inf)) minir[i]<- maxi[i]-degen
        else if ((mini[i]!=-Inf) && (maxi[i]==Inf)) maxir[i]<- mini[i]+degen
        }
    startsimplex <- t(apply(startsimplex,1,restrict))
    
    simplex<-startsimplex
    if (output) {
      if (simAnn) {
        if (schedule==1) cat("Performing Simulated Annealing with `exponential' cooling schedule",paste("(halflife=",halflife,").\n", sep=""))
        else {cat("Performing Simulated Annealing with `polynomial' cooling schedule",paste("(alpha=",alpha,").\n",sep=""))
              cat("Zero temperature after",K,"iterations.\n")}
      }
      else cat("Performing Nelder-Mead minimization.\n")
      cat("Calculations for starting simplex...\n")
      if(.Platform$OS.type == "windows") flush.console()
    }
    f <- apply(simplex,1,function(x){func(link(x))})  #  vector of `func'-values, corresponding to the simplex 
    min.index <- min.sample(f)
    if (output) {
      cat("best/worst value:",min(f),"/",max(f),"\n")
      if(.Platform$OS.type == "windows") flush.console()
    }
    best.ever <- c(f[min.index], simplex[min.index,])
    close.enough <- ((!simAnn) && ((sd(f) <= e) || any(f <= best.possible)))
    if (simAnn) {
      if (is.null(T.start)) T.start <- min(max(f)-best.possible, max(f))
      temp <- T.start
      if (schedule==1) faktor <- 0.5^(1/halflife)
      if ((output) && (schedule==1)){ 
        cat("Zero temperature after",ceiling(log(zero.temp/T.start,base=faktor)),"iterations.\n")
        if(.Platform$OS.type == "windows") flush.console()
      }
    }
    i <- 1
    while ((!close.enough) && (i<=max.iter)) {  # START iterations 
      if ((simAnn) && (temp>0)) f.hot <- f + temp*rlog(p+1)    # the "annealing" function values 
      else f.hot <- f
      min.index <- min.sample(f.hot)
      max.index <- max.sample(f.hot)
      centroid <- apply(simplex[-max.index,], 2, mean)  # centroid of simplex except worst point 
     # REFLEXION: 
      x.reflex <- restrict(centroid + a*(centroid-simplex[max.index,]))
      if (output) {
        cat("Reflexion\n")
        if(.Platform$OS.type == "windows") flush.console()
      }
      f.reflex <- func(link(x.reflex))
      if (f.reflex<best.ever[1]) best.ever <- c(f.reflex, x.reflex)
      if ((simAnn) && (temp>0)) f.reflex.hot <- f.reflex - temp*rlog(1)
      else f.reflex.hot <- f.reflex
     # If reflexion yielded improvement: EXPANSION. 
      if (f.reflex.hot <= f.hot[min.index]) {
        x.expand <- restrict(centroid + g*(x.reflex - centroid))
        if (output) { 
          cat("Expansion\n")
          if(.Platform$OS.type == "windows") flush.console()
        }
        f.expand <- func(link(x.expand))
        if (f.expand<best.ever[1]) best.ever <- c(f.expand, x.expand)        
        if ((simAnn) && (temp>0)) f.expand.hot <- f.expand - temp*rlog(1)
        else f.expand.hot <- f.expand
       # If `expanded point' is best, it is saved, 
        if (f.expand.hot < f.hot[min.index]) {
          simplex[max.index,] <- x.expand
          f[max.index] <- f.expand
          f.hot[max.index] <- f.expand.hot
          min.index <- max.index
          max.index <- max.sample(f.hot)
        }
       # otherwise `reflected point' is saved. 
        else {
          simplex[max.index,] <- x.reflex
          f[max.index] <- f.reflex
          f.hot[max.index] <- f.reflex.hot
          min.index <- max.index
          max.index <- max.sample(f.hot)
        }
      }
     # If `reflected point' wasn't best, but better than/equal to any other point (except maximum), it is saved. 
      else if (any(f.reflex.hot <= f.hot[-max.index])) {
             simplex[max.index,] <- x.reflex
             f[max.index] <- f.reflex
             f.hot[max.index] <- f.reflex.hot
             max.index <- max.sample(f.hot)
           }
           else { # If `reflected point' was only better than worst point (but worse than others), it is saved as well. 
             if (f.reflex.hot <= f.hot[max.index]) {
               simplex[max.index,] <- x.reflex
               f[max.index] <- f.reflex
               f.hot[max.index] <- f.reflex.hot
             }
            # Since reflexion didn't improve: CONTRACTION. 
             x.contract <- centroid + b*(simplex[max.index,]-centroid)
             if (output) {
               cat("Contraction\n")
               if(.Platform$OS.type == "windows") flush.console()
             }
             f.contract <- func(link(x.contract))
             if (f.contract<best.ever[1]) best.ever <- c(f.contract, x.contract)
             if ((simAnn) && (temp>0)) f.contract.hot <- f.contract - temp*rlog(1)
             else f.contract.hot <- f.contract
            # If `contracted' is better than (or equal to) maximum, the maximum is replaced 
             if (f.contract.hot <= f.hot[max.index]) {
               simplex[max.index,] <- x.contract
               f[max.index] <- f.contract
               f.hot[max.index] <- f.contract.hot
               min.index <- min.sample(f.hot)
               max.index <- max.sample(f.hot)
             }
             else {
              # If `contracted' was worse (greater) than maximum, the whole simplex is REDUCED: 
               for (j in (1:(p+1))[-min.index]) 
                 simplex[j,] <- simplex[min.index,]+r*(simplex[j,]-simplex[min.index,])
               if (output) {
                 cat("Reduction\n")
                 if(.Platform$OS.type == "windows") flush.console()
               }
               f.reduce <- apply(simplex[-min.index,],1,function(x){func(link(x))})
               if (any(f.reduce<best.ever[1])) {
                 best.ever <- c(f.reduce[order(f.reduce)[1]], simplex[-min.index,][order(f.reduce)[1],])
               }
               if ((simAnn) && (temp>0)) f.reduce.hot <- f.reduce + rlog(p)  # RV is ADDED!
               else f.reduce.hot <- f.reduce
               min.index <- min.sample(f.hot)
             }
           }
           
      # Check if simplex is degeneraed?
      doubles<-apply(simplex,2,duplicated)
      doublesdc<-apply(doubles,2,sum)
      if (any(doublesdc==p)){
       # If degenerated simplex replace by random vectors
        simplex[apply(doubles,1,any),apply(doubles,2,any)]<-rmunif(sum(apply(doubles,1,any)),minir[apply(doubles,2,any)],maxir[apply(doubles,2,any)])
        f <- apply(simplex,1,function(x){func(link(x))})  
        }
      close.enough <-((sd(f) <= e) || any(f <= best.possible))
      if (output) {
        if (simAnn) cat(paste(i,".",sep=""),"iteration; ",paste("temperature: ",as.character(signif(temp,4)),"; ", sep=""),"best/worst value:",best.ever[1],"/",max(f),"\n")
        else cat(paste(i,".",sep=""),"iteration; best/worst value:",f[min.index],"/",f[max.index],"\n")
        if(.Platform$OS.type == "windows") flush.console()
      }
      i <- i+1
      if (simAnn) { # adjust temperature 
        if (schedule==1) {
          temp <- T.start * faktor^i
          if (temp < zero.temp) temp <- 0
        }
        else {
          if (i<K) temp <- T.start * (1-i/K)^alpha
          else temp <- 0
        }
        if (temp == 0) simAnn <- FALSE
      }
    }
    if (output) {
      if (close.enough) if (any(f <= best.possible)) cat("Best possible value reached.\n")
                        else cat("Converged.\n")
      else cat("Stopped after",i-1,"iterations.\n")
    }
    if ((close.enough) | any(f <= best.possible)) converged <- TRUE
    else converged <- FALSE
    simplex<-t(apply(simplex,1,link))
    result <- list(minimum=link(best.ever[-1]), value=best.ever[1], iter=i-1, converged=converged, 
                   epsilon=sd(f), startsimplex=startsimplex, finalsimplex=simplex)
    return(result)
  }
  
  crossval.sample <- function(grouping, fold=10)
  # returns more or less equally sized cross-validation-samples. 
  {
    grouping <- factor(grouping)
    g <- length(levels(grouping)) #number of groups
    if (fold > length(grouping)) fold <- length(grouping)
    cv.groups <- rep(0,length(grouping))
    groupsizes <- c(0,cumsum(summary(grouping)))
    numbers <- c(rep(1:fold, length(grouping) %/% fold), 
                 sample(fold, length(grouping) %% fold)) # group numbers to be assigned 
    for (lev in 1:g) {
      index <- which(grouping==factor(levels(grouping)[lev], levels=levels(grouping))) # indices of class "lev"
      cv.groups[index] <- sample(numbers[(groupsizes[lev]+1):groupsizes[lev+1]])
    }
    return(cv.groups)
  }

#                                                      #
#  Beginning of  _M_A_I_N_ _P_R_O_C_E_D_U_R_E_  (RDA)  #
#                                                      #
  data <- x
  rm(x)
  # if `grouping' vector not given, first data column is taken. 
  if (is.null(grouping)) {
    grouping <- data[,1]
    data <- data[,-1]
  }
  data <- as.matrix(data)
  stopifnot(dim(data)[1]==length(grouping))
  grouping <- factor(grouping)
  classes <- levels(grouping)
  if (!is.null(dimnames(data))) varnames <- dimnames(data)[[2]]
  else varnames <- NULL
  grouping <- as.integer(grouping)
  if (is.null(prior)) { 
    prior <- tabulate(grouping)
    prior <- prior / sum(prior) 
  } # frequencies as prior
  else if (all(prior == 1)) 
    prior <- rep(1 / length(classes), length(classes))   # uniform prior
  names(prior) <- classes
  dimnames(data) <- NULL
  g <- max(grouping)             # number of groups 
  p <- ncol(data)                # number of variables 
  n <- length(grouping)          # number of observations 
  if (all(is.finite(regularization)))  # no optimization 
    opti <- list(minimum=regularization, conv=FALSE, iter=0) 
  else { # optimization 
    bothpar <- (!any(is.finite(regularization)))
    n.i <- rep(round(train.fraction*n), fold) # number of observations in bootstrap (training-)samples 
    if (output) {
      cat(" - RDA -\n")
      cat(n, "observations of", p, "variables in", g, "classes,\n")
      if (crossval) cat(paste(fold, "-fold cross-validation.", sep=""),"\n")
      else cat(fold, "bootstrap samples of", n.i[1], "observations each.\n")
      cat("Class names: ", paste(classes[1:(length(classes)-1)], col=",", sep=""),
          classes[length(classes)], "\n")
      if(.Platform$OS.type == "windows") flush.console()
    }  
    # draw bootstrap/crossval samples (row indices): 
    train <- NULL
    test.freq <- NULL
      if (crossval) { # cross-validation 
        indi <- crossval.sample(grouping, fold)
        tabu <- length(indi)-tabulate(indi)
        train <- matrix(0, nrow=max(tabu), ncol=fold)
        for (i in 1:fold) {
          train[1:tabu[i],i] <- which(indi != i)
          test.freq <- cbind(test.freq, tabulate(grouping[-train[,i]], g))
        }
        n.i <- apply(train, 2, function(x) sum(x > 0))
      }
    else { # no cross-validation, but bootstrapping 
      for (i in 1:fold) {
        new <- NULL
        for (j in 1:g) new <- c(new, sample(which(grouping == j), 2))
        new <- c(new, sample((1:n)[-new], n.i - 2 * g))
        train <- cbind(train, new)
        test.freq <- cbind(test.freq, tabulate(grouping[-train[,i]], g))
        # each sample now contains at least 2 elements from each group. 
        # (samples = columns) 
      }
      remove("new")
    }
    # compute parameter estimates (mu & Sigma) for each group in each training sample: 
    means <- covariances <- covpooled <- NULL
    for (i in 1:fold) {
      means <- array(c(means, 
                       as.vector(t(as.matrix(aggregate(data[train[,i],], 
                            by = list(grouping[train[,i]]), mean)[,-1])))),
                     c(p, g, i))
      # (p x g x i)-Array ... each "slice" contains g group means as column vectors. 
      new.covar <- array(unlist(by(data[train[,i],], grouping[train[,i]], var)), c(p,p,g))
      # (p x p x g)-Array, each slice is covariance matrix of one group. 
      covariances <- array(c(covariances,new.covar), c(p,p,g,i))
      # (p x p x g x i)-Array, each "hyperslice" contains g covariance matrices, as above. 
      new.cp <- array(new.covar, c(p*p,g))
      weights <- (tabulate(grouping[train[,i]])-1) / (n.i[i]-g) 
      # weights proportional to fraction of group in sample 
      covpooled <- array(c(covpooled, matrix(new.cp %*% weights, p, p)), c(p, p, i))
      # (p x p x i)-Array, each slice contains pooled covariance for i-th bootstrap sample. 
    }
    remove(list=c("new.covar","new.cp","weights"))
    if (bothpar) { # optimization over both parameters 
      if (is.null(startsimplex)) {
        #startsimplex <- matrix(rbeta(6,1,1),ncol=2) # Beta-RVs for Startsimplex
        startsimplex <- cbind(c(runif(1,1/11,4/11),runif(1,4/11,7/11),runif(1,7/11,10/11)),
                              c(runif(1,1/11,4/11),runif(1,4/11,7/11),runif(1,7/11,10/11)))
        perm <- cbind(c(1,3,2), c(2,1,3), c(3,1,2), c(2,3,1))
        startsimplex[,2] <- startsimplex[perm[,sample(4,1)],2]
      }
      if (trafo) {  # use transformation in Nelder-Mead. 
        linkfunc <- function(x)
            return(1/(1+exp(-x))) # sigmoidal transformation function (for both parameters) 
        linkinverse <- function(x)
            return(-log(1/x-1))   # inverse of link function 
        #startsimplex <- matrix(rnorm(6,0,1),ncol=2) # Normal-RVs for Startsimplex 
        startsimplex <- linkinverse(startsimplex) # transform startsimplex
        mini <- c(-Inf, -Inf)
        maxi <- c(Inf, Inf)
      }
      else {        # do not use transformation. 
        linkfunc <- function(x)
            return(x) # identity function 
        #linkinverse <- linkfunc
        mini=c(0,0)
        maxi=c(1,1)
      }
      dimnames(startsimplex) <- list(NULL, c("gamma","lambda"))
      # NELDER-MEAD #
      opti <- nelder.mead(goalfunc, startsimplex, mini=mini, maxi=maxi, 
        link=linkfunc, max.iter=max.iter, best.possible=0, out=output, 
        simAnn=simAnn, schedule=schedule, T.start=T.start, 
        halflife=halflife, zero.temp=zero.temp, alpha=alpha, K=K)
    }
    else { # optimization over single parameter 
      logit <- function(x)  return(1/(1+exp(-x)))
      tryval <- logit(seq(-4,4,le=12)) # values to try first 
      if (is.na(regularization[1])) {
        if (output) {
          cat("Optimizing gamma...\n")
          if(.Platform$OS.type == "windows") flush.console()
        }
        goalfu2 <- function(x)
            return(goalfunc(c(x, regularization[2])))
        err <- apply(matrix(tryval, ncol=1), 1, goalfu2)
        fromto <- c(0, tryval, 1)[which.min(err) + c(0,2)]
        minimize <- optimize(goalfu2, fromto)
        opti <- list(minimum=c(minimize$minimum, regularization[2]), 
                     value=minimize$objective, conv=TRUE, iter=-1)
      }
      else {
        if (output) {
          cat("Optimizing lambda...\n")
          if(.Platform$OS.type == "windows") flush.console()
        }
        goalfu2 <- function(x)
        {return(goalfunc(c(regularization[1], x)))}
        err <- apply(matrix(tryval,ncol=1),1,goalfu2)
        fromto <- c(0,tryval,1)[which.min(err)+c(0,2)]
        minimize <- optimize(goalfu2,fromto)
        opti <- list(minimum=c(regularization[1], minimize$minimum), 
                     value=minimize$objective, conv=TRUE, iter=-2)     
      }
    }
  }
  opt.par <- opti$minimum; names(opt.par) <- c("gamma","lambda")
  if (output) {
    cat("Regularization parameters:\n gamma:", round(opt.par[1],5), 
        "  lambda:", round(opt.par[2],5), "\n")
    if(.Platform$OS.type == "windows") flush.console()
  }
  # compute parameters for complete data: 
  means <- t(as.matrix(aggregate(data,by=list(grouping),mean)[,-1]))
  dimnames(means) <- list(varnames,classes)
 # covariances <- array(unlist(by(data[train[,i],],grouping[train[,i]],var)),
 #                      c(p,p,g), dimnames=list(varnames,varnames,classes))
  covariances <- array(unlist(by(data,grouping,var)),
                       c(p,p,g), dimnames=list(varnames,varnames,classes))
 # weights <- (tabulate(grouping)-1) / (n.i-g)
  weights <- (tabulate(grouping)-1) / (n-g)
  covpooled <- matrix(array(covariances, c(p*p,g)) %*% weights, p, p,
                      dimnames = list(varnames, varnames))
  # predict training data, compute apparent error rate: 
  if (estimate.error) {
    errors <- classify(data, means, covariances, covpooled, opt.par[1], opt.par[2], g=g, p=p, n=n, prior=prior) != grouping
    group.rates <- tabulate(grouping[errors],g) / tabulate(grouping)
    APER <- as.vector(t(prior) %*% group.rates)
    if (output) 
        cat("Apparent error rate (APER) for training data:", 
            round(APER * 100, 3), "%\n")  
  }
  else APER <- NA
  if (crossval) err <- c("APER"=APER, "crossval"=opti$value)
  else err <- c("APER"=APER, "bootstrap"=opti$value)
  result <- list(call=match.call(), 
                 regularization=opt.par, classes=classes, prior=prior, error.rate=err,
                 varnames=varnames,
                 means=means, covariances=covariances, covpooled=covpooled, 
                 converged=opti$conv, iter=opti$iter)
  class(result) <- "rda"
  return(result)
}

predict.rda <- function(object, newdata, posterior=TRUE, aslist=TRUE, ...)
{
  classify <- function(dataset, mu, sigma, pooled, gamma, lambda, g, p, n, prior)
  # difference to `classify'-function above is that LIKELIHOODS are returned instead of classifications
  # mu            : matrix with group means as columns 
  # sigma         : (p x p x g)-Array of group covariances 
  # pooled        : pooled covariance matrix 
  # gamma, lambda : regularization parameters
  # g, p, n       : numbers of classes, variables & observations
  {
    # compute likelihood of each datum for each group: 
    likelihood <- matrix(nrow=n, ncol=g)
    for (i in 1:g){ 
      reg.cov <- (1-lambda) * sigma[,,i] + lambda * pooled               # shift towards pooled cov. 
      reg.cov <- if (is.matrix(reg.cov)) # cov is a matrix
                   (1-gamma)*reg.cov + gamma*mean(diag(reg.cov))*diag(p) # shift towards identity 
                 else # one variable  =>  cov is 1x1-matrix:
                   (1-gamma)*reg.cov + gamma*reg.cov                     # shift towards identity 
      likelihood[,i] <- prior[i] * .dmvnorm(dataset, mu = mu[,i], inv.sigma = solve(reg.cov))
    }
    return(likelihood)
  }

  p <- dim(object$means)[1]
  g <- dim(object$means)[2]
  
      if (!inherits(object, "rda")) 
        stop("object not of class rda")
        if (!is.null(Terms <- object$terms)) {
        if (missing(newdata)) 
            newdata <- model.frame(object)
            else {
                newdata <- model.frame(as.formula(delete.response(Terms)), 
                newdata, na.action = function(x) x, xlev = object$xlevels)
            }
            x <- model.matrix(delete.response(Terms), newdata, contrasts = object$contrasts)
            xint <- match("(Intercept)", colnames(x), nomatch = 0)
            if (xint > 0) 
                x <- x[, -xint, drop = FALSE]
            }
        else {
        if (missing(newdata)) {
             if (!is.null(sub <- object$call$subset)) 
                    newdata <- eval.parent(parse(text = paste(deparse(object$call$x, 
                      backtick = TRUE), "[", deparse(sub, backtick = TRUE), 
                    ",]")))
                else newdata <- eval.parent(object$call$x)
                if (!is.null(nas <- object$call$na.action)) 
                 newdata <- eval(call(nas, newdata))
            }
            if (is.null(dim(newdata))) 
             dim(newdata) <- c(1, length(newdata))
            x <- as.matrix(newdata)
        }

  
  
  #if (!any(is.null(colnames(newdata)),is.null(object$varnames))) { # (both colnames & varnames are given) 
  #  if(all(is.element(object$varnames,colnames(newdata)))){        # (varnames is a subset of colnames)   
  #    newdata <- as.matrix(newdata[,object$varnames])
  #  }
  #}
  #if(is.vector(newdata)) newdata <- matrix(newdata, ncol=1)
  n <- dim(newdata)[1]
  newdata<-x
  likeli <- classify(newdata, object$means, object$covariances, object$covpooled,
                     object$regul[1], object$regul[2], g=g, p=p, n=dim(newdata)[1], 
                     prior=object$prior)
  colnames(likeli) <- object$classes
  classi <- apply(likeli, 1, function(x) order(x)[g])
  classi <- factor(object$classes[classi], levels = object$classes)
  postmat <- if (posterior) likeli / rowSums(likeli)
             else NULL
  if (aslist) result <- list("class" = classi, "posterior" = postmat)
  else { 
    result <- classi
    attr(result, "posterior") <- postmat
  }
  return(result)
}


print.rda <- function(x,...)
{
  #cat(" - RDA - \n")
  cat("Call:\n")
  print(x$call)
  cat("\nRegularization parameters:\n")
  #cat("gamma:", round(x$regu[1],5), " lambda:", round(x$regu[2],5), "\n")
  print(x$regu)
  #cat("\nClass prior:\n")
  cat("\nPrior probabilities of groups:\n")
  print(x$prior)
  cat("\nMisclassification rate:\n")
  cat("       apparent:",
    ifelse(is.na(x$error.rate[1]), "--", as.character(round(x$error.rate[1] * 100, 3))), "%\n")
  if (length(x$error.rate) > 1)
    cat(ifelse(names(x$error.rate)[2] == "crossval", 
        "cross-validated:", "   bootstrapped:"), 
        as.character(round(x$error.rate[2] * 100, 3)), "%\n")
  invisible(x)
}


plot.rda <- function(x, textplot=FALSE, ...)
{
  parpty <- par("pty")
  par(pty="s")
  if(textplot) {
    plot(c(0,1),c(0,1), type="n", axes=FALSE, xlab="groups", ylab="covariances",...)
    textcol <- "darkgrey"
    textsize <- 1.5
    textshift <- 0.02
    text(0+textshift, 0+textshift, "QDA", adj=c(0,0), cex=textsize, col=textcol)
    text(1-textshift, 0+textshift, "LDA", adj=c(1,0), cex=textsize, col=textcol)
    #text(0+textshift, 1-textshift, "cond. indep.", adj=c(0,1), cex=textsize, col=textcol)
    #text(1-textshift, 1-textshift, "nearest mean", adj=c(1,1), cex=textsize, col=textcol)
    text(0.5, 1-textshift, "i.i.d. variables", adj=c(0.5,1), cex=textsize, col=textcol)
    axis(1, at=c(0,1), labels=c("unequal", "equal"))
    axis(2, at=c(0,1), labels=c("correlated", "diagonal"))
  }
  else{
    plot(c(0,1),c(0,1), type="n", axes=FALSE, xlab=expression(lambda), ylab=expression(gamma),...)
    axis(1)
    axis(2)
  }
  lines(c(0,1,1,0,0), c(0,0,1,1,0), col="grey")
  lines(c(0,1), rep(x$regu[1],2), col="red1", lty="dotted")
  lines(rep(x$regu[2],2), c(0,1), col="red1", lty="dotted")
  points(x$regu[2], x$regu[1], pch=18, col="red2")
  par(pty=parpty)
  invisible(x$regu)
}


.dmvnorm <- function(x, mu=NA, inv.sigma=NA, logDetCov=NA)    
# Density of a Multivariate Normal Distribution 
# works for matrices as well as for vectors     
# (matrices are evaluated row-wise)             
#   !!   supply inverse of covariance   !!      
# `logDetCov' = log(det(Cov)) = log(1/det(inv.Cov)) = -1*log(det(inv.Cov))
{
  if (is.vector(x)) x <- t(x) # x is treated as 1 obsevation of (length(x)) variables
  if (is.na(logDetCov)) logDetCov <- -log(det(inv.sigma))
  singledens <- function(x, M=mu, IS=inv.sigma, ldc=logDetCov)  # density function for a single vector 
  {
    xm <- x - M
    return(as.numeric(exp(- 0.5 * 
           (length(M) * log(2*pi) 
            + ldc
            + (t(xm) %*% IS %*% xm)))))
  } 
  return(apply(x, 1, singledens))#, M=mu, IS=inv.sigma, ldc=logDetCov))
}  
shardsplot <- function(object, plot.type = c("eight", "four", "points", "n"),
    expand = 1, stck = TRUE, grd = FALSE, standardize = FALSE, data.or = NA,
    label = FALSE, plot = TRUE, classes = 0, vertices = TRUE,
    classcolors = "rainbow", wghts = 0, xlab = "Dimension 1", ylab = "Dimension 2", xaxs = "i", yaxs = "i", ...){
    diss <- FALSE
    plot.type <- match.arg(plot.type) # "delaunay" not yet implemented
    if (plot.type=="delaunay") grd <- FALSE
    if (class(object)=="som"){
        preimages <- object$code
        nzx <- object$x
        nzy <- object$y
        Z0 <- t(matrix(1:(nzx*nzy),nzx,nzy))
        nobs <- object$code.sum$nobs
        nobs.col <- nobs+1
        maxn <- max(nobs.col)
        cl.ord <- nobs.col
    }
    else if (class(object)=="EDAM"){
        preimages <- object$preimages
        Z0 <- object$Z
        if (classes[1]==0) cl.ord <- as.numeric(object$cl.ord)
        if (classes[1]!=0) cl.ord <- as.numeric(classes[as.vector(t(object$Z.old.terms))])
        if (ncol(preimages)==nrow(preimages) && preimages==t(preimages)){
            diss <- TRUE
            EV.dist <- preimages
            preimages <- cbind(c(1:nrow(preimages)),c(1:nrow(preimages)))
            rownames(preimages) <- rownames(object$preimages)
        }
    }
    if (standardize){
        EV.sd <- sqrt(apply(preimages,2,var))
        EV.sd[!EV.sd] <- 1
        preimages <- t(t(preimages)/EV.sd)
    }
    if (wghts[1]!=0) preimages <- preimages*kronecker(t(wghts),rep(1,nrow(preimages)))
    nzx <- ncol(Z0)
    nzy <- nrow(Z0)
    nzx.ex <- round(nzx*expand)
    nzy.ex <- round(nzy*expand)
    Cells0 <- cbind(kronecker(c(1:nzy),rep(1,nzx)),rep(c(1:nzx),nzy))
    Cells.ex <- Cells0*expand
    Cells.exldru <- Cells.ex
    Cells.exlurd <- Cells.ex
    EV <- preimages

    if (!diss) EV.dist <- distmirr(dist(EV, method = "euclidean"))
    dist.max <- 0
    dist.min <- max(EV.dist)
    for (i in 1:(nzy*nzx)){
        sub.mat <-  t(matrix(rep(Cells0[i,], nzx*nzy), 2, nzx*nzy))
        rel.dists <- EV.dist[which(diag((Cells0-sub.mat)%*%t(Cells0-sub.mat)) <= 2),i]
        rel.dists <- rel.dists[which(rel.dists > 0)]
        if (max(rel.dists) > dist.max) dist.max <- max(rel.dists)
        if (min(rel.dists) < dist.min) dist.min <- min(rel.dists)
    }
    Cells.ny <- (nzy*expand)
    Cells.nx <- (nzx*expand)
    if (any(expand!=1) || stck){
        for(j in 1:nzx){
            nb.dists <- rep(0, nzy-1)
            for (i in (2:nzy)){
              nb.dists[i-1] <- EV.dist[Z0[i,j],Z0[i-1,j]]
            }
            nb.dc <- sum(nb.dists)
            nb.pos <- c(expand,rep(0,nzy-1))
            for (l in 1:(nzy-1)){
                nb.pos[l+1] <- expand+(sum(nb.dists[1:l])*(Cells.ny-expand))/nb.dc
            }
            Cells.ex[Z0[,j],1] <- nb.pos
        }
        for (i in 1:nzy){
            nb.dists <- rep(0, nzx-1)
            for (j in (2:nzx)){
                nb.dists[j-1] <- EV.dist[Z0[i,j],Z0[i,j-1]]
            }
            nb.dc <- sum(nb.dists)
            nb.pos <- c(expand,rep(0, nzx-1))
            for (l in 1:(nzx-1)){
                nb.pos[l+1] <- expand+(sum(nb.dists[1:l])*(Cells.nx-expand))/nb.dc
            }
            Cells.ex[Z0[i,],2] <- nb.pos
        }
        if (min(nzx,nzy) > 2){
            for (j in 1:(nzy-2)){
                len <- min(nzx,nzy-j+1)
                nb.dists <- rep(0, len-1)
                rel.numbers <- rep(0,len)
                rel.numbers[1] <- Z0[j,1]
                for (i in (2:len)){
                    rel.numbers[i] <- Z0[i+j-1,i]
                    nb.dists[i-1] <- EV.dist[rel.numbers[i],rel.numbers[i-1]]
                }
                nb.dc <- sum(nb.dists)
                nb.pos <- cbind(c(j*expand,rep(0,len-1)),c(expand,rep(0,len-1)))
                for (l in 1 :(len-1)){
                    nb.pos[l+1,1] <- j*expand+(sum(nb.dists[1:l])*(len*expand-expand))/nb.dc
                    nb.pos[l+1,2] <- expand+(sum(nb.dists[1:l])*(len*expand-expand))/nb.dc
                }
                Cells.exldru[rel.numbers,] <- nb.pos
            }
            for (i in 2:(nzx-2)){
                len <- min(nzx-i+1,nzy)
                nb.dists <- rep(0, len-1)
                rel.numbers <- rep(0,len)
                rel.numbers[1] <- Z0[1,i]
                for (j in (2:len)){
                    rel.numbers[j] <- Z0[j,j+i-1]
                    nb.dists[j-1] <- EV.dist[rel.numbers[j],rel.numbers[j-1]]
                }
                nb.dc <- sum(nb.dists)
                nb.pos <- cbind(c(expand,rep(0,len-1)),c(i*expand,rep(0,len-1)))
                for (l in 1 :(len-1)){
                    nb.pos[l+1,1] <- expand+(sum(nb.dists[1:l])*(len*expand-expand))/nb.dc
                    nb.pos[l+1,2] <- i*expand+(sum(nb.dists[1:l])*(len*expand-expand))/nb.dc
                }
                Cells.exldru[rel.numbers,] <- nb.pos
            }
            for (j in 1:(nzy-2)){
                len <- min(nzx,nzy-j+1)
                nb.dists <- rep(0, len-1)
                rel.numbers <- rep(0,len)
                rel.numbers[1] <- Z0[j,nzx]
                for (i in (2:len)){
                    rel.numbers[i] <- Z0[i+j-1,nzx-i+1]
                    nb.dists[i-1] <- EV.dist[rel.numbers[i],rel.numbers[i-1]]
                }
                nb.dc <- sum(nb.dists)
                nb.pos <- cbind(c(j*expand,rep(0,len-1)),c(Cells.nx,rep(0,len-1)))
                for (l in 1 :(len-1)){
                    nb.pos[l+1,1] <- j*expand+(sum(nb.dists[1:l])*(len*expand-expand))/nb.dc
                    nb.pos[l+1,2] <- Cells.nx-(sum(nb.dists[1:l])*(len*expand-expand))/nb.dc
                }
                Cells.exlurd[rel.numbers,] <- nb.pos
            }
            for (i in 2:(nzx-1)){
                len <- min(i,nzy)
                nb.dists <- rep(0, len-1)
                rel.numbers <- rep(0,len)
                rel.numbers[1] <- Z0[1,i]
                for (j in (2:len)){
                    rel.numbers[j] <- Z0[j,i-j+1]
                    nb.dists[j-1] <- EV.dist[rel.numbers[j],rel.numbers[j-1]]
                }
                nb.dc <- sum(nb.dists)
                nb.pos <- cbind(c(expand,rep(0,len-1)),c(i*expand,rep(0,len-1)))
                for (l in 1:(len-1)){
                    nb.pos[l+1,1] <- expand+(sum(nb.dists[1:l])*(len*expand-expand))/nb.dc
                    nb.pos[l+1,2] <- i*expand-(sum(nb.dists[1:l])*(len*expand-expand))/nb.dc
                }
                Cells.exlurd[rel.numbers,] <- nb.pos
            }
        }
        if (plot.type!="four") Cells.ex <- (Cells.ex+Cells.exldru+Cells.exlurd)/3
        if (grd) Cells.ex <- round(Cells.ex)
    }
    if (plot){
        plot(c(expand,expand*nzy), c(expand,expand*nzx), type = "n", 
            xlab = xlab, ylab = ylab, xaxs = xaxs, yaxs = yaxs, ...)
        if (vertices && plot.type!="delaunay"){
            for (i in (1:nzx)){
                for (j in (1:nzy)){
                    if (j>1) lines(Cells.ex[Z0[c(j,j-1),c(i,i)],], col="gray")
                    if (i>1) lines(Cells.ex[Z0[c(j,j),c(i,i-1)],], col="gray")
                }
            }
        }
    }
    if (class(object)!="som"){
        vec.col.ord <-
            if(is.character(classcolors) && length(classcolors)==1)
                switch(classcolors,
                    "rainbow" = rainbow(max(cl.ord))[cl.ord],
                    "topo"    = topo.colors(max(cl.ord))[cl.ord],
                    "gray"    = gray(1:max(cl.ord)/max(cl.ord))[cl.ord],
                    stop("argument classcolors only support 'rainbow', 'topo', and 'gray'.")
                )
            else classcolors[cl.ord]
    }
    if (class(object)=="som"){
        if (any(is.na(data.or)))
            vec.col.ord <- c(rgb(1,1,1),
                hsv(1,1, ((maxn*1.5):1) / (maxn*1.5)))[nobs.col]
        if (!any(is.na(data.or))){
            code.classes <- rep(1,nzx*nzy)
            i <- 0
            nob <- nrow(data.or)
            dimen <- ncol(data.or)-1
            obj.codes <- object$visual$x+object$visual$y*nzx+1
            class.vec <- as.factor(data.or[,dimen+1])
            for (i in 1:(nzx*nzy)){
                class.table <- table(class.vec[(obj.codes==i)])
                if (length(class.table)>0) code.classes[i] <- which(class.table==max(class.table))[1]
            }
            vec.col.ord1 <-
                if(is.character(classcolors) && length(classcolors)==1)
                    switch(classcolors,
                        "rainbow" = rainbow(max(code.classes))[code.classes],
                        "topo"    = topo.colors(max(code.classes))[code.classes],
                        "gray"    = gray(1:max(code.classes)/max(code.classes))[code.classes],
                        stop("argument classcolors only support 'rainbow', 'topo', and 'gray'.")
                    )
                else classcolors[cl.ord]
            vec.col.ord <- rep(0,nzx*nzy)
            rgb.mat <- col2rgb(vec.col.ord1)/255
            for (i in 1:(nzx*nzy)){
                rgb.act <- rgb.mat[,i]
                V.act <- max(rgb.act)
                rgb.min <- min(rgb.act)
                S.act <- 0
                if (V.act>0) S.act <- (V.act-rgb.min)/V.act
                H.act <- 0
                if (S.act>0){
                    if (V.act==rgb.act[1]) H.act <- (rgb.act[2]-rgb.act[3])/(V.act-rgb.min)
                    else if (V.act==rgb.act[2]) H.act <- 2 + (rgb.act[3]-rgb.act[1])/(V.act-rgb.min)
                    else if (V.act==rgb.act[3]) H.act <- 4 + (rgb.act[1]-rgb.act[2])/(V.act-rgb.min)
                    H.act <- H.act*60
                    if (H.act<0) H.act <- H.act+360
                    H.act <- H.act/360
                }
                S.act <- (maxn*1.5-nobs[i])/(maxn*1.5)
                vec.col.ord[i] <- hsv(H.act, V.act, S.act)
            }
            vec.col.ord[!nobs] <- "#FFFFFF"
        }
    }
    for (i in 1:(nzy*nzx)){
        sub.mat <-  t(matrix(rep(Cells0[i,], nzx*nzy), 2, nzx*nzy))
        Cells0msubmat <- Cells0 - sub.mat
        rel.numbers <- which(diag(Cells0msubmat %*% t(Cells0msubmat)) <= 2)
        rel.dists <- EV.dist[rel.numbers, i]
        rel.numbers <- rel.numbers[rel.dists > 0]
        rel.dists <- rel.dists[rel.dists > 0]
        rel.points <- matrix(rep(Cells.ex[i,], length(rel.dists)),
            length(rel.dists), 2, byrow = TRUE)
        rel.points <- rel.points + Cells0msubmat[rel.numbers,] *
            (0.5-((rel.dists-dist.min) / (2*dist.max)))
        rel.coords <- Cells0msubmat[rel.numbers,]
        if(length(rel.dists) == 3){
            rel.points <- rbind(rel.points, Cells0[i,])
            rel.coords <- rbind(rel.coords, c(0,0))
            rel.numbers <- c(rel.numbers, i)
        }
        rel.count <- 0
        circle.mat <- matrix(c(-1, -1, -1, 0, -1, 1, 0, 1, 1, 1, 1, 0, 1, -1, 0, -1),
            8, 2, byrow = TRUE)
        draw.points <- matrix(rep(Cells.ex[i,], 8), 8, 2, byrow = TRUE)
        for(l in 1:8){
            for(j in 1:nrow(rel.coords)){
                if (t(rel.coords[j,]-circle.mat[l,]) %*% (rel.coords[j,]-circle.mat[l,]) == 0){
                    rel.count <- rel.count+1
                    draw.points[l,] <- rel.points[j,]
                }
            }
        }
        if (cl.ord[1] && class(object)!="som"){
            if (Cells0[i,1] > 1 && Cells0[i,2] > 1){
                if (cl.ord[Z0[Cells0[i,1]-1, Cells0[i,2]]] == cl.ord[i] &&
                    cl.ord[Z0[Cells0[i,1], Cells0[i,2]-1]] == cl.ord[i]){
                    if (cl.ord[Z0[Cells0[i,1]-1,Cells0[i,2]-1]] != cl.ord[i]){
                        draw.points[1,] <- c(draw.points[2,1], draw.points[8,2])
                    }
                }
            }
            if (Cells0[i,1] > 1 && Cells0[i,2] < nzx){
                if (cl.ord[Z0[Cells0[i,1]-1,Cells0[i,2]]] == cl.ord[i] &&
                    cl.ord[Z0[Cells0[i,1],Cells0[i,2]+1]] == cl.ord[i]){
                    if (cl.ord[Z0[Cells0[i,1]-1,Cells0[i,2]+1]] != cl.ord[i]){
                        draw.points[3,] <- c(draw.points[2,1], draw.points[4,2])
                    }
                }
            }
            if (Cells0[i,1] < nzy && Cells0[i,2] > 1){
                if (cl.ord[Z0[Cells0[i,1]+1,Cells0[i,2]]] == cl.ord[i] &&
                    cl.ord[Z0[Cells0[i,1],Cells0[i,2]-1]] == cl.ord[i]){
                    if (cl.ord[Z0[Cells0[i,1]+1,Cells0[i,2]-1]] != cl.ord[i]){
                        draw.points[7,] <- c(draw.points[6,1], draw.points[8,2])
                    }
                }
            }
            if (Cells0[i,1] < nzy && Cells0[i,2] < nzx){
                if (cl.ord[Z0[Cells0[i,1]+1, Cells0[i,2]]] == cl.ord[i] &&
                    cl.ord[Z0[Cells0[i,1],Cells0[i,2]+1]] == cl.ord[i]){
                        if (cl.ord[Z0[Cells0[i,1]+1, Cells0[i,2]+1]] != cl.ord[i]){
                            draw.points[5,] <- c(draw.points[6,1], draw.points[4,2])
                    }
                }
            }
        }
        if (plot.type == "four"){
            draw.points[1,] <- c(draw.points[2,1], draw.points[8,2])
            draw.points[3,] <- c(draw.points[2,1], draw.points[4,2])
            draw.points[7,] <- c(draw.points[6,1], draw.points[8,2])
            draw.points[5,] <- c(draw.points[6,1], draw.points[4,2])
        }
        draw.points[(draw.points[,2] < 1),2] <- expand
        draw.points[(draw.points[,1] < 1),1] <- expand
        if (plot && plot.type!="delaunay" && plot.type!="n"){
            if (plot.type!="points") {
                polygon(draw.points, col = vec.col.ord[i], ...)
                if (label)
                    text(mean(draw.points[,1]), mean(draw.points[,2]),
                        rownames(preimages)[i], ...)
            }
            if (plot.type=="points")
                points(Cells.ex[i,1], Cells.ex[i,2], col=vec.col.ord[i], pch = 19, ...)
        }
        if (label && plot.type!="delaunay")
            text(mean(draw.points[,1]), mean(draw.points[,2]), rownames(preimages)[i], ...)
    }
    if (plot && plot.type=="delaunay"){
        lnames <- 0
        if (label) lnames <- rownames(preimages)
        cont.shardsplot(coords=Cells.ex, classes=cl.ord, ydistmat=EV.dist,
            radius=radius, smallest=smallest, percentage=percentage, convexify=convexify,
            label=label, vertices=vertices, classcolors=classcolors, pnew=FALSE,
            vec.col.ord=vec.col.ord, lnames=lnames, ...)
    }
    Cells.ex.dist <- distmirr(dist(Cells.ex, method = "euclidean"))
    results <- list(Cells.ex = Cells.ex, S = as.numeric(TopoS(EV.dist,Cells.ex.dist)))
    class(results) <- "EDAM.ex"
    invisible(results)
}
sknn<-function (x, ...) 
    UseMethod("sknn")

sknn.default <- function(x,grouping,k=3,...)
{
cl <- match.call()
cl[[1]] <- as.name("sknn")
structure(list(learn = x, grouping = grouping, lev = levels(grouping), k=k, call = cl), class = "sknn")
}


### sknn bei verschiedenen Eingabeformaten:
sknn.formula<-function (formula, data = NULL, ..., subset, na.action = na.fail) 
{
    m <- match.call(expand.dots = FALSE)
    if (is.matrix(eval.parent(m$data))) 
        m$data <- as.data.frame(data)
    m$... <- NULL
    m[[1]] <- as.name("model.frame")
    m <- eval.parent(m)
    Terms <- attr(m, "terms")
    grouping <- model.response(m)
    x <- model.matrix(Terms, m)
    xvars <- as.character(attr(Terms, "variables"))[-1]
    if ((yvar <- attr(Terms, "response")) > 0) 
        xvars <- xvars[-yvar]
    xlev <- if (length(xvars) > 0) {
        xlev <- lapply(m[xvars], levels)
        xlev[!sapply(xlev, is.null)]
    }
    xint <- match("(Intercept)", colnames(x), nomatch = 0)
    if (xint > 0) 
        x <- x[, -xint, drop = FALSE]
    res <- sknn.default(x, grouping, ...)
    res$terms <- Terms
    cl <- match.call()
    cl[[1]] <- as.name("sknn")
    res$call <- cl
    res$contrasts <- attr(x, "contrasts")
    res$xlevels <- xlev
    attr(res, "na.message") <- attr(m, "na.message")
    if (!is.null(attr(m, "na.action"))) 
        res$na.action <- attr(m, "na.action")
    res
}

sknn.matrix<-function (x, grouping, ..., subset, na.action = na.fail) 
{
    if (!missing(subset)) {
        x <- x[subset, , drop = FALSE]
        grouping <- grouping[subset]
    }
    if (!missing(na.action)) {
        dfr <- na.action(structure(list(g = grouping, x = x), 
            class = "data.frame"))
        grouping <- dfr$g
        x <- dfr$x
    }
    res <- sknn(x, grouping, ...)
    cl <- match.call()
    cl[[1]] <- as.name("sknn")
    res$call <- cl
    res
}

sknn.data.frame<-function (x, ...) 
{
   res <- sknn.matrix(structure(data.matrix(x), class = "matrix"), 
        ...)
    cl <- match.call()
    cl[[1]] <- as.name("sknn")
    res$call <- cl
    res
}


predict.sknn<-function(object, newdata,...)
{
spsknn <- function(x,object)
    {
        abstand<-apply(object$learn, 1, function(y) sum((y-x)^2))
        kdach<-(object$grouping[order(abstand)][1:object$k])
        return(table(kdach))
    } 

if (!inherits(object, "sknn")) 
        stop("object not of class sknn")
    if (!is.null(Terms <- object$terms)) {
        if (missing(newdata)) 
            newdata <- model.frame(object)
        else {
            newdata <- model.frame(as.formula(delete.response(Terms)), 
                newdata, na.action = function(x) x, xlev = object$xlevels)
        }
        x <- model.matrix(delete.response(Terms), newdata, contrasts = object$contrasts)
        xint <- match("(Intercept)", colnames(x), nomatch = 0)
        if (xint > 0) 
            x <- x[, -xint, drop = FALSE]
    }
    else {
        if (missing(newdata)) {
            if (!is.null(sub <- object$call$subset)) 
                newdataa <- eval.parent(parse(text = paste(deparse(object$call$x, 
                  backtick = TRUE), "[", deparse(sub, backtick = TRUE), 
                  ",]")))
            else newdata <- eval.parent(object$call$x)
            if (!is.null(nas <- object$call$na.action)) 
                newdata <- eval(call(nas, newdata))
        }
        if (is.null(dim(newdata))) 
            dim(newdata) <- c(1, length(newdata))
        x <- as.matrix(newdata)
    }

werte<-t(apply(x,1,spsknn,object=object))
werte<-werte/object$k
classes <- factor(max.col(werte), levels = seq(along = object$lev), 
        labels = object$lev)
result<-list(posterior=werte,class=classes)
return(result)
}
stepclass <- function(x, ...)
{
  UseMethod("stepclass")
}


stepclass.formula <- function(formula, data, method, ...)
{
  variables <- dimnames(attributes(terms(formula))$factors)[[1]]
  response <- variables[1]
  discriminators <- variables[-1]
  if(any(discriminators == ".")) {
    exclude <- c(response, discriminators[discriminators != "."])
    discriminators <- colnames(data)[!is.element(colnames(data), exclude)]
  }
  result <- stepclass(x=data[, discriminators], grouping=data[, response], 
                      method=method, ...)
  result$call <- match.call()
  return(result)
}
 

stepclass.default <-function (x, grouping, method, improvement = 0.05, 
    maxvar = Inf, start.vars = NULL, direction = c("both", "forward", "backward"), 
    criterion = "CR", fold = 10, cv.groups = NULL, output = TRUE, ...) 
{
    cr <- c("correctness rate", "accuracy", "abiltity to seperate", "confidence")
    switch(criterion,
        CR = cr <- cr[1],
        AC = cr <- cr[2], 
        AS = cr <- cr[3],
        CF = cr <- cr[4],
        {   criterion <- "CR"
            cr <- cr[1]
            cat("Unknown criterion. Changed to", cr, "\n")
        }
    )
    textoutput <- function(rate, variables, into.model = NA, 
        out.of.model = NA, variablenames = varnames) {
        if (is.na(into.model)) {
            if (!is.na(out.of.model)) 
                in.out <- paste("  out: \"", variablenames[out.of.model], 
                  "\"; ", sep = "")
            else in.out <- "  starting"
        }
        else in.out <- paste("  in: \"", variablenames[into.model], 
            "\"; ", sep = "")
        cat(paste(cr, ": ", 1-round(rate, 5), ";", in.out,
            " ", sep = ""))
        if (length(variables) == 1) 
            cat(paste("variables (1):", variablenames[variables], "\n"))
        else cat(paste("variables (", length(variables), "):", sep = ""), 
            paste(variablenames[variables[-length(variables)]], 
            coll = ",", sep = ""), variablenames[variables[length(variables)]], "\n")
        invisible()
    }

    null.rate <- function(grouping, fold, cv.groups, criterion)
    {
        goalfunc<-0
        for (f in 1:fold) {
            train <- (cv.groups != f)
            test <- which(!train)
            train <- which(train)
            training<-grouping[train]
            hg <- table(training)
            classi <- matrix(rep(hg/length(training),length(test)), 
                nrow = length(test), ncol = length(levels(grouping)), byrow = TRUE)
            goalfunc <- goalfunc + (1 - ucpm(classi, grouping[test])[[criterion]])
        }
        return(goalfunc / fold)
    }
        
    
    
    cv.rate <- function(vars, data = data, grouping = grouping, 
        method = method, fold = fold, cv.groups = cv.groups, 
        criterion = criterion, ...) {
        predicted <- numeric(length(grouping))
 
        goalfunc <- 0
        for (f in 1:fold) {
            train <- (cv.groups != f)
            test <- which(!train)
            train <- which(train)
            traindat <- data[train, vars]
            if (is.vector(traindat)) 
                if (is.data.frame(data)) {
                  traindat <- as.data.frame(matrix(traindat, 
                    ncol = length(vars)))
                  names(traindat) <- names(data)[vars]
                }
                else {
                  traindat <- matrix(traindat, ncol = length(vars))
                }
            object <- try(do.call(method, list(traindat, grouping[train], ...)), 
                silent = TRUE)
            if (class(object) != "try-error") {
                testdat <- data[test, vars]
                if (is.vector(testdat)) 
                  if (is.data.frame(data)) {
                    testdat <- as.data.frame(matrix(testdat, 
                      ncol = length(vars)))
                    names(testdat) <- names(data)[vars]
                  }
                  else {
                    testdat <- matrix(testdat, ncol = length(vars))
                  }
                classi <- try(predict(object, testdat), silent = TRUE)
                if (class(classi) != "try-error") {
                  if (is.list(classi)) 
                    classi <- classi$posterior
                  predicted[test] <- as.numeric(max.col(classi))
                  goalfunc <- goalfunc + (1 - ucpm(classi, grouping[test])[[criterion]])
                }
     
            }
        #goalfunc<-goalfunc+(1-ucpm(classi,)$criterion)
        }
        #group.rates <- NULL
        #for (lev in 1:length(levels(grouping))) {
        #    group.rates <- c(group.rates, (1-ucpm(classi,)$criterion))
        #}
        
        if (any(predicted == 0)){
            goalfunc <- fold
            warning(" Error(s) in modeling/prediction step.\n")
        }
        return(goalfunc / fold)
    }
    min.sec <- function(seconds) {
        seconds <- if (length(seconds) >= 3) 
            seconds[3]
        else seconds[1]
        result <- c(seconds%/%3600)
        result <- c(hr = result, min = (seconds - result * 3600) %/% 60, 
            sec = (seconds - result * 3600) %% 60)
        return(result)
    }
    data <- x
    rm("x")
    switch(method, 
        lda = require("MASS"), 
        qda = require("MASS"), 
        rpart = require("rpart"), 
        naiveBayes = require("e1071"))
    stopifnot(dim(as.data.frame(data))[1] == length(grouping))
    runtime <- proc.time()[3]
    direction <- match.arg(direction)
    grouping <- factor(grouping)
    g <- length(levels(grouping))
    if (is.finite(maxvar)) 
        improvement <- 0
    fwd <- (direction == "forward") || (direction == "both")
    bwd <- (direction == "backward") || (direction == "both")
    if (!is.null(dimnames(data))) 
        varnames <- dimnames(data)[[2]]
    else varnames <- paste("var", as.character(1:dim(data)[2]), 
        sep = ".")
    if ((direction == "backward") && is.null(start.vars)) 
        start.vars <- 1:dim(data)[2]
    if (is.character(start.vars)) {
        model <- seq(along = varnames)[is.element(varnames, start.vars)]
        start.vars <- model
    }
    else model <- start.vars
    out <- 1:length(varnames)
    out <- out[!is.element(out, model)]
    finished <- FALSE
    if (is.null(cv.groups)) {
        if (fold > length(grouping)) 
            fold <- length(grouping)
        cv.groups <- rep(0, dim(data)[1])
        groupsizes <- c(0, cumsum(summary(grouping)))
        numbers <- c(rep(1:fold, length(grouping) %/% fold), 
            sample(fold, length(grouping) %% fold))
        for (lev in 1:g) {
            index <- which(grouping == factor(levels(grouping)[lev], 
                levels = levels(grouping)))
            cv.groups[index] <- sample(numbers[(groupsizes[lev] + 1):groupsizes[lev + 1]])
        }
    }
    else {
        cv.groups <- as.numeric(factor(cv.groups[1:length(grouping)]))
        fold <- max(cv.groups)
    }
    if (output) {
        cat(paste(" `stepwise classification', using ", fold, 
            "-fold cross-validated ", cr, " of method ", method,
            "'.\n", sep = ""))
        cat(dim(data)[1], "observations of", dim(data)[2], "variables in", 
            g, "classes; direction:", direction, "\n")
        if (!is.finite(maxvar)) 
            cat("stop criterion: improvement less than", 
                round(improvement * 100, 2), "%.\n")
        else cat("stop criterion: assemble", maxvar, "best variables.\n")
        if (.Platform$OS.type == "windows") 
            flush.console()
    }
    if (is.null(start.vars)){
        old.rate <- null.rate(grouping = grouping, fold = fold, 
            cv.groups = cv.groups, criterion = criterion)   
    }
    else {
        old.rate <- cv.rate(vars = start.vars, data = data, grouping = grouping, 
            method = method, fold = fold, cv.groups = cv.groups, 
            criterion = criterion, ...)
        if (output) {
            textoutput(old.rate, model)
            if (.Platform$OS.type == "windows") 
                flush.console()
        }
    }
    result.e <- old.rate
    result.v <- matrix(c("start", "0"), ncol = 2)
    last.changed <- NA
    while (!finished) {
        error.rates <- NULL
        if (fwd && (length(out[!is.element(out, last.changed)]) >= 1)) 
            for (tryvar in out[!is.element(out, last.changed)]) {
                newrate <- cv.rate(vars = c(model, tryvar), data = data, 
                  grouping = grouping, method = method, fold = fold, 
                  cv.groups = cv.groups, criterion = criterion, ...)
                error.rates <- rbind(error.rates, c(var = tryvar, rate = newrate))
            }
        else if (fwd && (length(out[!is.element(out, last.changed)]) == 0))
            error.rates <- rbind(error.rates, c(var = 1, rate = 1))
        if (bwd && (length(model[!is.element(model, last.changed)]) >= 2)) 
            for (outvar in model[!is.element(model, last.changed)]) {
                trymodel <- model[!is.element(model, outvar)]
                newrate <- cv.rate(trymodel, data = data, grouping = grouping, 
                  method = method, fold = fold, cv.groups = cv.groups, 
                  criterion = criterion, ...)
                error.rates <- rbind(error.rates, c(var = outvar, rate = newrate))
            }
        else if (bwd && (length(model[!is.element(model, last.changed)]) == 1)) 
            error.rates <- rbind(error.rates, c(var = 1, rate = 1))
        best <- order(error.rates[, 2])[1]
        last.changed <- error.rates[best, 1]
        #if (error.rates[best, 2] >= improvement * old.rate) {
        if ((old.rate - error.rates[best, 2]) < improvement) {
            if (error.rates[best, 2] >= old.rate) 
                finished <- TRUE
            else {
                if (is.element(error.rates[best, 1], model)) {
                  index <- is.element(model, error.rates[best, 1])
                  out <- c(out, model[index])
                  model <- model[!index]
                  result.v <- rbind(result.v, c("out", out[length(out)]))
                  result.e <- c(result.e, error.rates[best, 2])
                  if (output) {
                    textoutput(error.rates[best, 2], model, out.of.model = error.rates[best, 1])
                    if (.Platform$OS.type == "windows") 
                      flush.console()
                  }
                  old.rate <- error.rates[best, 2]
                }
                else finished <- TRUE
            }
        }
        else {
            if (is.element(error.rates[best, 1], model)) {
                index <- is.element(model, error.rates[best, 1])
                out <- c(out, model[index])
                model <- model[!index]
                result.v <- rbind(result.v, c("out", out[length(out)]))
                result.e <- c(result.e, error.rates[best, 2])
                if (output) {
                  textoutput(error.rates[best, 2], model, 
                    out.of.model = error.rates[best, 1])
                  if (.Platform$OS.type == "windows") 
                    flush.console()
                }
            }
            else {
                model <- c(model, error.rates[best, 1])
                result.v <- rbind(result.v, c("in", model[length(model)]))
                result.e <- c(result.e, error.rates[best, 2])
                out <- out[!is.element(out, error.rates[best, 1])]
                if (output) {
                  textoutput(error.rates[best, 2], model, 
                    into.model = error.rates[best, 1])
                  if (.Platform$OS.type == "windows") 
                    flush.console()
                }
            }
            old.rate <- error.rates[best, 2]
        }
        if (is.finite(maxvar)) 
            finished <- (length(model) >= maxvar)
    }
    runtime <- min.sec(proc.time()[3] - runtime)
    if (output) {
        cat("\n")
        print(runtime)
        cat("\n")
    }
    object <- try(do.call(method, list(data[, model], grouping, ...)), 
        silent = TRUE)
    if (class(object) != "try-error") {
        classi <- try(predict(object, data[, model]), silent = TRUE)
        if (class(classi) != "try-error") {
            if (is.list(classi)) 
                classi <- classi$posterior
        aper <- (1-ucpm(classi,grouping)[[criterion]])
        }
        else aper <- NA
    }
    else aper <- NA
    model <- sort(model)
  result <- list("call" = match.call(), "method" = method, 
        "start.variables" = start.vars, 
        "process" = cbind.data.frame("step" = result.v[,1], 
                                     "var" = as.numeric(result.v[,2]), 
                                     "varname" = c("--", varnames[as.numeric(result.v[-1, 2])]),
                                     "result.pm" = 1 - result.e),
        "model" = cbind.data.frame(  "nr" = model, 
                                     "name" = I(varnames[model])),
        "result.pm" = c("crossval" = 1 - result.e[length(result.e)], "apparent" = aper),
        "runtime" = runtime, "performance.measure" = cr) #, "cv.groups"=cv.groups)
  rownames(result$process) <- as.character(0:(length(result.v[ ,1]) - 1))
  if(length(result$start.variables)) 
    names(result$start.variables) <- varnames[result$start.variables]
  class(result) <- "stepclass"
  return(result)
}


print.stepclass <- function(x,...)
{
  kommalist <- function(charvec)
  {
    if (length(charvec) == 1) cat(charvec)
    else cat(paste(charvec[-length(charvec)], ",", sep = ""), charvec[length(charvec)])
  }
  cat("method      :", x$method, "\n")
  cat("final model : ")
  kommalist(x$model$name)
  cat("\n")
  cat(x$performance.measure, "=", as.character(signif(x$result.pm[1],4)), "\n")
  invisible(x)
}

plot.stepclass <- function(x, ...)
{
  signum <- rep("-", length(x$process$var)-1)
  signum[as.character(x$process$step[-1]) == "in"] <- "+"
  change <- c("START", paste(signum, x$process$varname[-1]))
  par(mar = c(10, 4, 4, 2) + 0.1)
  plot(seq(along = x$process[,1]), x$process$result.pm, type = "b", 
       xlab = "", ylab = paste("estimated", x$performance.measure), xaxt = "n", ...)
  axis(1, at = seq(along = x$process$result.pm), labels = change, las = 3, ...)
  invisible(x$progress)
}
svmlight<-function (x, ...) 
    UseMethod("svmlight")


svmlight.default <- function(x, grouping, temp.dir=NULL, pathsvm=NULL, del=TRUE, 
    type="C", class.type = "oaa" ,svm.options=NULL,  prior=NULL, out=FALSE, ...)
{
y <- grouping
 pick<-NULL     
 if (type!="R")
   {
    ### Construct Dummymatrix for 1-a-a classification: 
    ### ncol= no. of classes, (i,j)=+1 iff obj i comes from class j, -1 else 
    if(is.null(prior)) prior <- table(y) / length(y)
    ys <- as.factor(y)
    tys <- table(ys)
    lev <- levels(ys)
    if (class.type !="oao")
      {
       class.type<-"oaa"
       ymat <- matrix(-1, nrow = nrow(x), ncol = length(tys))
       ymat[cbind(seq(along = ys), sapply(ys, function(x) which(x == lev)))] <- 1
      }
     else
       { ## Classification: one against one
        nclass <- length(table(levels(y)))
        m <- (nclass - 1)
        minus <- nclass + 1 - sequence(m:1)
        plus <- rep(1:m, m:1)
        pick <- rbind(plus, minus)
        xsplit <- split(data.frame(x), ys)
        ymat <- list()
        xlist <- list()
        for(k in 1:ncol(pick)){
            ymat[[k]] <- c(rep(1, nrow(xsplit[[ pick[1, k] ]])), rep(-1, nrow(xsplit[[ pick[2, k] ]])))
            xlist[[k]] <- rbind(xsplit[[ pick[1, k] ]], xsplit[[ pick[2, k] ]])
            }
        } 
    counts <- as.vector(tys)     
  }
 else
   {
   ### Regression: So transform vector to 1-col matrix and set options, J ...
     ymat <- matrix(y,ncol=1)
     lev <- NULL
     counts <- 1
     svm.options <- paste("-z r ",svm.options)
     J <- 1 
     }
    svm.model <- list()

### call svm_learn
    cmd <- if (is.null(pathsvm)) 
        "svm_learn"
    else file.path(pathsvm, "svm_learn")
    
    if(is.matrix(ymat)) J <- 1:ncol(ymat)
    if(is.list(ymat)) J <- 1:length(ymat)
    
### "paste" file names for training data and model
    train.filename <- paste(temp.dir, "_train_", J, ".dat", sep = "")
    model.filename <- paste(temp.dir, "_model_", J, ".txt", sep = "")
    PWin <- .Platform$OS.type == "windows"
 for (j in J){
      ### construct matrix to learn svm in a format, so that svmlight can read it
      if (class.type !="oao")
         {
          train <- svmlight.file(cbind(ymat[,j], x), train = TRUE)    
          }
        else
          {
           train <- svmlight.file(cbind(ymat[[j]], xlist[[j]]), train = TRUE)
           }  
      ### save to disk
          write.table(train, file = train.filename[j], row.names = FALSE, 
          col.names = FALSE, quote = FALSE)
          if (PWin) 
              system(paste(cmd, svm.options, train.filename[j], model.filename[j]), 
                  show.output.on.console = out)
          else 
              system(paste(cmd,svm.options, train.filename[j], model.filename[j]))
      ### store learned model
          svm.model[[j]] <- readLines(model.filename[j])
      }
    if (del) 
       file.remove(c(train.filename, model.filename))

    cl <- match.call()
    cl[[1]] <- as.name("svmlight")
    structure(list(prior = prior, counts = counts, lev = lev, 
        temp.dir = temp.dir, pathsvm = pathsvm, del = del, 
        type=type, class.type=class.type, pick=pick, J=J,
        svm.model = svm.model, svm.options = svm.options, call = cl), class = "svmlight")
}


#### svmlight interface for different calls. Copied and adapted from lda.xxx
svmlight.formula <- function(formula, data = NULL, ..., subset, na.action = na.fail) 
{
    m <- match.call(expand.dots = FALSE)
    if (is.matrix(eval.parent(m$data))) 
        m$data <- as.data.frame(data)
    m$... <- NULL
    m[[1]] <- as.name("model.frame")
    m <- eval.parent(m)
    Terms <- attr(m, "terms")
    grouping <- model.response(m)
    x <- model.matrix(Terms, m)
    xvars <- as.character(attr(Terms, "variables"))[-1]
    if ((yvar <- attr(Terms, "response")) > 0) 
        xvars <- xvars[-yvar]
    xlev <- if (length(xvars) > 0) {
        xlev <- lapply(m[xvars], levels)
        xlev[!sapply(xlev, is.null)]
    }
    xint <- match("(Intercept)", colnames(x), nomatch = 0)
    if (xint > 0) 
        x <- x[, -xint, drop = FALSE]
    res <- svmlight.default(x, grouping, ...)
    res$terms <- Terms
    cl <- match.call()
    cl[[1]] <- as.name("svmlight")
    res$call <- cl
    res$contrasts <- attr(x, "contrasts")
    res$xlevels <- xlev
    attr(res, "na.message") <- attr(m, "na.message")
    if (!is.null(attr(m, "na.action"))) 
        res$na.action <- attr(m, "na.action")
    res
}

svmlight.matrix <- function(x, grouping, ..., subset, na.action = na.fail) 
{
    if (!missing(subset)) {
        x <- x[subset, , drop = FALSE]
        grouping <- grouping[subset]
    }
    if (!missing(na.action)) {
        dfr <- na.action(structure(list(g = grouping, x = x), 
            class = "data.frame"))
        grouping <- dfr$g
        x <- dfr$x
    }
    res <- svmlight.default(x, grouping, ...)
    cl <- match.call()
    cl[[1]] <- as.name("svmlight")
    res$call <- cl
    res
}

svmlight.data.frame <- function (x, ...) 
{
   res <- svmlight.matrix(structure(data.matrix(x), class = "matrix"), 
        ...)
    cl <- match.call()
    cl[[1]] <- as.name("svmlight")
    res$call <- cl
    res
}


### predict method for svmlight
predict.svmlight <- function(object, newdata, scal=TRUE,...)
{
### utility function for oao classification
sf<-function(x,pick)
{
erg<-numeric(max(pick))
names(erg)<-1:max(pick)
dummy<-table(diag(pick[(x<0)+1,]))
erg[names(dummy)]<-dummy
return(erg)
}




### copied and adapted form predict.lda
    if (!inherits(object, "svmlight")) 
        stop("object not of class svmlight")
    if (!is.null(Terms <- object$terms)) {
        if (missing(newdata)) 
            newdata <- model.frame(object)
        else {
            newdata <- model.frame(as.formula(delete.response(Terms)), 
                newdata, na.action = function(x) x, xlev = object$xlevels)
        }
        x <- model.matrix(delete.response(Terms), newdata, contrasts = object$contrasts)
        xint <- match("(Intercept)", colnames(x), nomatch = 0)
        if (xint > 0) 
            x <- x[, -xint, drop = FALSE]
    }
    else {
        if (missing(newdata)) {
            if (!is.null(sub <- object$call$subset)) 
                newdata <- eval.parent(parse(text = paste(deparse(object$call$x, 
                  backtick = TRUE), "[", deparse(sub, backtick = TRUE), 
                  ",]")))
            else newdata <- eval.parent(object$call$x)
            if (!is.null(nas <- object$call$na.action)) 
                newdata <- eval(call(nas, newdata))
        }
        if (is.null(dim(newdata))) 
            dim(newdata) <- c(1, length(newdata))
        x <- as.matrix(newdata)
    }
#######################################

### save test data on disk
    x <- svmlight.file(cbind(rep(0,nrow(x)),x), train = TRUE)
    #x <- svmlight.file(x, train = FALSE)
    test.filename <- paste(object$temp.dir, "_test_.dat", sep = "")
    write.table(x, file = test.filename, row.names = FALSE, 
        col.names = FALSE, quote = FALSE)
    werte <- NULL
    
    cmd <- if (is.null(object$pathsvm)) 
        "svm_classify"
    else file.path(object$pathsvm, "svm_classify")
    J <- object$J
    model.filename <- paste(object$temp.dir, "_model_", J, ".txt", sep = "")
    pred.filename <- paste(object$temp.dir, "_pred_", J, ".txt", sep = "")
### for all learned models predict call svm_classify
    for (j in J){
        writeLines(object$svm.model[[j]], model.filename[j])
        system(paste(cmd, test.filename, model.filename[j], pred.filename[j]))
    ### read predicted values
        prognose <- read.table(pred.filename[j], header = FALSE)[ , 1]
        werte <- cbind(werte, prognose)
    }
    if (object$del) 
        file.remove(c(test.filename, pred.filename, model.filename))
    if (object$type=="C")
    {
    ### Classification: choose class with highest decision value f(x)
    if  (object$class.type=="oao")   
        {
        werte2<-apply(werte,1,sf,pick=object$pick) 
        werte<-t(werte2)
        }
    
    
    
    classes <- factor(max.col(werte), levels = seq(along = object$lev), 
        labels = object$lev)
    colnames(werte)<-object$lev 
    
    if (scal) werte<-e.scal(werte)$sv
    return(list(class = classes, posterior = werte))
  }
   else {return(as.vector(werte))} ### Regression
}


svmlight.file <- function(x, train = FALSE,...)
{
    if(is.vector(x)) x <- t(x)
    erg <- x
    sn <- 1:nrow(x) 
    if(!train) erg[sn, 1] <- paste("1:", x[sn, 1], sep = "")
    if(ncol(x) > 1){
        j <- 2:ncol(x)
        erg[ , -1] <- matrix(paste(j - train, t(x[,j]), sep = ":"), ncol = ncol(x)-1, byrow = TRUE)
    }
    return(erg)
} 
tritrafo <- function(x, y=NULL, z=NULL, check=TRUE, tolerance=0.0001)
# projects 3D-mixture onto 2D-triangle 
{
  projector <- cbind( "8.am"  =c(cos((2*pi)*(7/12)), sin((2*pi)*(7/12)))*(2/3),
                     "12.noon"=c(0,2/3),
                      "4.pm"  =c(cos((2*pi)*(11/12)), sin((2*pi)*(11/12)))*(2/3))
  trafo <- function(mix, projct=projector)
  {
    return(projct %*% (mix - rep(1/3, 3)))
  }
  if (is.matrix(x)) dat <- x
  else if (is.null(y)) dat <- t(x)
       else {
         if (is.null(z)) z <- 1-(x+y)
         dat <- cbind(x,y,z)
       }
  result <- t(apply(dat,1,trafo))
  colnames(result) <- c("x","y")
  if (check) {
    valid <- apply(dat,1,function(x)all(is.finite(x)))  # ignore NAs etc. 
    if (any(dat[valid,]<0)) warning("negative components")
    if (any(valid) && !identical(all.equal(rowSums(matrix(dat[valid,],ncol=3)), rep(1,sum(valid)), tolerance=tolerance), TRUE)) 
      warning("components do not sum to one")
    }
  return(result)
}

trilines <- function(x, y=NULL, z=NULL, ...)
{
  result <- tritrafo(x,y,z)
  lines(result[,1], result[,2], ...)
  invisible(result)
}

tripoints <- function(x, y=NULL, z=NULL, ...)
{
  result <- tritrafo(x,y,z)
  points(result[,1], result[,2], ...)
  invisible(result)
}

trigrid <- function(x=seq(0.1,0.9,by=0.1), y=NULL, z=NULL, lty="dashed", col="grey", ...)
{
  makegrid<-function(val, dim=1)
  {
    if (length(val)>0) {
      col2 <- 1- val
      col1 <- rep(val, rep(2,length(val)))
      col2 <- rep(col2, rep(2,length(col2)))
      col3 <- col2
      col2[c(TRUE,FALSE)] <- 0
      col3[c(FALSE,TRUE)] <- 0
      permu <- cbind(1:3, c(2,1,3), c(2,3,1))
      coords <- cbind(col1, col2, col3)[,permu[,dim]]
      result <- matrix(NA, ncol=3, nrow=length(val)*3)
      result[c(TRUE,TRUE,FALSE),] <- coords
    }
    else result <- NULL
    return(result)
  }
  if (is.null(x)) grx <- rep(NA,3)
  else {
    x <- sort(unique(x[(x>=0) & (x<1)]))
    grx <- makegrid(x,1)
  }
  if (is.null(y)) {
    gridlines <- rbind(grx, rep(NA,3), grx[,c(2,1,3)], rep(NA,3), grx[,c(2,3,1)])
  }
  else {
    y <- sort(unique(y[(y>=0) & (y<1)]))
    z <- sort(unique(z[(z>=0) & (z<1)]))
    gridlines <- rbind(grx, rep(NA,3), makegrid(y,2), rep(NA,3), makegrid(z,3))
  }
  trilines(gridlines, lty=lty, col=col, ...)
  invisible(gridlines)
}

triframe <- function(label=1:3, label.col=1, cex=1,...)
{
  shift <- 1.1
  corners <- tritrafo(x=diag(3))
  if(length(label)==3) {
    text(corners[1,1]*shift, corners[1,2]*shift, label[1], adj=c(1/3,1), col=label.col, cex=1)
    text(corners[2,1]*shift, corners[2,2]*shift, label[2], adj=c(0.5,0), col=label.col, cex=1)
    text(corners[3,1]*shift, corners[3,2]*shift, label[3], adj=c(2/3,1), col=label.col, cex=1)
  }
  invisible(trilines(diag(3)[c(1,2,3,1),],...))
}

triplot <- function(x=NULL, y=NULL, z=NULL, main="",
                    frame=TRUE, label=1:3,
                    grid=seq(0.1,0.9,by=0.1), center=FALSE, set.par=TRUE, ...)
{
  margin <- c(0.1, 0.1, 0.1, 0.1) # bottom, left, top, right 
  corners <- tritrafo(x=diag(3))
  rownames(corners) <- paste("corner", 1:3, sep="")
  if (set.par) {
    if (main!="") newmar <- c(0, 0, 4, 0) + 0.1
    else newmar <- rep(0, 4) + 0.1
    par(mar=newmar)
  }
  plot.new()
  plot.window(xlim=c(corners[1,1]-margin[2], corners[3,1]+margin[4]),
              ylim=c(corners[1,2]-margin[1], corners[2,2]+margin[3]), asp=1)
  # grid: 
  if (!(is.logical(grid))) trigrid(x=grid)
  else if (grid) trigrid(x=seq(0.1,0.9,by=0.1))
  # centerlines 
  if (center) trilines(centerlines(3))
  # outer triangle & corner labels: 
  if (frame) triframe(label=label)
  # main title: 
  if (main!="") title(main=main)
  # points (if supplied): 
  if (!is.null(x))
    tripoints(x,y,z,...)
  invisible(corners)
}

triperplines <- function(x, y=NULL, z=NULL, lcol="red", pch=17, ...)
{
  if (all(is.null(c(y,z)))) point <- x[1:3]
  else if (all(is.null(z))) point <- c(x[1], y[1], 1-x[1]-y[1])
  else point <- c(x[1], y[1], z[1])
  projektor <- cbind("10.am" = -c(cos((2*pi)*(7/12)), sin((2*pi)*(7/12))),
                      "6.pm" = -c(0,1),
                      "2.pm" = -c(cos((2*pi)*(11/12)), sin((2*pi)*(11/12))))
  tpoint <- tritrafo(point)[1,]
  footlines <- rbind(tpoint,
                     tpoint+projektor[,1]*point[1],
                     rep(NA,2),
                     tpoint,
                     tpoint+projektor[,2]*point[2],
                     rep(NA,2),
                     tpoint,
                     tpoint+projektor[,3]*point[3])
  rownames(footlines) <- NULL
  lines(footlines, col=lcol, ...)
  if(pch) tripoints(point, pch = pch, ...)
  invisible(footlines)
}
ucpm <- function(m, tc, ec = NULL)
{
# membership values
# true classes
# estimated classes
    klassen <- function(y){
        n <- length(y)
        m <- length(table(y))
        leer <- numeric(n*m)
        Y <- matrix(leer, nrow=n, ncol=m)
        for(i in 1:n) Y[i,y[i]] <- 1
        return(Y)
    }
    if(is.null(colnames(m))) colnames(m) <- levels(tc)
    membercheck(m)
    G <- ncol(m)
    N <- nrow(m)
    dummy <- matrix(0, nrow=N, ncol=G)
    if (is.null(ec)) 
        ec <- factor(max.col(m), levels = seq(along = colnames(m)), 
            labels = colnames(m))
    CR <- mean(ec == tc)
    e.ec <- klassen(ec)
    e.c <- klassen(tc)
    dummy <- apply((e.ec - m)^2, 1, sum)
    AS <- 1 - sqrt(G) / (N * sqrt(G-1)) * sum(sqrt(dummy))
    dummy <- apply((e.c - m)^2, 1, sum)
    AC <- 1 - sqrt(G) / (N * sqrt(G-1)) * sum(sqrt(dummy))
    CF <- mean(apply(m,1,max))
    CFvec <- apply(m, 1, max)
    CFvec <- unlist(lapply(split(CFvec,tc), mean))
    return(list(CR=CR, AC=AC, AS=AS, CF=CF, CFvec=CFvec))
}
