bincombinations <- function(p) {

  retval <- matrix(0, nrow=2^p, ncol=p)
  
  for(n in 1:p){
    retval[,n] <- rep(c(rep(0, (2^p/2^n)), rep(1, (2^p/2^n))),
                      length=2^p)
  }
  retval
}



cmeanscl <- function (x, centers, iter.max = 100, verbose = FALSE,
                    dist = "euclidean", method = "cmeans",
                    m=2, rate.par = NULL) 
{
  xrows <- dim(x)[1]
  xcols <- dim(x)[2]
  xold <- x
  perm <- sample(xrows)
  x <- x[perm, ]
  ## initial values are given
  if (is.matrix(centers)) 
    ncenters <- dim(centers)[1]
  else {
   ## take centers random vectors as initial values
    ncenters <- centers
    centers <- x[rank(runif(xrows))[1:ncenters], ]+0.001
  }
  
  dist <- pmatch(dist, c("euclidean", "manhattan"))
  if (is.na(dist)) 
    stop("invalid distance")
  if (dist == -1) 
    stop("ambiguous distance")
  
  method <- pmatch(method, c("cmeans", "ufcl"))
  if (is.na(method)) 
    stop("invalid clustering method")
  if (method == -1) 
    stop("ambiguous clustering method")
  
  if (method == 2) {
    if (missing(rate.par)) {
      rate.par <- 0.3
    }
  }
            
  initcenters <- centers
  ## dist <- matrix(0, xrows, ncenters)
  ## necessary for empty clusters
  pos <- as.factor(1:ncenters)
  rownames(centers) <- pos
  iter <- integer(1)
  
  if (method == 1) {
    retval <- .C("cmeans", xrows = as.integer(xrows),
                 xcols = as.integer(xcols), 
                 x = as.double(x), ncenters = as.integer(ncenters), 
                 centers = as.double(centers), 
                 iter.max = as.integer(iter.max), iter = as.integer(iter), 
                 verbose = as.integer(verbose), dist = as.integer(dist-1), 
                 U=double(xrows*ncenters), UANT=double(xrows*ncenters),
                 m=as.double(m), ermin=double(1))
  }
  else if (method == 2) {
    retval <- .C("ufcl", xrows = as.integer(xrows),
                 xcols = as.integer(xcols), 
                 x = as.double(x), ncenters = as.integer(ncenters), 
                 centers = as.double(centers), 
                 iter.max = as.integer(iter.max), iter = as.integer(iter), 
                 verbose = as.integer(verbose), dist = as.integer(dist-1), 
                 U=double(xrows*ncenters), m=as.double(m),
                 rate.par = as.double(rate.par), ermin=double(1))
  }
  
  centers <- matrix(retval$centers, ncol = xcols, dimnames = dimnames(initcenters))
  
  U <- retval$U
  U <- matrix(U, ncol=ncenters)
#  clusterU <- max.col(U)
   clusterU <- apply(U,1,which.max)
  clusterU <- clusterU[order(perm)]
  U <- U[order(perm),]
  
  clustersize <- as.integer(table(clusterU))
  
  retval <- list(centers = centers, cluster = clusterU,
                 size = clustersize, dist = dist,
                 m=retval$m, member=U, withinsd = retval$ermin,
                 learning = list(ncenters = ncenters,
                   initcenters = initcenters, iter = retval$iter - 1,
                   rate.par = rate.par), call = match.call())
  
  class(retval) <- c("fclust")
  return(retval)
}





countpattern <- function(x, matching=FALSE)
{
     nvar <- dim(x)[2]
     n <- dim(x)[1]

     ## build matrix of all possible binary vectors
     b <- matrix(0, 2^nvar, nvar)
     for (i in 1:nvar)
         b[, nvar+1-i] <- rep(rep(c(0,1),c(2^(i-1),2^(i-1))),2^(nvar-i))

     namespat <- b[,1]
     for (i in 2:nvar)
         namespat <- paste(namespat, b[,i], sep="")

     xpat <- x[,1]
     for (i in 2:nvar)
         xpat <- 2*xpat+x[,i]
     xpat <- xpat+1

     pat <- tabulate(xpat, nb=2^nvar)
     names(pat) <- namespat

     if (matching)
         return(list(pat=pat, matching=xpat))
     else
         return(pat)
 }
cshell <- function (x, centers, iter.max = 100, verbose = FALSE,
                    dist = "euclidean", method = "cshell",
                    m=2, radius= NULL) 
{
  xrows <- dim(x)[1]
  xcols <- dim(x)[2]
  xold <- x
  perm <- sample(xrows)
  x <- x[perm, ]
  ## initial values are given
  if (is.matrix(centers)) 
    ncenters <- dim(centers)[1]
  else {
   ## take centers random vectors as initial values
    ncenters <- centers
    centers <- x[rank(runif(xrows))[1:ncenters], ]+0.001
  }

  ##initialize radius
  if (missing(radius))
    radius <- rep(0.2,ncenters)
  else
    radius <- as.double(radius)

  dist <- pmatch(dist, c("euclidean", "manhattan"))
  if (is.na(dist)) 
    stop("invalid distance")
  if (dist == -1) 
    stop("ambiguous distance")
  
  method <- pmatch(method, c("cshell"))
  if (is.na(method)) 
    stop("invalid clustering method")
  if (method == -1) 
    stop("ambiguous clustering method")
  
  initcenters <- centers
  ## dist <- matrix(0, xrows, ncenters)
  ## necessary for empty clusters
  pos <- as.factor(1:ncenters)
  rownames(centers) <- pos
  iter <- integer(1)

  flag <- integer(1)

  
  retval <- .C("cshell", xrows = as.integer(xrows),
               xcols = as.integer(xcols), 
               x = as.double(x), ncenters = as.integer(ncenters), 
               centers = as.double(centers), 
               iter.max = as.integer(iter.max), iter = as.integer(iter), 
               verbose = as.integer(verbose), dist = as.integer(dist-1), 
               U=double(xrows*ncenters), UANT=double(xrows*ncenters),
               m=as.double(m), ermin=double(1), radius=as.double(radius),
               flag=as.integer(flag))

  centers <- matrix(retval$centers, ncol = xcols, dimnames = dimnames(initcenters))
  
  
  radius <- as.double(retval$radius)
  U <- retval$U
  U <- matrix(U, ncol=ncenters)
  UANT <- retval$UANT
  UANT <- matrix(UANT, ncol=ncenters)

  
  iter <- retval$iter
  flag <- as.integer(retval$flag)
  
  
  ##Optimization part
  while (((flag == 1) || (flag==4)) && (iter<=iter.max)){
    
    flag <- 3
    
    
    system <- function (spar=c(centers,radius), x, U, m, i){
      k <- dim(x)[1]
      d <- dim(x)[2]
      nparam<-length(spar)
      
      v<-spar[1:(nparam-1)]
      r<-spar[nparam]
      
      ##distance matrix x_k - v_i
      distmat <- t(t(x)-v)
      
      ##norm from x_k - v_i
      normdist <- distmat[,1]^2
      for (j in 2:d)
        normdist<-normdist+distmat[,j]^2
      normdist <- sqrt(normdist)
      
      ##equation 5
      op <- sum( (U[,i]^m) * (normdist-r) )^2
      ##equation 4
      equationmatrix <- ((U[,i]^m) * (1-r/normdist))*distmat
      op<- op+apply(equationmatrix, 2, sum)^2
      
    }
    
    
    for (i in 1:ncenters){
      spar <- c(centers[i,],radius[i])
      npar <- length(spar)
      
      optimres <- optim(spar ,system, method="CG", x=x, U=U, m=m, i=i)
      centers[i,] <- optimres$par[1:(npar-1)]
      radius[i] <- optimres$par[npar]
    }
    
    
    retval <- .C("cshell", xrows = as.integer(xrows),
                 xcols = as.integer(xcols), 
                 x = as.double(x), ncenters = as.integer(ncenters), 
                 centers = as.double(centers), 
                 iter.max = as.integer(iter.max), iter = as.integer(iter-1), 
                 verbose = as.integer(verbose), dist = as.integer(dist-1), 
                 U=as.double(U), UANT=as.double(UANT),
                 m=as.double(m), ermin=double(1), radius=as.double(radius),
                 flag=as.integer(flag))
    
    flag<-retval$flag
    if (retval$flag!=2)
      flag<-1
    
    
    centers <- matrix(retval$centers, ncol = xcols, dimnames = dimnames(initcenters))
    
    radius <- as.double(retval$radius)
    U <- retval$U
    U <- matrix(U, ncol=ncenters)
    UANT <- retval$UANT
    UANT <- matrix(UANT, ncol=ncenters)
    
    
    
    iter <- retval$iter
    
  }
  
  centers <- matrix(retval$centers, ncol = xcols, dimnames = dimnames(initcenters))
  
  U <- retval$U
  U <- matrix(U, ncol=ncenters)
  
  
  clusterU <- apply(U,1,which.max)
  clusterU <- clusterU[order(perm)]
  U <- U[order(perm),]
  
  clustersize <- as.integer(table(clusterU))
  radius <- as.double(retval$radius)
  
  retval <- list(centers = centers, cluster = clusterU,
                 size = clustersize, dist = dist,
                 m=retval$m, member=U, withinsd = retval$ermin,
                 learning = list(ncenters = ncenters,
                   initcenters = initcenters, iter = retval$iter - 1,
                   radius=radius), call = match.call())
  
  class(retval) <- c("fclust")
  return(retval)
}
  














rdiscrete <- function (n, probs, values = 1:length(probs), method="inverse",
                       aliasmatrix = NULL)
{
    
    if (length(probs) != length(values))
        stop("rdiscrete: probs and values must have the same length.")
    if (sum(probs < 0) > 0)
        stop("rdiscrete: probs must not contain negative values.")
    
    if (n == 1)
        return (values[sum(runif(1) > p) + 1])
    else
    {
        method <- pmatch(method, c("inverse", "alias"))
        if (is.na(method))
            stop("rdiscrete: unknown method.")
        if (method == 1)
        {
            p <- cumsum(probs)/sum(probs)
            l <- length(probs)
            m <- numeric(n)
            a <- runif(n)
            for (i in 1:n)
                m[i] <- sum(a[i] > p)
            return(values[m + 1])
        }
        else
        {
            if (missing(aliasmatrix))
                aliasmatrix <- aliasmat(probs)
            return(rdiscrete.alias(n, probs, aliasmatrix))
        }
    }
}

aliasmat <- function(p)
{
    p <- p / sum(p)
    q <- p * (pn <- length(p))
    r <- matrix(0, nrow=pn, ncol=3)

    eps <- .Machine$double.eps
    while (sum(!is.na(q))>1)
    {
        qklein <- min((1:pn)[q<=1+eps],na.rm=TRUE)
        qgross <- max((1:pn)[q>=1-eps],na.rm=TRUE)
        r[qklein,] <- c(q[qklein],qklein,qgross)
        q[qgross] <- q[qgross] + q[qklein] - 1
        q[qklein] <- NA
    }
    qmittel <- (1:pn)[!is.na(q)]
    r[qmittel,] <-  c(q[qmittel],qmittel,qmittel)

    return(r)
}

rdiscrete.alias <- function(n, p, aliasmatrix = aliasmat(p))
{
    x <- sample(1:nrow(aliasmatrix), n, replace=TRUE)
    y <- runif(n)
    
    retval <- rep(0, length=n)
    
    eins <- (y <= aliasmatrix[x,1])
    
    retval[eins] <- (1:length(p))[aliasmatrix[x[eins],2]]
    retval[!eins] <- (1:length(p))[aliasmatrix[x[!eins],3]]
    
    return(retval)
}

    
    
aliasmat2prob <- function(r)
{
    p <- rep(0, length = length(unique(r[,2:3])))
    names(p) <- (pnames <- sort(unique(r[,2:3])))
    
    for(n in pnames){
        if(any(r[,2]==n)){
            p[n] <- r[r[,2]==n,1]
        }
        if(any(r[,3]==n)){
            p[n] <- p[n] + sum(1 - r[r[,3]==n,1])
        }
    }
    p<- p/length(p)
    p
}



ddiscrete <- function (x, probs, values = 1:length(probs))
{
    
    if (length(probs) != length(values))
        stop("ddiscrete: probs and values must have the same length.")
    if (sum(probs < 0) > 0)
        stop("ddiscrete: probs must not contain negative values.")
    if (!is.array(x) && !is.vector(x) && !is.factor(x))
        stop("ddiscrete: x must be an array or a vector or a factor.")
    
    p <- probs/sum(probs)
    
    y <- as.vector(x)
    l <- length(y)
    z <- rep(0,l)
    
    for (i in 1:l)
        if (any(values == y[i]))
            z[i] <- p[values == y[i]]
    
    z <- as.numeric(z)
    if (is.array(x))
        dim(z) <- dim(x)
    
    return(z)
}


pdiscrete <- function (q, probs, values = 1:length(probs))
{
    
    if (length(probs) != length(values))
        stop("pdiscrete: probs and values must have the same length.")
    if (sum(probs < 0) > 0)
        stop("pdiscrete: probs must not contain negative values.")
    if (!is.array(q) & !is.vector(q))
        stop("pdiscrete: q must be an array or a vector")
    
    p <- probs/sum(probs)
    
    y <- as.vector(q)
    l <- length(y)
    z <- rep(0,l)
    
    for (i in 1:l)
        z[i] <- sum(p[values <= y[i]])
    
    z <- as.numeric(z)
    if (is.array(q))
        dim(z) <- dim(q)
    
    return(z)
}

qdiscrete <- function (p, probs, values = 1:length(probs))
{
    
    if (length(probs) != length(values))
        stop("qdiscrete: probs and values must have the same length.")
    if (sum(probs < 0) > 0)
        stop("qdiscrete: probs must not contain negative values.")
    if (!is.array(p) & !is.vector(p))
        stop("qdiscrete: p must be an array or a vector")
    
    probs <- cumsum(probs)/sum(probs)
    
    y <- as.vector(p)
    l <- length(y)
    z <- rep(0,l)
    
    for (i in 1:l)
        z[i] <- length(values) - sum(y[i] <= probs) + 1
    
    z <- as.numeric(z)
    z <- values[z]
    if (is.array(p))
        dim(z) <- dim(p)
    
    return(z)
  }



element <- function(x, i) {

  if(!is.array(x))
    stop("x is not an array")
  
  ni <- length(i)
  dx <- dim(x)
  
  if(length(i)!=length(dx))
    stop("Wrong number of subscripts")

  if(ni==1){
    return(x[i])
  }
  else{
    m1 <- c(i[1], i[2:ni]-1)
    m2 <- c(1,cumprod(dx)[1:(ni-1)])
    return(x[sum(m1*m2)])
  }
}

hamming.distance <- function(x,y){
  
  z<-NULL
  
  if(is.vector(x) && is.vector(y)){
    z <- sum(as.logical(x) != as.logical(y))
  }
  else{
    z <- matrix(0,nrow=nrow(x),ncol=nrow(x))
    for(k in 1:(nrow(x)-1)){
      for(l in (k+1):nrow(x)){
	z[k,l] <- hamming.distance(x[k,], x[l,])
	z[l,k] <- z[k,l]
      }
    }
    dimnames(z) <- list(dimnames(x)[[1]], dimnames(x)[[1]])
  }
  z
}



















hamming.window <- function (n)
  {
    if (n == 1)
      c <- 1
    else
      {
	n <- n-1
	c <- 0.54 - 0.46*cos(2*pi*(0:n)/n)
      }
    return(c)
  }
hanning.window <- function (n)
  {
    if (n == 1)
      c <- 1
    else
      {
	n <- n-1
	c <- 0.5 - 0.5*cos(2*pi*(0:n)/n)
      }
    return(c)
  }
ica <- function(X, lrate, epochs=100, ncomp=dim(X)[2], 
                      fun="negative")
  {
    if (!is.matrix(X))
      {
        if (is.data.frame(X))
          X <- as.matrix(X)
        else
          stop("ica: X must be a matrix or a data frame")
      }
    if (!is.numeric(X))
      stop("ica: X contains non numeric elements")
            
    m <- dim(X)[1]
    n <- dim(X)[2]

    Winit <- matrix(rnorm(n*ncomp), ncomp, n)
    W <- Winit

    if (!is.function(fun))
      {
        funlist <- c("negative kurtosis", "positive kurtosis",
                     "4th moment")
        p <- pmatch(fun, funlist)
        if (is.na(p))
          stop("ica: invalid fun")
        funname <- funlist[p]
        if (p == 1) fun <- tanh
        else if (p == 2) fun <- function(x) {x - tanh(x)}
        else if (p == 3) fun <- function(x) {sign(x)*x^2}
      }
    else funname <- as.character(substitute(fun))
    
    for (i in 1:epochs)
      for (j in 1:m)
        {
          x <- X[j,, drop=F]
          y <- W%*%t(x)
          gy <- fun(y)
          W <- W + lrate*gy%*%(x-t(gy)%*%W)
        }
    colnames(W) <- NULL
    pr <- X%*%t(W)
    retval <- list(weights = W, projection = pr, epochs = epochs,
                fun = funname, lrate = lrate, initweights = Winit)
    class(retval) <- "ica"
    return(retval)
  }


print.ica <- function(x)
  {
    cat(x$epochs, "Trainingssteps with a learning rate of", x$lrate, "\n")
    cat("Function used:", x$fun,"\n\n")
    cat("Weightmatrix\n")
    print(x$weights)
  }

plot.ica <- function(x, ...) pairs(x$pr, ...)

impute <- function(x, what=c("median", "mean")){

    what <- match.arg(what)
    
    if(what == "median"){
        retval <-
            apply(x, 2,
                  function(z) {z[is.na(z)] <- median(z, na.rm=TRUE); z})
    }
    else if(what == "mean"){
        retval <-
            apply(x, 2,
                  function(z) {z[is.na(z)] <- mean(z, na.rm=TRUE); z})
    }
    else{
        stop("`what' invalid")
    }
    retval
}
interpolate <- function(x, a, adims=lapply(dimnames(a), as.numeric),
                        method="linear"){

  if(is.vector(x)) x<- matrix(x, ncol=length(x))

  if(!is.array(a))
    stop("a is not an array")

  ad <- length(dim(a))
  
  method <- pmatch(method, c("linear", "constant"))
  if (is.na(method)) 
    stop("invalid interpolation method")
  
  if(any(unlist(lapply(adims, diff))<0))
    stop("dimensions of a not ordered")

  retval <- rep(0, nrow(x))
  bincombi <- bincombinations(ad)
  
  for(n in 1:nrow(x)){

    ## the `leftmost' corner of the enclosing hypercube 
    leftidx <- rep(0, ad)
    xabstand <- rep(0, ad)
    aabstand <- rep(0, ad)

    for(k in 1:ad){
      if(x[n,k] < min(adims[[k]]) || x[n,k] > max(adims[[k]]))
        stop("No extrapolation allowed")
      else{
        leftidx[k] <- max(seq(adims[[k]])[adims[[k]] <= x[n,k]])
        ## if at the right border, go one step to the left
        if(leftidx[k] == length(adims[[k]]))
          leftidx[k] <- leftidx[k] - 1
        
        xabstand[k] <- x[n,k] - adims[[k]][leftidx[k]]
        aabstand[k] <- adims[[k]][leftidx[k]+1] - adims[[k]][leftidx[k]]
      }
    }

    coefs <- list()
    if(method==1){
      for(k in 1:(2^ad)){
        retval[n] <- retval[n] +
          element(a, leftidx+bincombi[k,]) *
            prod((aabstand-
                  convexcoeff(xabstand, aabstand*bincombi[k,]))/aabstand)
      }
    }
    else if(method==2){
      retval[n] <- element(a, leftidx)
    }
  }

  names(retval) <- rownames(x)
  retval
}
  
  

          

convexcoeff <- function(x, y) {
  for(k in 1:length(x)){
    if(y[k]>0)
      x[k] <- y[k]-x[k]
  }
  x
}
        
kurtosis <- function (x, na.rm = FALSE)
{
  if (na.rm) 
    x <- x[!is.na(x)]
  sum((x-mean(x))^4)/(length(x)*var(x)^2) - 3
}


lca <- function(x, k, niter=100, matchdata=FALSE, verbose=FALSE)
{

    ## if x is a data matrix -> create patterns
    if (is.matrix(x))
    {
        if (matchdata)
        {
            x <- countpattern(x, matching=TRUE)
            xmat <- x$matching
            x <- x$pat
        }
        else
            x <- countpattern(x, matching=FALSE)
    }
    else   ## if no data ist given, matchdata must be FALSE
        matchdata <- FALSE

    n <- sum(x)
    npat <- length(x)
    nvar <- round(log(npat)/log(2))
    
    ## build matrix of all possible binary vectors
    b <- matrix(0, 2^nvar, nvar)
    for (i in 1:nvar)
        b[, nvar+1-i] <- rep(rep(c(0,1),c(2^(i-1),2^(i-1))),2^(nvar-i))

    ## initialize probabilities
    classprob <- runif(k)
    classprob <- classprob/sum(classprob)
    names(classprob) <- 1:k
    p <- matrix(runif(nvar*k), k)

    
    pas <- matrix(0, k, npat)
    classsize <- numeric(k)
    
    for (i in 1:niter)
    {
        for (j in 1:k)
        {
            ## P(pattern|class)
            mp <- t(b)*p[j,]+(1-t(b))*(1-p[j,])
            pas[j,] <- drop(exp(rep(1,nvar)%*%log(mp))) # column product
        }
        ##  P(pattern|class)*P(class)
        pas <- t(t(pas)*classprob)        
        
        ## P(class|pattern)
        sump <- drop(rep(1,k)%*%pas)  # column sums
        pas <- t(t(pas)/sump)

        spas <- t(t(pas)*x)
        classsize <- drop(spas%*%rep(1,npat))  # row sums
        classprob <- classsize/n
        p <- pas%*%(x*b)/classsize
        if (verbose)
            cat("Iteration:", i, "\n")
    }

    for (j in 1:k)
    {
        mp <- t(b)*p[j,]+(1-t(b))*(1-p[j,])
        pas[j,] <- drop(exp(rep(1,nvar)%*%log(mp)))*classprob[j]
                                        # column product
    }

    ## LogLikelihood
    pmust <-  drop(rep(1,k)%*%pas)  # column sums
    ll <- sum(x*log(pmust))

    ## Likelihoodquotient
    xg0 <- x[x>0]
    ll0 <- sum(xg0*log(xg0/n))
    lq <- 2*(ll0-ll)

    ## bic
    bic <- -2*ll+log(n)*(k*(nvar+1)-1)
    bicsat <- -2*ll0+log(n)*(2^nvar-1)
    
    ## chisq
    ch <- sum((x-n*pmust)^2/(n*pmust))
    
    ## P(class|pattern)
    sump <- drop(rep(1,k)%*%pas)  # column sums
    pas <- t(t(pas)/sump)

    mat <- max.col(t(pas))
    if (matchdata)
        mat <- mat[xmat]

    colnames(p) <- 1:nvar
    rownames(p) <- 1:k
    
    lcaresult <- list(classprob=classprob, p=p, matching=mat,
                      logl=ll, loglsat=ll0,
                      chisq=ch, lhquot=lq, bic=bic, bicsat=bicsat, n=n,
                      np=(k*(nvar+1)-1), matchdata=matchdata)

    class(lcaresult) <- "lca"
    return(lcaresult)
}


print.lca <- function(l)
{
    cat("LCA-Result\n")
    cat("----------\n\n")

    cat("Datapoints:", l$n, "\n")
    cat("Classes:   ", length(l$classprob), "\n")
    cat("Probability of classes\n")
    print(round(l$classprob,3))

    cat("Itemprobabilities\n")
    print(round(l$p,2))
}

summary.lca <- function(l)
{
    nvar <- ncol(l$p)
    l$npsat <- 2^nvar-1
    l$df <- 2^nvar-1-l$np
    l$pvallhquot <- 1-pchisq(l$lhquot,l$df)
    l$pvalchisq <- 1-pchisq(l$chisq,l$df)
    l$k <- length(l$classprob)

    ## remove unnecessary list elements
    l$classprob <- NULL
    l$p <- NULL
    l$matching <- NULL
    
    class(l) <- "summary.lca"
    return(l)
}


print.summary.lca <- function(l)
{
    cat("LCA-Result\n")
    cat("----------\n\n")

    cat("Datapoints:", l$n, "\n")
    cat("Classes:   ", l$k, "\n")

    cat("\nGoodness of fit statistics:\n\n")
    cat("Number of parameters, estimated model:", l$np, "\n")
    cat("Number of parameters, saturated model:", l$npsat, "\n")
    cat("Log-Likelihood, estimated model:      ", l$logl, "\n")
    cat("Log-Likelihood, saturated model:      ", l$loglsat, "\n")

    cat("\nInformation Criteria:\n\n")
    cat("BIC, estimated model:", l$bic, "\n")
    cat("BIC, saturated model:", l$bicsat, "\n")

    cat("\nTestStatistics:\n\n")
    cat("Likelihood ratio:  ", l$lhquot,
        "  p-val:", l$pvallhquot, "\n")
    cat("Pearson Chi^2:     ", l$chisq,
        "  p-val:", l$pvalchisq, "\n")
    cat("Degress of freedom:", l$df, "\n")
}


bootstrap.lca <- function(l, nsamples=10, lcaiter=30, verbose=FALSE)
{

    n <- l$n
    classprob <- l$classprob
    nclass <- length(l$classprob)
    prob <- l$p
    nvar <- ncol(l$p)
    npat <- 2^nvar

    ## build matrix of all possible binary vectors
    b <- matrix(0, npat, nvar)
    for (i in 1:nvar)
        b[, nvar+1-i] <- rep(rep(c(0,1),c(2^(i-1),2^(i-1))),2^(nvar-i))

    ll <- lq <- ll0 <- ch <- numeric(nsamples)
    
    for (i in 1:nsamples)
    {
        ## generate data
        cm <- sample(1:nclass, size=n, replace=TRUE, prob=classprob)
        x <- matrix(runif(n*nvar), nrow=n)
        x <- (x<prob[cm,])+0
        x <- countpattern(x, matching=FALSE)

        if (verbose)
            cat ("Start of Bootstrap Sample", i, "\n")
        lc <- lca(x, nclass, niter=lcaiter, verbose=verbose)
        
        ll[i] <- lc$logl
        ll0[i] <- lc$loglsat
        lq[i] <- lc$lhquot
        ch[i] <- lc$chisq

        if (verbose)
            cat("LL:", ll[i], " LLsat:", ll0[i],
                " Ratio:", lq[i], " Chisq:", ch[i], "\n")
    }


    lratm <- mean(lq)
    lrats <- sd(lq)
    chisqm <- mean(ch)
    chisqs <- sd(ch)

    zl <- (l$lhquot-lratm)/lrats
    zc <- (l$ch-chisqm)/chisqs
    pzl <- 1-pnorm(zl)
    pzc <- 1-pnorm(zc)
    pl <- sum(l$lhquot<lq)/length(lq)
    pc <- sum(l$ch<ch)/length(ch)
    
    lcaboot <- list(logl=ll, loglsat=ll0, lratio=lq, 
                    lratiomean=lratm, lratiosd=lrats,
                    lratioorg=l$lhquot, zratio=zl,
                    pvalzratio=pzl, pvalratio=pl,
                    chisq=ch, chisqmean=chisqm, chisqsd=chisqs,
                    chisqorg=l$ch, zchisq=zc,
                    pvalzchisq=pzc, pvalchisq=pc,
                    nsamples=nsamples, lcaiter=lcaiter)
    class(lcaboot) <- "bootstrap.lca"
    return(lcaboot)
}

print.bootstrap.lca <- function(bl)
{
    cat("Bootstrap of LCA\n")
    cat("----------------\n\n")

    cat ("Number of Bootstrap Samples:    ", bl$nsamples, "\n")
    cat ("Number of LCA Iterations/Sample:", bl$lcaiter, "\n")
    
    cat("Likelihood Ratio\n\n")
    cat("Mean:", bl$lratiomean, "\n")
    cat("SDev:", bl$lratiosd, "\n")
    cat("Value in Data Set:", bl$lratioorg, "\n")
    cat("Z-Statistics:     ", bl$zratio, "\n")
    cat("P(Z>X):           ", bl$pvalzratio, "\n")
    cat("P-Val:            ", bl$pvalratio, "\n\n")

    cat("Pearson's Chisquare\n\n")
    cat("Mean:", bl$chisqmean, "\n")
    cat("SDev:", bl$chisqsd, "\n")
    cat("Value in Data Set:", bl$chisqorg, "\n")
    cat("Z-Statistics:     ", bl$zchisq, "\n")
    cat("P(Z>X):           ", bl$pvalzchisq, "\n")
    cat("P-Val:            ", bl$pvalchisq, "\n\n")
}

predict.lca <- function(lcares, x)
{
    if (lcares$matchdata)
        stop("predict.lca: only possible, if lca has been called with matchdata=FALSE")
    else
    {
        x <- countpattern(x, matching=TRUE)
        return(lcares$matching[x$matching])
    }
}
classAgreement <- function(tab){

    n <- sum(tab)
    ni <- apply(tab, 1, sum)  # row sums
    nj <- apply(tab, 2, sum)  # column sums

    m <- min(length(ni), length(nj))
    p0 <-  sum(diag(tab[1:m,1:m]))/n
    pc <- sum(ni[1:m]*nj[1:m])/n^2
    
    n2 <- choose(n,2)
    rand <- 1 + ( sum(tab^2) - (sum(ni^2)+sum(nj^2))/2 )/n2

    nis2 <- sum(choose(ni[ni>1],2))
    njs2 <- sum(choose(nj[nj>1],2))
    crand <- (sum(choose(tab[tab>1],2)) - (nis2*njs2)/n2)/
        ((nis2+njs2)/2 - (nis2*njs2)/n2)
    
    list(diag = p0, kappa = (p0-pc)/(1-pc), rand=rand, crand=crand)
}


matchClasses <- function(tab, method = "rowmax", iter=1, verbose=TRUE,
                         maxexact=9){

    methods <- c("rowmax", "greedy", "exact")
    method <- pmatch(method, methods)
    
    rmax <- apply(tab,1,which.max)
    myseq <- 1:ncol(tab)
    cn <- colnames(tab)
    rn <- rownames(tab)
    if(is.null(cn)){
        cn <- myseq
    }
    if(is.null(rn)){
        rn <- myseq
    }
    
    if(method==1){
        retval <- rmax
    }
    if(method==2 | method==3){
        if(ncol(tab)!=nrow(tab)){
            stop("Unique matching only for square tables.")
        }

        
        dimnames(tab) <- list(myseq, myseq)
        cmax <- apply(tab,2,which.max)
        retval <- rep(NA, ncol(tab))
        names(retval) <- colnames(tab)

        baseok <- cmax[rmax]==myseq
        for(k in myseq[baseok]){
            therow <- (tab[k,])[-rmax[k]]
            thecol <- (tab[, rmax[k]])[-k]            
            if(max(outer(therow, thecol, "+")) < tab[k, rmax[k]]){
                retval[k] <- rmax[k]
            }
            else{
                baseok[k] <- FALSE
            }
        }
        
                   

        if(verbose){
            cat("Direct agreement:", sum(baseok),
                "of", ncol(tab), "pairs\n")
        }

        if(!all(baseok)){
            if(method==3){
                if(sum(!baseok)>maxexact){
                    method <- 2
                    warning(paste("Would need permutation of", sum(!baseok),
                                  "numbers, resetting to greedy search\n"))
                }
                else{
                    iter <- gamma(ncol(tab)-sum(baseok)+1)
                    if(verbose){
                        cat("Iterations for permutation matching:", iter, "\n")
                    }
                    perm <- permutations(ncol(tab)-sum(baseok))
                }
            }
            
            ## rest for permute matching
            if(any(baseok)){
                rest <- myseq[-retval[baseok]]
            }
            else{
                rest <- myseq
            }
            
            for(l in 1:iter){
                newretval <- retval
                if(method == 2){
                    ok <- baseok
                    while(sum(!ok)>1){
                        rest <- myseq[!ok]
                        k <- sample(rest, 1)
                        if(any(ok)){
                            rmax <- tab[k, -newretval[ok]]
                        }
                        else{
                            rmax <- tab[k,]
                        }
                        newretval[k] <- as.numeric(names(rmax)[which.max(rmax)])
                        ok[k] <- TRUE
                    }
                    newretval[!ok] <- myseq[-newretval[ok]]
                }
                else{
                    newretval[!baseok] <- rest[perm[l,]]
                }
                
                if(l>1){
                    agree <- sum(diag(tab[,newretval]))/sum(tab)
                    if(agree>oldagree){
                        retval <- newretval
                        oldagree <- agree
                    }
                }
                else{
                    retval <- newretval
                    agree <- oldagree <- sum(diag(tab[,newretval]))/sum(tab)
                }
            }
        }
    }

    if(verbose){
        cat("Cases in matched pairs:",
            round(100*sum(diag(tab[,retval]))/sum(tab), 2), "%\n")
    }
            
    if(any(as.character(myseq)!=cn)){
        retval <- cn[retval]
    }
    names(retval) <- rn
    
    retval
}

compareMatchedClasses <- function(x, y, maxexact=9,
                                   method="rowmax", iter=1, verbose=FALSE)
{
    if(missing(y)){
        retval <- list(diag=matrix(NA, nrow=ncol(x), ncol=ncol(x)),
                       kappa=matrix(NA, nrow=ncol(x), ncol=ncol(x)),
                       rand=matrix(NA, nrow=ncol(x), ncol=ncol(x)),
                       crand=matrix(NA, nrow=ncol(x), ncol=ncol(x)))
        for(k in 1:(ncol(x)-1)){
            for(l in (k+1):ncol(x)){
                tab <- table(x[,k], x[,l])
                m <- matchClasses(tab, method=method, iter=iter,
                                  verbose=verbose, maxexact=maxexact)
                a <- classAgreement(tab[,m])
                retval$diag[k,l] <- a$diag
                retval$kappa[k,l] <- a$kappa
                retval$rand[k,l] <- a$rand
                retval$crand[k,l] <- a$crand
            }
        }
    }
    else{
        x <- as.matrix(x)
        y <- as.matrix(y)
        retval <- list(diag=matrix(NA, nrow=ncol(x), ncol=ncol(y)),
                       kappa=matrix(NA, nrow=ncol(x), ncol=ncol(y)),
                       rand=matrix(NA, nrow=ncol(x), ncol=ncol(y)),
                       crand=matrix(NA, nrow=ncol(x), ncol=ncol(y)))
        for(k in 1:ncol(x)){
            for(l in 1:ncol(y)){
                tab <- table(x[,k], y[,l])
                m <- matchClasses(tab, method=method, iter=iter,
                                  verbose=verbose, maxexact=maxexact)
                a <- classAgreement(tab[,m])
                retval$diag[k,l] <- a$diag
                retval$kappa[k,l] <- a$kappa
                retval$rand[k,l] <- a$rand
                retval$crand[k,l] <- a$crand
            }
        }
    }
    
    retval
}


permutations <- function(n) { 
    z <- matrix(1)
    for (i in 2:n) { 
        x <- cbind(z, i)
        a <- c(1:i, 1:(i - 1))
        z <- matrix(0, ncol=ncol(x), nrow=i*nrow(x))
        z[1:nrow(x),] <- x 
        for (j in 2:i-1) { 
            z[j*nrow(x)+1:nrow(x),] <- x[, a[1:i+j]] 
        } 
    } 
    dimnames(z) <- NULL
    z
} 
moment <- function(x, order = 1, center = FALSE, absolute = FALSE,
		   na.rm = FALSE) {
  if (na.rm) 
    x <- x[!is.na(x)]
  if (center)
    x <- x - mean(x)
  if (absolute)
    x <- abs(x)
  sum(x ^ order) / length(x)
}
rmvnorm <- function(n, mu=rep(0, nrow(sigma)),
                      sigma=diag(length(mu))){

  if(nrow(sigma) != ncol(sigma)){
    stop("sigma must be a square matrix")
  }

  if(length(mu) != nrow(sigma)){
    stop("mu and sigma have non-conforming size")
  }
  
  sigsvd <- svd(sigma)
  retval <- t(sigsvd$v %*% (t(sigsvd$u) * sqrt(sigsvd$d)))
  retval <- matrix(rnorm(n * ncol(sigma)), nrow = n) %*% retval
  retval <- sweep(retval, 2, mu, "+")
  retval
}


dmvnorm <- function(x, mu, sigma){

  if(is.vector(x)){
    x <- matrix(x, ncol=length(x))
  }

  if(missing(mu)){
    mu <- rep(0, length=ncol(x))
  }
  
  if(missing(sigma)){
    sigma <- diag(ncol(x))
  }

  if(ncol(x) != ncol(sigma)){
    stop("x and sigma have non-conforming size")
  }
  
  if(nrow(sigma) != ncol(sigma)){
    stop("sigma must be a square matrix")
  }
  if(length(mu) != nrow(sigma)){
    stop("mu and sigma have non-conforming size")
  }

  retval <- exp(-mahalanobis(x, center=mu, cov=sigma)/2)
  det <- prod(eigen(sigma, sym=TRUE)$values)
  retval<- retval / (sqrt(det) * sqrt(2*pi)^ncol(x))

  retval
}
  
plot.fclust <- function(clobj, x, centers=TRUE, initcenters=FALSE,
                         color=rainbow(clobj$learning$ncenters),...){
  
  x <- as.matrix(x)

  
  if(dim(x)[2]>2){
    pairs(x, col=color[cl$cluster], ...)
  }
  else{
    plot(x, col=color[cl$cluster], ...)
    if(centers)
      points(cl$centers, pch=4,col=color,cex=2)
    if(initcenters)
      points(clobj$learning$initcenters, pch=3,col=color,cex=2)
  }
}

plot.stft <- function (Y, col = gray (63:0/63))
  {
    x <- Y$values
    image(x=1:dim(x)[1], y=1:dim(x)[2], z=x, col=col)
}
# Generated automatically from pnm.R.in by configure.

plot.pnm <- function(pnmobj, xlab="", ylab="",
                     axes=FALSE, ...) {

  d <- dim(pnmobj)
  maxval <- attr(pnmobj, "maxval")

  if(attr(pnmobj, "type") =="ppm"){
    col <- rgb(pnmobj[1,,]/maxval,pnmobj[2,,]/maxval,pnmobj[3,,]/maxval)
    z <- matrix(1:length(col), nrow=d[2])
    image(x=1:d[3], y=1:d[2], z=t(z[d[2]:1,]),
        col=col, xlab=xlab, ylab=ylab, axes=axes, ...)
  }
  else{
    image(x=1:d[2], y=1:d[1], z=t(pnmobj[d[1]:1,]),
          col=gray((0:maxval)/maxval),
          xlab=xlab, ylab=ylab, axes=axes, ...)
  }
}


read.pnm <- function(file){

    if(! TRUE) stop("Sorry, no pnm library available")


    pnmhead <- .C("readpnminit",
                  file = as.character(file),
                  nc = as.integer(1),
                  nr = as.integer(1),
                  maxval = as.integer(1),
                  type="XXX")
    
    red <- integer(pnmhead$nc * pnmhead$nr)
    if(pnmhead$type == "ppm"){
        green <- integer(pnmhead$nc * pnmhead$nr)
        blue <- integer(pnmhead$nc * pnmhead$nr)
        retval <- array(dim=c(3, pnmhead$nr, pnmhead$nc),
                        dimnames=c(c("red", "gtreen", "blue"), NULL, NULL))
    }
    else {
        green <- as.integer(1)
        blue <- as.integer(1)
        retval <- matrix(ncol = pnmhead$nc, nrow=pnmhead$nr)
    }
    
    
    pnm <- .C("readpnm",
              file = as.character(file),
              red=red, green=green, blue=blue);
    
    
    if(pnmhead$type == "ppm"){
        retval[1,,] <- matrix(pnm$red, ncol = pnmhead$nc, byrow=TRUE)
        retval[2,,] <- matrix(pnm$green, ncol = pnmhead$nc, byrow=TRUE)
        retval[3,,] <- matrix(pnm$blue, ncol = pnmhead$nc, byrow=TRUE)
    }
    else{
        retval <- matrix(pnm$red, ncol = pnmhead$nc, byrow=TRUE)
    }
    attr(retval, "maxval") <- pnmhead$maxval
    attr(retval, "type") <- pnmhead$type
    class(retval) <- "pnm"
    retval
}


write.pgm <- function(pgmobj, file="Rimage.pgm",
                      forceplain=FALSE){

    if(! TRUE) stop("Sorry, no pnm library available")

    retval <- .C("writepgm",
                 file = as.character(file),
                 image = as.integer(t(pgmobj)),
                 nc = as.integer(ncol(pgmobj)),
                 nr = as.integer(nrow(pgmobj)),
                 maxval = as.integer(attr(pgmobj, "maxval")),
                 forceplain = as.integer(forceplain))
}


channel.pnm <- function(pnmobj, chan="red"){
    
    chan <- pmatch(chan, c("red", "green", "blue"))
    
    retval <- pnmobj[chan,,]
    attr(retval, "maxval") <- attr(pnmobj, "maxval")
    attr(retval, "type") <- "pgm"
    class(retval) <- "pnm"
    retval
}
  
print.fclust <- function (clobj)
  {
    
    if (!is.null(clobj$learning$iter))
      cat("\n                            Clustering on Training Set\n\n\n")
    else
      cat("\n                              Clustering on Test Set\n\n\n")
    
    cat("Number of Clusters: ", clobj$learning$ncenters, "\n")
    cat("Sizes  of Clusters: ", clobj$size, "\n\n")
    
  }

rbridge <- function(end=1, frequency=1000) {

  z <- rwiener(end=end, frequency=frequency)
  ts(z - time(z)*as.vector(z)[frequency],
     start=1/frequency, frequency=frequency)
}

read.octave <- function (file, quiet=FALSE) {

  nr <- 0
  nc <- 0

  if(!quiet)
    cat("Header: ")
  
  head <- scan(file=file,what=character(),nlines=4, sep=":", quiet=quiet)
  if(length(head) != 8){
    stop("Header seem to be corrupt")
  }
  for(k in 1:4){
    if(head[2*k-1] == "# rows"){
      nr <- as.integer(head[2*k])
    }
      else if(head[2*k-1] == "# columns"){
	nc <- as.integer(head[2*k])
      }
  }

  if(!quiet)
    cat("Data  : ")

  z <- scan(file=file,skip=4,quiet=quiet)
  if(length(z) != nc*nr){
    stop("Wrong number of data elements")
  }

  if((nr>1) && (nc>1)){
    if(!quiet)
      cat(paste("Matrix:", nr, "rows,", nc, "columns\n"))
    
    z<-matrix(z, nrow=nr, ncol=nc, byrow=TRUE)
  }
    else if(!quiet){
      cat("Vector:", nr*nc, "elements\n")
    }
  z
}
	      
rectangle.window <- function (n)
  rep (1, n)
rwiener <- function(end=1, frequency=1000) {

  z<-cumsum(rnorm(end*frequency)/sqrt(frequency))
  ts(z, start=1/frequency, frequency=frequency)
}

scaclust <- function (x, centers, iter.max = 100, verbose = FALSE,
                    method = "ad", theta= NULL) 
{
  xrows <- dim(x)[1]
  xcols <- dim(x)[2]
  xold <- x
  perm <- sample(xrows)
  x <- x[perm, ]
  ## initial values are given
  if (is.matrix(centers)) 
    ncenters <- dim(centers)[1]
  else {
   ## take centers random vectors as initial values
    ncenters <- centers
    centers <- x[rank(runif(xrows))[1:ncenters], ]+0.001
  }
##  method <- pmatch(method, c("ad", "mtv", "sand","ml"))
  method <- pmatch(method, c("ad", "mtv", "sand","mlm"))
  if (is.na(method)) 
    stop("invalid clustering method")
  if (method == -1) 
    stop("ambiguous clustering method")
  
  if (method == 1) {
    beta <- 1/xcols
    taf <- 0 }
  if (method == 2) {
    beta <- 0.5
    taf <- xcols/2 }
  if (method == 3) {
    beta <- 1/xcols
    taf <- 1 }
  if (method == 4){
    beta <- 0.0
    taf <- -1 }

  
  ##initialize theta
 ## if (method != 4){
  if (missing(theta))
    theta <- rep(1.0,ncenters)
  else
    theta <- as.double(theta)
  
  ##}
  
  
  initcenters <- centers
  ## dist <- matrix(0, xrows, ncenters)
  ## necessary for empty clusters
  pos <- as.factor(1:ncenters)
  rownames(centers) <- pos
  iter <- integer(1)
  
##  if ((method == 1) || (method == 2) || (method == 3)){
    retval <- .C("common", xrows = as.integer(xrows),
                 xcols = as.integer(xcols), 
                 x = as.double(x), ncenters = as.integer(ncenters), 
                 centers = as.double(centers), 
                 itermax = as.integer(iter.max), iter = as.integer(iter), 
                 verbose = as.integer(verbose), U=double(xrows*ncenters),
                 UANT=double(xrows*ncenters), beta=as.double(beta),
                 taf=as.double(taf), theta=as.double(theta),ermin=double(1))

  centers <- matrix(retval$centers, ncol = xcols, dimnames = dimnames(initcenters))
  
  U <- retval$U
  U <- matrix(U, ncol=ncenters)
  ##  clusterU <- max.col(U)
  clusterU <- apply(U,1,which.max)
  clusterU <- clusterU[order(perm)]
  U <- U[order(perm),]
  
  clustersize <- as.integer(table(clusterU))
  
  retval <- list(centers = centers, cluster = clusterU,
                 size = clustersize, member=U, error = retval$ermin,
                 learning = list(ncenters = ncenters,
                   initcenters = initcenters, iter = retval$iter - 1,
                   theta = theta), call = match.call())
  
  class(retval) <- c("fclust")
  return(retval)
}



sigmoid <- function(x) 1/(1 + exp(-x))

dsigmoid <- function(x) sigmoid(x) * (1 - sigmoid(x))

d2sigmoid <- function(x) dsigmoid(x) * (1 - 2 * sigmoid(x))
skewness <- function (x, na.rm = FALSE)
{
  if (na.rm) 
    x <- x[!is.na(x)]
  sum((x-mean(x))^3)/(length(x)*sd(x)^3)
}

stft <- function(X, win=min(80,floor(length(X)/10)), 
                 inc=min(24, floor(length(X)/30)), coef=64, 
		 wtype="hanning.window")
  {
    numcoef <- 2*coef
    if (win > numcoef)
      {
	win <- numcoef
	cat ("stft: window size adjusted to", win, ".\n")
      }
    numwin <- trunc ((length(X) - win) / inc)

    ## compute the windows coefficients
    wincoef <- eval(parse(text=wtype))(win)

    ## create a matrix Z whose columns contain the windowed time-slices
    z <- matrix (0, numwin + 1, numcoef)
    y <- z
    st <- 1
    for (i in 0:numwin)
      {
	z[i+1, 1:win] <- X[st:(st+win-1)] * wincoef
	y[i+1,] <- fft(z[i+1,])
	st <- st + inc
      }

    Y<- list (values = Mod(y[,1:coef]), windowsize=win, increment=inc,
		  windowtype=wtype)
    class(Y) <- "stft"
    return(Y)
  }



svm <-
function (x,
          y,
          svm.type    = NULL,
          kernel.type = "radial",
          degree      = 3,
          gamma       = 1/ncol(as.matrix(x)),
          coef0       = 0,
          cost        = 1,
          nu          = 0.5,
          cachesize   = 40,
          tolerance   = 0.001,
          epsilon     = 0.5,
          shrinking   = TRUE)
{
  if (is.vector(x))
    x <- t(t(x))
  else
    x <- as.matrix(x)

  if (is.null (svm.type)) svm.type <-
      if (is.factor(y)) "C-classification" else "regression"

  svm.type    <- pmatch (svm.type, c("C-classification",
                                     "nu-classification",
                                     "one-classification",
                                     "regression"),1) - 1
  
  kernel.type <- pmatch (kernel.type, c("linear",
                                        "polynomial",
                                        "radial",
                                        "sigmoid"),3) - 1

  if (!is.vector(y) && !is.factor (y)) stop ("y must be a vector or a factor.")
  if (length(y) != nrow(x)) stop ("x and y don't match.")

  if (cachesize < 0.1) cachesize <- 0.1
  
  lev <- NULL
  # in case of classification: map levels into {-1,1}
  if (is.factor(y)) {
    lev <- levels (y)
    y <- codes (y) * 2 - 3
  } else if (svm.type < 3) {
    lev <- levels (as.ordered (y))
    y <- codes (as.ordered(y)) * 2 - 3
  }

  if (length (lev) > 2) stop ("sorry, can't handle more than 2 classes !")
    
  cret <- .C ("svmtrain",
              as.double  (t(x)), as.integer(nrow (x)), as.integer(ncol(x)),
              as.double  (y),
              as.integer (svm.type),
              as.integer (kernel.type),
              as.double  (degree),
              as.double  (gamma),
              as.double  (coef0),
              as.double  (cost),
              as.double  (nu),
              as.double  (cachesize),
              as.double  (tolerance),
              as.double  (epsilon),
              as.integer (shrinking),
              nr    = integer (1),
              index = integer (nrow(x)),
              coefs = double  (nrow(x)),
              rho   = double  (1)
             )
  
  ret <- list (
               call        = match.call(),
               svm.type    = svm.type,
               kernel.type = kernel.type,
               cost        = cost,
               degree      = degree,
               gamma       = gamma,
               coef0       = coef0,
               nu          = nu,
               epsilon     = epsilon,
               
               levels      = lev,
               nr          = cret$nr,                  #number of sv
               sv          = t(t(x[cret$index==1,])),        #copy of sv
               index       = which (cret$index==1),    #indexes of sv in x
               rho         = cret$rho,                 #constant in decision function
               coefs       = cret$coefs[cret$index==1] #coefficiants of sv
              )
  class (ret) <- "svm"
  ret
} 

predict.svm <- function (model, x, type = "class") {

  type <- pmatch (type, c("class","raw"), 1)

  if (is.vector (x))
    x <- t(t(x))
  else
    x <- as.matrix(x)
    
  if (ncol(x) != ncol(model$sv)) stop ("test vector does not match model !")

  ret <- .C ("svmclassify",
             #model
             as.double  (t(model$sv)),
             as.integer (nrow(model$sv)), as.integer(ncol(model$sv)),
             as.double  (model$coefs),
             as.double  (model$rho),
             
             #parameter
             as.integer (model$svm.type),
             as.integer (model$kernel.type),
             as.double  (model$degree),
             as.double  (model$gamma),
             as.double  (model$coef0),

             #test matrix
             as.double (t(x)), as.integer (nrow(x)),
             
             #decision-values
             ret = double  (nrow(x))
            )$ret

  if ((type == 2) || is.null(model$levels))
    ret
  else {
    ret2 <- rep (model$levels[2],nrow(x))
    ret2 [which(ret < 0)] <- model$levels[1]
    as.factor (ret2)
  }
}

print.svm <- function (model) {
  cat ("\nCall:\n",deparse (model$call),"\n\n")
  cat ("Parameters:\n")
  cat ("   SVM-Type: ",c("C-classification",
                         "nu-classification",
                         "one-classification",
                         "regression")[model$svm.type+1],"\n")
  cat (" SVM-Kernel: ",c("linear",
                         "polynomial",
                         "radial",
                         "sigmoid")[model$kernel.type+1],"\n")
  cat ("       cost: ",model$cost,"\n")
  cat ("     degree: ",model$degree,"\n")
  cat ("      gamma: ",model$gamma,"\n")
  cat ("     coef.0: ",model$coef0,"\n")
  cat ("         nu: ",model$nu,"\n")
  cat ("    epsilon: ",model$epsilon,"\n")
  cat ("       cost: ",model$cost,"\n\n")
  
  cat ("\nNumber of Support Vectors:\n",model$nr,"\n\n")
  cat ("Rho:\n",model$rho,"\n\n")
  
}

summary.svm <- function (x) {
  print (x)
  cat ("Support Vectors:\n")
  print (x$sv)
  cat ("\n\nCoefficiants:\n")
  print (x$coefs)
  cat ("\n\n")
}
  
.First.lib <- function(lib, pkg){

    library.dynam("e1071", pkg, lib)
}
